Extend config file to generate tiled LCZ and sprawl areas

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

@@ -578,6 +578,18 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils {
                         if (grid_rowCol && grid_rowCol in Boolean) {
                             grid_indicators_tmp.rowCol = grid_rowCol
                         }
+                        def lcz_lod = grid_indicators.lcz_lod
+                        if (lcz_lod && lcz_lod in Integer) {
+                            if (lcz_lod < 0 && lcz_lod >10) {
+                                error "The number of level of details to aggregate the LCZ must be between 0 and 10"
+                                return
+                            }
+                            grid_indicators_tmp.put("lcz_lod", lcz_lod)
+                        }
+                        def sprawl_areas = grid_indicators.sprawl_areas
+                        if (sprawl_areas && sprawl_areas in Boolean) {
+                            grid_indicators_tmp.put("sprawl_areas", sprawl_areas)
+                        }
 
                         defaultParameters.put("grid_indicators", grid_indicators_tmp)
                     } else {
@@ -634,7 +646,8 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils {
                 "grid_indicators",
                 "road_traffic",
                 "population",
-                "ground_acoustic"]
+                "ground_acoustic",
+                "sprawl_areas"]
     }
 
 
@@ -681,8 +694,9 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils {
         }
 
         def tmp_results = [:]
+        def rows =  h2gis_datasource.rows("SELECT ST_CollectionExtract(THE_GEOM, 3) as the_geom,code_insee  FROM COMMUNE")
         //We process each zones because the input zone can overlap several communes
-        h2gis_datasource.eachRow("SELECT ST_CollectionExtract(THE_GEOM, 3) as the_geom,code_insee  FROM COMMUNE") { row ->
+        rows.each { row ->
             Geometry geom = row.the_geom
             def code_insee = row.code_insee
             info "Processing the commune with the code insee :  ${code_insee}"
@@ -835,20 +849,28 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils {
         def grid_indicators_params = processing_parameters.grid_indicators
         //Compute the grid indicators
         if (grid_indicators_params) {
+            info("Start computing grid_indicators")
             def x_size = grid_indicators_params.x_size
             def y_size = grid_indicators_params.y_size
             Geometry geomEnv = h2gis_datasource.getSpatialTable(results.zone).getExtent()
             String gridTableName = Geoindicators.WorkflowGeoIndicators.createGrid(h2gis_datasource, geomEnv,
                     x_size, y_size, srid, grid_indicators_params.rowCol)
             if (gridTableName) {
                 String rasterizedIndicators = Geoindicators.WorkflowGeoIndicators.rasterizeIndicators(h2gis_datasource, gridTableName,
-                        grid_indicators_params.indicators,
+                        grid_indicators_params.indicators,grid_indicators_params.lcz_lod,
                         results.building, results.road, results.vegetation,
                         results.water, results.impervious,
                         results.rsu_lcz, results.rsu_utrf_area, "", "",
                         processing_parameters.prefixName)
                 if (rasterizedIndicators) {
                     results.put("grid_indicators", rasterizedIndicators)
+                    if(grid_indicators_params.sprawl_areas){
+                       String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2gis_datasource, rasterizedIndicators, Math.max(x_size,y_size))
+                        if(sprawl_areas){
+                            results.put("sprawl_areas", sprawl_areas)
+                        }
+                    }
+                    info("End computing grid_indicators")
                 }
             } else {
                 info "Cannot create a grid to aggregate the indicators"

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

@@ -248,9 +248,10 @@ Map formatLayers(JdbcDataSource datasource, Map layers, float distance, float hL
         hLevMin = 3
     }
     if (!datasource) {
-        error "The database to store the BD Topo data doesn't exist"
+        error "The database to store the BDTopo data doesn't exist"
         return
     }
+    info "Formating BDTopo GIS layers"
     //Prepare the existing bdtopo data in the local database
     def importPreprocess = BDTopo.InputDataLoading.loadV2(datasource, layers, distance)
 
@@ -292,6 +293,8 @@ Map formatLayers(JdbcDataSource datasource, Map layers, float distance, float hL
 
     debug "End of the BDTopo extract transform process."
 
+    info "All layers have been formatted"
+
     return ["building"  : finalBuildings, "road": finalRoads, "rail": finalRails, "water": finalHydro,
             "vegetation": finalVeget, "impervious": finalImpervious, "urban_areas": urbanAreas, "zone": zoneTable]
 

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

@@ -447,6 +447,8 @@ Map formatLayers(JdbcDataSource datasource, Map layers, float distance, float hL
         error "The database to store the BD Topo data doesn't exist"
         return
     }
+    info "Formating BDTopo GIS layers"
+
     //Prepare the existing bdtopo data in the local database
     def importPreprocess = BDTopo.InputDataLoading.loadV3(datasource, layers, distance)
 
@@ -487,7 +489,7 @@ Map formatLayers(JdbcDataSource datasource, Map layers, float distance, float hL
             importPreprocess.water,
             zoneTable)
 
-    debug "End of the BDTopo extract transform process."
+    info "All layers have been formatted"
 
     return ["building"  : finalBuildings, "road": finalRoads, "rail": finalRails, "water": finalHydro,
             "vegetation": finalVeget, "impervious": finalImpervious, "urban_areas": urbanAreas, "zone": zoneTable]

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

@@ -333,7 +333,7 @@ String formatBuildingLayer(JdbcDataSource datasource, String building, String zo
             }
         }
     }
-    debug 'Buildings transformation finishes'
+    info "Building formatted"
     return outputTableName
 }
 
@@ -567,7 +567,7 @@ String formatRoadLayer(JdbcDataSource datasource, String road, String zone = "")
             }
         }
     }
