Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add all sprawl indicators : areas, distances, cool distance #950

Merged
merged 4 commits into from
Apr 5, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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}"
Expand Down Expand Up @@ -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")
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
]
]
]
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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()
}
}
}
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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<<grid_indicators
grid_indicators = Geoindicators.DataUtils.joinTables(datasource, tablesToJoin, postfix("grid_indicators"))
}
datasource.dropTable(tablesToDrop)

return ["sprawl_areas":sprawl_areas, "grid_indicators": grid_indicators]
}
Original file line number Diff line number Diff line change
Expand Up @@ -658,11 +658,11 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d
if (rasterizedIndicators) {
h2gis_datasource.dropTable(grid)
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)
}
def 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")
}
Expand Down Expand Up @@ -911,7 +911,8 @@ def extractProcessingParameters(def processing_parameters) {
"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"]
"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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -649,7 +649,7 @@ class WorflowOSMTest extends WorkflowAbstractTest {
File dirFile = new File(directory)
dirFile.delete()
dirFile.mkdir()
def location = "Redon"
def location = "Kervignac"
def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(location)
def grid_size = 100
location = nominatim.bbox
Expand Down Expand Up @@ -690,9 +690,10 @@ class WorflowOSMTest extends WorkflowAbstractTest {
"SEA_LAND_FRACTION",
"ASPECT_RATIO",
//"SVF",
"HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"],
"lcz_lod":2,
"sprawl_areas":true
"HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS",
"SPRAWL_AREAS",
"SPRAWL_DISTANCES", "SPRAWL_COOL_DISTANCE"],
//"lcz_lod":2
]/*, "worldpop_indicators": true,

"road_traffic" : true,
Expand Down
Loading