Add new exact indicators at grid scale

bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy

@@ -542,7 +542,9 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils {
                     }
                     def allowed_grid_indicators = ["BUILDING_FRACTION", "BUILDING_HEIGHT", "BUILDING_POP", "BUILDING_TYPE_FRACTION", "WATER_FRACTION", "VEGETATION_FRACTION",
                                                    "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION", "UTRF_FLOOR_AREA_FRACTION", "LCZ_FRACTION", "LCZ_PRIMARY", "FREE_EXTERNAL_FACADE_DENSITY",
-                                                   "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION"]
+                                                   "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY",
+                                                   "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION", "ASPECT_RATIO",
+                                                   "SVF", "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"]
                     def allowedOutputIndicators = allowed_grid_indicators.intersect(list_indicators*.toUpperCase())
                     if (allowedOutputIndicators) {
                         //Update the RSU indicators list according the grid indicators

geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/DataUtils.groovy

@@ -73,7 +73,7 @@ String joinTables(JdbcDataSource datasource, Map inputTableNamesWithId, String o
             }
             leftQuery += " LEFT JOIN $key as $alias ON $alias.$value = $columnKey "
         }
-        indexes += "CREATE INDEX IF NOT EXISTS ${key}_ids ON $key USING BTREE($value);"
+        indexes += "CREATE INDEX IF NOT EXISTS ${key}_ids ON $key ($value);"
         alias++
     }
 

geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy

@@ -447,14 +447,19 @@ String distributionCharacterization(JdbcDataSource datasource, String distribTab
         def outputTableName = prefix prefixName, BASENAME
 
         // Get all columns from the distribution table and remove the geometry column if exists
-        def allColumns = datasource."$distribTableName".columns
+        def allColumns = datasource.getColumnNames(distribTableName)
         if (allColumns.contains(GEOMETRY_FIELD)) {
             allColumns -= GEOMETRY_FIELD
         }
         // Get the distribution columns and the number of columns
         def distribColumns = allColumns.minus(inputId.toUpperCase())
         def nbDistCol = distribColumns.size()
 
+        if(distribColumns.size()==0){
+            error("Any columns to compute the distribution from the table $distribTableName".toString())
+            return
+        }
+
         def idxExtrem = nbDistCol - 1
         def idxExtrem_1 = nbDistCol - 2
         if (extremum.toUpperCase() == "LEAST") {

geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy

@@ -126,7 +126,7 @@ String createTSU(JdbcDataSource datasource, String inputTableName, String inputz
                     CREATE TABLE $outputTableName AS 
                         SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(a.the_geom, $epsg) AS the_geom
                         FROM ST_EXPLODE('(
-                                SELECT ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))) AS the_geom 
+                                SELECT ST_BUFFER(ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))), -0.01) AS the_geom 
                                 FROM $inputTableName)') AS a,
                             $inputzone AS b
                         WHERE a.the_geom && b.the_geom 
@@ -136,9 +136,9 @@ String createTSU(JdbcDataSource datasource, String inputTableName, String inputz
         datasource """
                     DROP TABLE IF EXISTS $outputTableName;
                     CREATE TABLE $outputTableName AS 
-                        SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(ST_FORCE2D(the_geom), $epsg) AS the_geom 
+                        SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(st_buffer(the_geom, -0.01), $epsg) AS the_geom 
                         FROM ST_EXPLODE('(
-                                SELECT ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))) AS the_geom 
+                                SELECT ST_BUFFER(ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))),-0.01) AS the_geom 
                                 FROM $inputTableName)') where st_area(the_geom) > $area""".toString()
     }
 

geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy

@@ -104,6 +104,30 @@ class RsuIndicatorsTests {
         assertEquals 0d, h2GIS.firstRow("SELECT * FROM ${p} WHERE id_rsu = 5").FREE_EXTERNAL_FACADE_DENSITY
     }
 