-    debug('Roads transformation finishes')
+    info "Road formatted"
     return outputTableName
 }
 
@@ -648,7 +648,7 @@ String formatHydroLayer(JdbcDataSource datasource, String water, String zone = "
             }
         }
     }
-    debug('Hydro transformation finishes')
+    info "Water formatted"
     return outputTableName
 }
 
@@ -744,7 +744,7 @@ String formatRailsLayer(JdbcDataSource datasource, String rail, String zone = ""
             }
         }
     }
-    debug('Rails transformation finishes')
+    info "Rail formatted"
     return outputTableName
 }
 
@@ -851,7 +851,7 @@ String formatVegetationLayer(JdbcDataSource datasource, String vegetation, Strin
             }
         }
     }
-    debug('Vegetation transformation finishes')
+    info "Vegetation formatted"
     return outputTableName
 }
 
@@ -908,7 +908,7 @@ String formatImperviousLayer(H2GIS datasource, String impervious) {
         }
     }
     datasource.execute("DROP TABLE IF EXISTS $polygonizedTable".toString())
-    debug('Impervious areas transformation finishes')
+    info "Impervious areas formatted"
     return outputTableName
 }
 

bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowDebugTest.groovy

@@ -137,7 +137,7 @@ class WorkflowDebugTest {
     @Test
     void testIntegrationFolderInput() {
         def input_data = "/media/ebocher/Extreme SSD/bdtopo/bdtopo2/BDTOPO_2-2_TOUSTHEMES_SHP_LAMB93_D035_2018-09-25/BDTOPO/1_DONNEES_LIVRAISON_2018-11-00144/BDT_2-2_SHP_LAMB93_D035-ED182"
-        def locations = [ "Redon"]
+        def locations = [ "Bordeaux"]
         String directory = "/tmp/bdtopo2"
         File dirFile = new File(directory)
         dirFile.delete()
@@ -146,7 +146,7 @@ class WorkflowDebugTest {
                 "description" : "Example of configuration file to run the BDTopo workflow and store the results in a folder",
                 "geoclimatedb": [
                         "folder": dirFile.absolutePath,
-                        "name"  : "bdtopo_workflow_db;AUTO_SERVER=TRUE",
+                        "name"  : "bdtopo_workflow_db",
                         "delete": true
                 ],
                 "input"       : [
@@ -158,25 +158,31 @@ class WorkflowDebugTest {
                 "parameters"  :
                         ["distance"       : 0,
                          rsu_indicators   : [
-                                 "indicatorUse": ["LCZ",  "UTRF"]
+                                 "indicatorUse": ["LCZ",  "UTRF",  "TEB"]
+
                          ],
                          "grid_indicators": [
-                                 "x_size"    : 1000,
-                                 "y_size"    : 1000,
+                                 "x_size"    : 100,
+                                 "y_size"    : 100,
                                  "indicators" :["BUILDING_FRACTION", "BUILDING_HEIGHT", "BUILDING_POP",
                                  "BUILDING_TYPE_FRACTION", "WATER_FRACTION", "VEGETATION_FRACTION",
-                                 "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "FREE_EXTERNAL_FACADE_DENSITY",
+                                 "ROAD_FRACTION", "IMPERVIOUS_FRACTION",
+                                                //"FREE_EXTERNAL_FACADE_DENSITY",
                                  "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY",
-                                 "SEA_LAND_FRACTION", "ASPECT_RATIO", "SVF",
+                                 "SEA_LAND_FRACTION", "ASPECT_RATIO",
+                                                //"SVF",
                                  "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS",
                                  "UTRF_AREA_FRACTION", "UTRF_FLOOR_AREA_FRACTION",
                                  "LCZ_PRIMARY"]
+                                 ,
+                                 "lcz_lod":2,
+                                 "sprawl_areas":true
                          ]
                         ]
         ]
         //BDTopo.v2(bdTopoParameters)
 
-        input_data = "/media/ebocher/Extreme SSD/bdtopo/bdtopo3/BDTOPO_3-0_TOUSTHEMES_SHP_LAMB93_D035_2022-09-15/BDTOPO"
+        input_data = "/home/ebocher/Téléchargements/BDTOPO_3-3_TOUSTHEMES_SHP_LAMB93_D033_2023-12-15/BDTOPO/1_DONNEES_LIVRAISON_2023-12-00099"
 
         directory = "/tmp/bdtopo3"
         dirFile = new File(directory)
@@ -189,4 +195,5 @@ class WorkflowDebugTest {
 
     }
 
+
 }

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

@@ -1002,12 +1002,11 @@ String linearRoadOperations(JdbcDataSource datasource, String rsuTable, String r
                     levelConsiderated.each { filtering += "$baseFiltering AND b.$Z_INDEX=$it OR " }
                     filtering = filtering[0..-4]
                 }
-                def selectionQuery = "DROP TABLE IF EXISTS $roadInter; " +
-                        "CREATE TABLE $roadInter AS SELECT a.$ID_COLUMN_RSU AS id_rsu, " +
-                        "ST_AREA(a.$GEOMETRIC_COLUMN_RSU) AS rsu_area, ST_INTERSECTION(a.$GEOMETRIC_COLUMN_RSU, " +
-                        "b.$GEOMETRIC_COLUMN_ROAD) AS the_geom $ifZindex FROM $rsuTable a, $roadTable b " +
-                        "WHERE $filtering;"
-                datasource selectionQuery.toString()
+                datasource.execute("""DROP TABLE IF EXISTS $roadInter; 
+                        CREATE TABLE $roadInter AS SELECT a.$ID_COLUMN_RSU AS id_rsu, 
+                        ST_AREA(a.$GEOMETRIC_COLUMN_RSU) AS rsu_area, ST_INTERSECTION(a.$GEOMETRIC_COLUMN_RSU, 
+                        b.$GEOMETRIC_COLUMN_ROAD) AS the_geom $ifZindex FROM $rsuTable a, $roadTable b 
+                        WHERE $filtering;""".toString())
 
                 // If all roads are considered at the same level...
                 if (!levelConsiderated) {

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

@@ -910,7 +910,7 @@ Map computeTypologyIndicators(JdbcDataSource datasource, String building_indicat
     def utrfFloorArea = "UTRF_RSU_FLOOR_AREA" + baseNameUtrfRsu
 
     if (indicatorUse.contains("LCZ")) {
-        info """ The LCZ classification will be performed """
+        info """The LCZ classification will be performed """
         def lczIndicNames = ["GEOM_AVG_HEIGHT_ROOF"              : "HEIGHT_OF_ROUGHNESS_ELEMENTS",
                              "BUILDING_FRACTION_LCZ"             : "BUILDING_SURFACE_FRACTION",
                              "ASPECT_RATIO"                      : "ASPECT_RATIO",
@@ -948,7 +948,7 @@ Map computeTypologyIndicators(JdbcDataSource datasource, String building_indicat
     // If the UTRF indicators should be calculated, we only affect a URBAN typo class
     // to each building and then to each RSU
     if (indicatorUse.contains("UTRF") && runUTRFTypology) {
-        info """ The URBAN TYPOLOGY classification is performed """
+        info """The URBAN TYPOLOGY classification is performed """
         def gatheredScales = Geoindicators.GenericIndicators.gatherScales(datasource,
                 building_indicators, block_indicators,
                 rsu_indicators, "BUILDING",
@@ -1104,7 +1104,7 @@ Map createUnitsOfAnalysis(JdbcDataSource datasource, String zone, String buildin
                           String rsu, double surface_vegetation,
                           double surface_hydro, double surface_urban_areas,
                           double snappingTolerance, List indicatorUse = ["LCZ", "UTRF", "TEB"], String prefixName = "") {
-    info "Create the units of analysis..."
+    info "Create the spatial units..."
     def idRsu = "id_rsu"
     def tablesToDrop = []
     if (!rsu) {
@@ -1744,6 +1744,7 @@ Map computeGeoclimateIndicators(JdbcDataSource datasource, String zone, String b
  * @param h2gis_datasource the local H2GIS database
  * @param gridTableName the name of the grid table to aggregate the data
  * @param list_indicators indicators names to compute
+ * @param lcz_lod level to aggregate the LCZ on the grid. Set by default to -1
  * @param buildingTable name
  * @param roadTable name
  * @param vegetationTable name
@@ -1757,7 +1758,7 @@ Map computeGeoclimateIndicators(JdbcDataSource datasource, String zone, String b
  * @return
  */
 String rasterizeIndicators(JdbcDataSource datasource,
-                           String grid, List list_indicators,
+                           String grid, List list_indicators, Integer lcz_lod,
                            String building, String road, String vegetation,
                            String water, String impervious, String rsu_lcz,
                            String rsu_utrf_area, String rsu_utrf_floor_area, String sea_land_mask,
@@ -1795,7 +1796,6 @@ String rasterizeIndicators(JdbcDataSource datasource,
                 false, "lcz")
         if (upperScaleAreaStatistics) {
             indicatorTablesToJoin.put(upperScaleAreaStatistics, grid_column_identifier)
-
             if (list_indicators_upper.contains("LCZ_PRIMARY")) {
                 def resultsDistrib = Geoindicators.GenericIndicators.distributionCharacterization(datasource,
                         upperScaleAreaStatistics, upperScaleAreaStatistics, grid_column_identifier,
@@ -2227,6 +2227,10 @@ String rasterizeIndicators(JdbcDataSource datasource,
 
     // Remove temporary tables
     datasource.dropTable(tablesToDrop)
+
+    if(lcz_lod){
+        return Geoindicators.GridIndicators.multiscaleLCZGrid(datasource, grid_indicators_table,grid_column_identifier,lcz_lod)
+    }
     return grid_indicators_table
 }
 

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

@@ -183,8 +183,6 @@ class WorkflowGeoIndicatorsTest {
                 inputTableNames.railTable, inputTableNames.vegetationTable,
                 inputTableNames.hydrographicTable, "", "", "", "", "",
                 ["indicatorUse": indicatorUse, svfSimplified: false], prefixName)
-        datasource.save(inputTableNames.vegetationTable, "/tmp/vegetation.fgb", true)
-        datasource.save(inputTableNames.roadTable, "/tmp/road.fgb", true)
         assertNotNull(geoIndicatorsCompute_i)
         checkRSUIndicators(datasource, geoIndicatorsCompute_i.rsu_indicators)
         assertEquals(listUrbTyp.Bu.sort(), datasource.getTable(geoIndicatorsCompute_i.building_indicators).columns.sort())
@@ -515,7 +513,7 @@ class WorkflowGeoIndicatorsTest {
         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, "building", null, null, null, null, null,
                 null, null, null)
         assertNull(grid_indicators)
         def list_indicators = ["BUILDING_FRACTION", "BUILDING_HEIGHT", "BUILDING_POP",
@@ -524,8 +522,8 @@ class WorkflowGeoIndicatorsTest {
                                "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,
+        grid_indicators = Geoindicators.WorkflowGeoIndicators.rasterizeIndicators(datasource, grid, list_indicators, null,
+                "building", null, null, null, null, null, null,
                 null, null, null)
         assertNotNull(grid_indicators)
         assertEquals(1, datasource.getRowCount(grid_indicators))
@@ -609,4 +607,19 @@ class WorkflowGeoIndicatorsTest {
         assertNotNull Geoindicators.buildNumber()
     }
 
+    @Disabled
+    @Test
+     void test(){
+        datasource.load("/tmp/road_inter.geojson", "road", true)
+        datasource.load("/tmp/rsu_table.geojson", "rsu", true)
+
+        datasource.createSpatialIndex("road")
+        datasource.createSpatialIndex("rsu")
+        datasource.execute("""
+        select st_intersection(a.the_geom, b.the_geom) as the_geom from road as a, rsu as b where a.the_geom && b.the_geom 
+        and st_intersects(a.the_geom, b.the_geom);
+        """.toString())
+    }
+
+
 }