diff --git a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy index aafc6d4fbe..f92bc8207c 100644 --- a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy +++ b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy @@ -550,7 +550,8 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { "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"] + "SVF", "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS","SPRAWL_AREAS", + "SPRAWL_DISTANCES", "SPRAWL_COOL_DISTANCE"] def allowedOutputIndicators = allowed_grid_indicators.intersect(list_indicators*.toUpperCase()) if (allowedOutputIndicators) { //Update the RSU indicators list according the grid indicators @@ -586,11 +587,6 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { } 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 { error "Please set a valid list of indicator names in ${allowed_grid_indicators}" @@ -863,12 +859,13 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { results.rsu_lcz, results.rsu_utrf_area, "", "", processing_parameters.prefixName) if (rasterizedIndicators) { + h2gis_datasource.dropTable(gridTableName) 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) - } + Map sprawl_indic = Geoindicators.WorkflowGeoIndicators.sprawlIndicators(h2gis_datasource,rasterizedIndicators, "id_grid", grid_indicators_params.indicators, + Math.max(x_size,y_size).floatValue()) + if(sprawl_indic){ + results.put("sprawl_areas", sprawl_indic.sprawl_areas) + results.put("grid_indicators", sprawl_indic.grid_indicators) } info("End computing grid_indicators") } diff --git a/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowDebugTest.groovy b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowDebugTest.groovy index c61bf1e282..bafbf0834b 100644 --- a/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowDebugTest.groovy +++ b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowDebugTest.groovy @@ -173,10 +173,9 @@ class WorkflowDebugTest { //"SVF", "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS", "UTRF_AREA_FRACTION", "UTRF_FLOOR_AREA_FRACTION", - "LCZ_PRIMARY"] - , - "lcz_lod":2, - "sprawl_areas":true + "LCZ_PRIMARY", "SPRAWL_AREAS", + "SPRAWL_DISTANCES", "SPRAWL_COOL_DISTANCE"] + //, "lcz_lod":2 ] ] ] diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy index 2ac325ad9d..1d25f68610 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy @@ -278,7 +278,7 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, Stri * @author Erwan Bocher (CNRS) * TODO : convert this method as a function table in H2GIS */ -String gridDistances(JdbcDataSource datasource, String input_polygons, String grid, String id_grid) { +String gridDistances(JdbcDataSource datasource, String input_polygons, String grid, String id_grid, boolean keep_geometry=true) { if (!input_polygons) { error("The input polygons cannot be null or empty") return @@ -294,24 +294,49 @@ String gridDistances(JdbcDataSource datasource, String input_polygons, String gr int epsg = datasource.getSrid(grid) def outputTableName = postfix("grid_distances") - datasource.execute(""" DROP TABLE IF EXISTS $outputTableName; - CREATE TABLE $outputTableName (THE_GEOM GEOMETRY,ID INT, DISTANCE FLOAT); + if(keep_geometry) { + datasource.execute(""" DROP TABLE IF EXISTS $outputTableName; + CREATE TABLE $outputTableName (THE_GEOM GEOMETRY,$id_grid INT, DISTANCE FLOAT); """.toString()) - datasource.createSpatialIndex(input_polygons) - datasource.createSpatialIndex(grid) + datasource.createSpatialIndex(input_polygons) + datasource.createSpatialIndex(grid) - datasource.withBatch(100) { stmt -> - datasource.eachRow("SELECT the_geom from $input_polygons".toString()) { row -> - Geometry geom = row.the_geom - if (geom) { - IndexedFacetDistance indexedFacetDistance = new IndexedFacetDistance(geom) - datasource.eachRow("""SELECT the_geom, ${id_grid} as id from $grid + datasource.withBatch(100) { stmt -> + datasource.eachRow("SELECT the_geom from $input_polygons".toString()) { row -> + Geometry geom = row.the_geom + if (geom) { + IndexedFacetDistance indexedFacetDistance = new IndexedFacetDistance(geom) + datasource.eachRow("""SELECT the_geom, ${id_grid} as id from $grid where ST_GEOMFROMTEXT('${geom}',$epsg) && the_geom and st_intersects(ST_GEOMFROMTEXT('${geom}',$epsg) , ST_POINTONSURFACE(the_geom))""".toString()) { cell -> - Geometry cell_geom = cell.the_geom - double distance = indexedFacetDistance.distance(cell_geom.getCentroid()) - stmt.addBatch "insert into $outputTableName values(ST_GEOMFROMTEXT('${cell_geom}',$epsg), ${cell.id},${distance})".toString() + Geometry cell_geom = cell.the_geom + double distance = indexedFacetDistance.distance(cell_geom.getCentroid()) + stmt.addBatch "insert into $outputTableName values(ST_GEOMFROMTEXT('${cell_geom}',$epsg), ${cell.id},${distance})".toString() + } + } + } + } + }else{ + datasource.execute(""" DROP TABLE IF EXISTS $outputTableName; + CREATE TABLE $outputTableName ($id_grid INT, DISTANCE FLOAT); + """.toString()) + + datasource.createSpatialIndex(input_polygons) + datasource.createSpatialIndex(grid) + + datasource.withBatch(100) { stmt -> + datasource.eachRow("SELECT the_geom from $input_polygons".toString()) { row -> + Geometry geom = row.the_geom + if (geom) { + IndexedFacetDistance indexedFacetDistance = new IndexedFacetDistance(geom) + datasource.eachRow("""SELECT the_geom, ${id_grid} as id from $grid + where ST_GEOMFROMTEXT('${geom}',$epsg) && the_geom and + st_intersects(ST_GEOMFROMTEXT('${geom}',$epsg) , ST_POINTONSURFACE(the_geom))""".toString()) { cell -> + Geometry cell_geom = cell.the_geom + double distance = indexedFacetDistance.distance(cell_geom.getCentroid()) + stmt.addBatch "insert into $outputTableName values(${cell.id},${distance})".toString() + } } } } diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index 33b3cb5241..7c8075d693 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -1803,7 +1803,7 @@ String rasterizeIndicators(JdbcDataSource datasource, /* * Make aggregation process with previous grid and current rsu lcz */ - if (list_indicators_upper.intersect(["LCZ_FRACTION", "LCZ_PRIMARY"]) && rsu_lcz) { + if (list_indicators_upper.intersect(["LCZ_FRACTION", "LCZ_PRIMARY", "SPRAWL_AREAS", "SPRAWL_DISTANCES", "SPRAWL_COOL_DISTANCE"]) && rsu_lcz) { def indicatorName = "LCZ_PRIMARY" String upperScaleAreaStatistics = Geoindicators.GenericIndicators.upperScaleAreaStatistics( datasource, grid, grid_column_identifier, @@ -2243,7 +2243,8 @@ String rasterizeIndicators(JdbcDataSource datasource, // Remove temporary tables datasource.dropTable(tablesToDrop) - if(lcz_lod){ + //We must compute the LOC indicators if the user also wants sprawl indicators + if(lcz_lod || list_indicators_upper.intersect(["SPRAWL_AREAS", "SPRAWL_DISTANCES", "SPRAWL_COOL_DISTANCE"])){ return Geoindicators.GridIndicators.multiscaleLCZGrid(datasource, grid_indicators_table,grid_column_identifier,lcz_lod) } return grid_indicators_table @@ -2452,4 +2453,79 @@ static boolean modelCheck(String modelName) { } } return true +} + +/*** + * This method is used to compute a set of sprawl_areas indicators as + * - the sprawl_areas layer + * - the distances inside and outside the sprawl_areas + * - the distance to cool areas (LCZ 101, 102, 103, 104,106,107) inside the sprawl_areas + * + * @param datasource + * @param grid_indicators + * @param list_indicators + * @param distance the erode and dilate the geometries + * @return the sprawl_areas layer plus new distance columns on the input grid_indicators + */ +Map sprawlIndicators(JdbcDataSource datasource, String grid_indicators, String id_grid, List list_indicators, float distance) { + if (!list_indicators) { + info "The list of indicator names cannot be null or empty" + return + } + + //Concert the list of indicators to upper case + allowed_indicators = ["SPRAWL_AREAS", "SPRAWL_DISTANCES", "SPRAWL_COOL_DISTANCE"] + def list_indicators_upper = list_indicators.collect { it.toUpperCase() } + + def tablesToDrop = [] + def tablesToJoin = [:] + tablesToJoin.put(grid_indicators,id_grid) + String sprawl_areas + if (list_indicators_upper.intersect(["SPRAWL_AREAS", "SPRAWL_DISTANCES", "SPRAWL_COOL_DISTANCE"]) && grid_indicators) { + sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(datasource, grid_indicators, distance) + } + if (sprawl_areas) { + //Compute the distances + if (list_indicators_upper.contains("SPRAWL_DISTANCES")) { + String inside_sprawl_areas = Geoindicators.GridIndicators.gridDistances(datasource, sprawl_areas, grid_indicators, id_grid, false) + if (inside_sprawl_areas) { + datasource.execute("""ALTER TABLE $inside_sprawl_areas RENAME COLUMN DISTANCE TO SPRAWL_INDIST""".toString()) + tablesToDrop << inside_sprawl_areas + String inverse_sprawl_areas = Geoindicators.SpatialUnits.inversePolygonsLayer(datasource, sprawl_areas) + if (inverse_sprawl_areas) { + tablesToDrop << inverse_sprawl_areas + tablesToJoin.put(inside_sprawl_areas, id_grid) + String outside_sprawl_areas = Geoindicators.GridIndicators.gridDistances(datasource, inverse_sprawl_areas, grid_indicators, id_grid, false) + if (outside_sprawl_areas) { + datasource.execute("""ALTER TABLE $outside_sprawl_areas RENAME COLUMN DISTANCE TO SPRAWL_OUTDIST""".toString()) + tablesToDrop << outside_sprawl_areas + tablesToJoin.put(outside_sprawl_areas, id_grid) + } + } + } + } + if (list_indicators_upper.contains("SPRAWL_COOL_DISTANCE")) { + String cool_areas = Geoindicators.SpatialUnits.extractCoolAreas(datasource, grid_indicators, distance) + if (cool_areas) { + tablesToDrop << cool_areas + String inverse_cool_areas = Geoindicators.SpatialUnits.inversePolygonsLayer(datasource, cool_areas) + if (inverse_cool_areas) { + tablesToDrop << inverse_cool_areas + String cool_distances = Geoindicators.GridIndicators.gridDistances(datasource, inverse_cool_areas, grid_indicators, id_grid, false) + if (cool_distances) { + tablesToDrop << cool_distances + datasource.execute("""ALTER TABLE $cool_distances RENAME COLUMN DISTANCE TO SPRAWL_COOL_INDIST""".toString()) + tablesToJoin.put(cool_distances, id_grid) + } + } + } + } + } + if (tablesToJoin.size() > 1) { + tablesToDrop<