+    @Test
+    void freeExternalFacadeDensityExactTest2() {
+        // Only the first 1 first created buildings are selected for the tests
+        h2GIS """
+                DROP TABLE IF EXISTS tempo_build, tempo_rsu; 
+                CREATE TABLE tempo_build(id_build int, the_geom geometry, height_wall double);
+                INSERT INTO tempo_build VALUES (1, 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'::GEOMETRY, 10),
+                                               (2, 'POLYGON ((10 0, 20 0, 20 20, 10 20, 10 0))'::GEOMETRY, 10),
+                                                (3, 'POLYGON ((30 30, 50 30, 50 50, 30 50, 30 30))'::GEOMETRY, 10);
+                CREATE TABLE tempo_rsu(id_rsu int, the_geom geometry);
+                INSERT INTO tempo_rsu VALUES    (1, 'POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))'::GEOMETRY);
+        """
+        // First calculate the correlation table between buildings and rsu
+        def buildingTableRelation = Geoindicators.SpatialUnits.spatialJoin(h2GIS,
+                "tempo_build", "tempo_rsu", "id_rsu", null, "test")
+
+        assertNotNull(buildingTableRelation)
+        def p = Geoindicators.RsuIndicators.freeExternalFacadeDensityExact(h2GIS,
+                buildingTableRelation, "tempo_rsu",
+                "id_rsu", "test")
+        assertNotNull(p)
+        assertEquals 0.16, h2GIS.firstRow("SELECT * FROM ${p} WHERE id_rsu = 1").FREE_EXTERNAL_FACADE_DENSITY
+    }
+
     @Test
     void groundSkyViewFactorTest() {
         // Only the first 1 first created buildings are selected for the tests
@@ -114,7 +138,7 @@ class RsuIndicatorsTests {
         h2GIS "DROP TABLE IF EXISTS corr_tempo; CREATE TABLE corr_tempo AS SELECT a.*, b.the_geom, b.height_wall " +
                 "FROM rsu_build_corr a, tempo_build b WHERE a.id_build = b.id_build"
 
-        def p = Geoindicators.RsuIndicators.groundSkyViewFactor(h2GIS, "rsu_test", "corr_tempo",
+        def p = Geoindicators.RsuIndicators.groundSkyViewFactor(h2GIS, "rsu_test","id_rsu", "corr_tempo",
                 0.008, 100, 60, "test")
         assertNotNull(p)
         assertEquals 0.54, h2GIS.firstRow("SELECT * FROM test_rsu_ground_sky_view_factor " +
@@ -155,7 +179,7 @@ class RsuIndicatorsTests {
         def listLayersBottom = [0, 10, 20, 30, 40, 50]
         def numberOfDirection = 4
         def rangeDeg = 360 / numberOfDirection
-        def p = Geoindicators.RsuIndicators.projectedFacadeAreaDistribution(h2GIS, "tempo_build", "rsu_test", listLayersBottom,
+        def p = Geoindicators.RsuIndicators.projectedFacadeAreaDistribution(h2GIS, "tempo_build", "rsu_test","id_rsu", listLayersBottom,
                 numberOfDirection, "test")
         assertNotNull(p)
         def concat = ""
@@ -189,7 +213,7 @@ class RsuIndicatorsTests {
         def listLayersBottom = [0, 10, 20, 30, 40, 50]
         def numberOfDirection = 4
         def rangeDeg = 360 / numberOfDirection
-        def p = Geoindicators.RsuIndicators.projectedFacadeAreaDistribution(h2GIS, "tempo_build", "rsu_test", listLayersBottom,
+        def p = Geoindicators.RsuIndicators.projectedFacadeAreaDistribution(h2GIS, "tempo_build", "rsu_test", "id_rsu",listLayersBottom,
                 numberOfDirection, "test")
         assertNotNull(p)
         def concat = ""
@@ -275,7 +299,7 @@ class RsuIndicatorsTests {
         def listLayersBottom = [0, 10, 20, 30, 40, 50]
         def numberOfDirection = 4
         def pFacadeDistrib = Geoindicators.RsuIndicators.projectedFacadeAreaDistribution(h2GIS, "tempo_build",
-                "rsu_test", listLayersBottom,
+                "rsu_test","id_rsu", listLayersBottom,
                 numberOfDirection, "test")
         assertNotNull(pFacadeDistrib)
         def pGeomAvg = Geoindicators.GenericIndicators.unweightedOperationFromLowerScale(h2GIS, "tempo_build",
@@ -292,7 +316,7 @@ class RsuIndicatorsTests {
         h2GIS "CREATE TABLE rsu_table AS SELECT a.*, b.geom_avg_height_roof, b.the_geom " +
                 "FROM test_rsu_projected_facade_area_distribution a, test_unweighted_operation_from_lower_scale b " +
                 "WHERE a.id_rsu = b.id_rsu"
-        def p = Geoindicators.RsuIndicators.effectiveTerrainRoughnessLength(h2GIS, "rsu_table",
+        def p = Geoindicators.RsuIndicators.effectiveTerrainRoughnessLength(h2GIS, "rsu_table","id_rsu",
                 "projected_facade_area_distribution", "geom_avg_height_roof",
                 listLayersBottom, numberOfDirection, "test")
         assertNotNull(p)
@@ -344,7 +368,7 @@ class RsuIndicatorsTests {
         h2GIS "DROP TABLE IF EXISTS rsu_tempo; CREATE TABLE rsu_tempo AS SELECT *, CASEWHEN(id_rsu = 1, 2.3," +
                 "CASEWHEN(id_rsu = 2, 0.1, null)) AS effective_terrain_roughness_length FROM rsu_test"
 
-        def p = Geoindicators.RsuIndicators.effectiveTerrainRoughnessClass(h2GIS, "rsu_tempo", "effective_terrain_roughness_length",
+        def p = Geoindicators.RsuIndicators.effectiveTerrainRoughnessClass(h2GIS, "rsu_tempo","id_rsu", "effective_terrain_roughness_length",
                 "test")
         assertNotNull(p)
         def concat = ""
@@ -699,7 +723,7 @@ class RsuIndicatorsTests {
                 "tempo_rsu",
                 "id_rsu",
                 [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50],
-                12,
+                12,true,
                 "test")
         assertNotNull(p)
         assertEquals 0.00566, h2GIS.firstRow("SELECT * FROM ${p} WHERE id_rsu = 1").FRONTAL_AREA_INDEX_H0_5_D30_60, 0.00001

geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy

@@ -177,6 +177,7 @@ class SpatialUnitsTests {
         assert h2GIS.getTable(pPointOnSurface).getRowCount() == 2
     }
 
+
     @Test
     void prepareGeometriesForRSUWithFilterTest() {
         h2GIS.load(SpatialUnitsTests.class.class.getResource("road_test.geojson"), true)

geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy

@@ -428,6 +428,56 @@ class WorkflowGeoIndicatorsTest {
         }
     }
 
+    @Test
+    void rasterizeIndicators1() {
+        datasource.execute("""
+        DROP TABLE IF EXISTS building;
+        CREATE TABLE BUILDING (id_build int, id_block int, id_rsu int, zindex int, the_geom geometry, height_wall float, height_roof float, type varchar, pop double precision );
+        INSERT INTO BUILDING VALUES (1, 1, 1, 0, 'POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))'::GEOMETRY, 5, 10, 'office', 100);
+        """.toString())
+        String grid = Geoindicators.WorkflowGeoIndicators.createGrid(datasource, datasource.getExtent("building"), 10, 10, 0)
+        assertNotNull(grid)
+        String grid_indicators = Geoindicators.WorkflowGeoIndicators.rasterizeIndicators(datasource, grid, [],
+                "building", null, null, null, null, null,
+                null, null, null)
+        assertNull(grid_indicators)
+        def list_indicators = ["BUILDING_FRACTION", "BUILDING_HEIGHT", "BUILDING_POP",
+                               "BUILDING_TYPE_FRACTION", "WATER_FRACTION", "VEGETATION_FRACTION",
+                               "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "FREE_EXTERNAL_FACADE_DENSITY",
+                               "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY",
+                               "SEA_LAND_FRACTION", "ASPECT_RATIO", "SVF",
+                               "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"]
+        grid_indicators = Geoindicators.WorkflowGeoIndicators.rasterizeIndicators(datasource, grid, list_indicators,
+                "building", null, null, null, null, null,
+                null, null, null)
+        assertNotNull(grid_indicators)
+        assertEquals(1, datasource.getRowCount(grid_indicators))
+        def rows = datasource.firstRow("select * from $grid_indicators".toString())
+        println(rows)
+        assertEquals(1d, rows.BUILDING_FRACTION)
+        assertEquals(0d, rows.HIGH_VEGETATION_FRACTION)
+        assertEquals(0d, rows.IMPERVIOUS_FRACTION)
+        assertEquals(0d, rows.LOW_VEGETATION_FRACTION)
+        assertEquals(0d, rows.ROAD_FRACTION)
+        assertEquals(0d, rows.WATER_FRACTION)
+        assertEquals(0d, rows.UNDEFINED_FRACTION)
+        assertEquals(100d, rows.SUM_POP)
+        assertEquals(10d, rows.AVG_HEIGHT_ROOF)
+        assertEquals(0d, rows.STD_HEIGHT_ROOF)
+        assertTrue(10d-rows.GEOM_AVG_HEIGHT_ROOF< 0.0001)
+        assertEquals(10d, rows.AVG_HEIGHT_ROOF_AREA_WEIGHTED)
+        assertEquals(0d, rows.STD_HEIGHT_ROOF_AREA_WEIGHTED)
+        assertEquals(2d, rows.FREE_EXTERNAL_FACADE_DENSITY)
+        assertEquals(3d, rows.BUILDING_SURFACE_DENSITY)
+        assertTrue(1.5 -rows.EFFECTIVE_TERRAIN_ROUGHNESS_LENGTH<0.001)
+        assertEquals(8, rows.EFFECTIVE_TERRAIN_ROUGHNESS_CLASS)
+        assertEquals(2d, rows.ASPECT_RATIO)
+        assertTrue(0.5 -rows.SVF<0.1)
+        assertEquals(1, rows.TYPE_OFFICE)
+
+    }
+
+
     /**
      * Method to check the result for the RSU indicators table
      * Please add new checks here

osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy

@@ -445,7 +445,7 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d
 
             def extract = OSMTools.Loader.extract(query)
             if (extract) {
-                Geometry geomArea = h2gis_datasource.firstRow("select st_extent(the_geom) as the_geom from ${zones.zone}").the_geom
+                Geometry geomArea = h2gis_datasource.firstRow("select st_extent(the_geom) as the_geom from ${zones.zone}".toString()).the_geom
                 geomArea.setSRID(srid)
                 Map gisLayersResults = OSM.InputDataLoading.createGISLayers(h2gis_datasource, extract, geomArea, srid)
                 if (gisLayersResults) {
@@ -846,8 +846,11 @@ def extractProcessingParameters(def processing_parameters) {
                     return
                 }
                 def allowed_grid_indicators = ["BUILDING_FRACTION", "BUILDING_HEIGHT", "BUILDING_POP", "BUILDING_TYPE_FRACTION", "WATER_FRACTION", "VEGETATION_FRACTION",
-                                               "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION","UTRF_FLOOR_AREA_FRACTION",  "LCZ_FRACTION", "LCZ_PRIMARY", "FREE_EXTERNAL_FACADE_DENSITY",
-                                               "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION"]
+                                               "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION","UTRF_FLOOR_AREA_FRACTION",
+                                               "LCZ_FRACTION", "LCZ_PRIMARY", "FREE_EXTERNAL_FACADE_DENSITY",
+                                               "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", "BUILDING_HEIGHT_DIST",
+                                               "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION", "ASPECT_RATIO", "SVF",
+                                               "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"]
                 def allowedOutputIndicators = allowed_grid_indicators.intersect(list_indicators*.toUpperCase())
                 if (allowedOutputIndicators) {
                     //Update the RSU indicators list according the grid indicators

osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy

@@ -652,6 +652,7 @@ class WorflowOSMTest extends WorkflowAbstractTest {
         dirFile.delete()
         dirFile.mkdir()
 
+
        //def nominatim = OSMTools.Utilities.getNominatimData("Lorient")
 
         def osm_parmeters = [
@@ -662,7 +663,7 @@ class WorflowOSMTest extends WorkflowAbstractTest {
                         "delete": false
                 ],
                 "input"       : [
-                        "locations": ["Toulouse"],//["Pont-de-Veyle"],//[nominatim["bbox"]],//["Lorient"],
+                        "locations": ["Redon"],//["Pont-de-Veyle"],//[nominatim["bbox"]],//["Lorient"],
                          "area": 2800,
                         /*"timeout":182,
                         "maxsize": 536870918,
@@ -673,22 +674,21 @@ class WorflowOSMTest extends WorkflowAbstractTest {
                 "parameters"  :
                         ["distance"                                             : 0,
                          "rsu_indicators"                                       : [
-                                 "indicatorUse": ["LCZ", "UTRF", "TEB"]
+                                 "indicatorUse": ["LCZ"] //, "UTRF", "TEB"]
                          ],"grid_indicators": [
                                 "x_size": 100,
                                 "y_size": 100,
                                 //"rowCol": true,
                                 "indicators":  ["BUILDING_FRACTION","BUILDING_HEIGHT", "BUILDING_POP",
                                                 "BUILDING_TYPE_FRACTION","WATER_FRACTION","VEGETATION_FRACTION",
-                                                "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION",
-                                                "UTRF_FLOOR_AREA_FRACTION",
-                                                "LCZ_FRACTION", "LCZ_PRIMARY", "FREE_EXTERNAL_FACADE_DENSITY",
-                                                "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY",
-                                                "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION"]
-                        ],    "worldpop_indicators": true,
-                         "road_traffic"                                         : true,
+                                                "ROAD_FRACTION", "IMPERVIOUS_FRACTION",
+                                                "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY",  "SEA_LAND_FRACTION",
+                                                "ASPECT_RATIO","SVF",
+                                                "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"]
+                        ],    "worldpop_indicators": false,
+                         "road_traffic"                                         : false,
                          "noise_indicators"                                     : [
-                                 "ground_acoustic": true
+                                 "ground_acoustic": false
                          ]
                         ]
         ]

pom.xml

@@ -36,7 +36,7 @@
         <orbisdata-version>2.1.0-SNAPSHOT</orbisdata-version>
         <slf4j-version>2.0.9</slf4j-version>
         <logback-version>1.4.11</logback-version>
-        <groovy-version>3.0.11</groovy-version>
+        <groovy-version>3.0.19</groovy-version>
         <picocli-version>4.6.3</picocli-version>
         <commons-io-version>2.11.0</commons-io-version>
         <commons-compress-version>1.21</commons-compress-version>