Skip to content

Commit

Permalink
Merge pull request #941 from ebocher/osm_sea_land_mask_nantes
Browse files Browse the repository at this point in the history 
Fix sea/land mask with grid indicators
  • Loading branch information
@ebocher
ebocher authored Feb 23, 2024
2 parents cf93684 + 43a04e8 commit af880c7
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 43 deletions.
80 changes: 46 additions & 34 deletions ...cators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1027,7 +1027,7 @@ Map computeTypologyIndicators(JdbcDataSource datasource, String building_indicat
*/
Map createUnitsOfAnalysis(JdbcDataSource datasource, String zone, String building,
String road, String rail, String vegetation,
String water, String sea_land_mask, String urban_areas,
String water, String sea_land_mask, String urban_areas,
String rsu, double surface_vegetation,
double surface_hydro, double surface_urban_areas,
double snappingTolerance, List indicatorUse = ["LCZ", "UTRF", "TEB"], String prefixName = "") {
Expand All @@ -1039,7 +1039,7 @@ Map createUnitsOfAnalysis(JdbcDataSource datasource, String zone, String buildin
rsu = Geoindicators.SpatialUnits.createTSU(datasource, zone, road, rail,
vegetation, water,
sea_land_mask, urban_areas, surface_vegetation,
surface_hydro,surface_urban_areas, prefixName)
surface_hydro, surface_urban_areas, prefixName)
if (!rsu) {
info "Cannot compute the RSU."
return
Expand Down Expand Up @@ -1112,7 +1112,7 @@ Map createUnitsOfAnalysis(JdbcDataSource datasource, String zone, String buildin
*/
Map getParameters() {
return [
"surface_vegetation" : 10000d, "surface_hydro": 2500d, "surface_urban_areas":10000d,
"surface_vegetation" : 10000d, "surface_hydro": 2500d, "surface_urban_areas": 10000d,
"snappingTolerance" : 0.01d, "indicatorUse": ["LCZ", "UTRF", "TEB"],
"mapOfWeights" : ["sky_view_factor" : 1, "aspect_ratio": 1, "building_surface_fraction": 1,
"impervious_surface_fraction" : 1, "pervious_surface_fraction": 1,
Expand Down Expand Up @@ -1210,7 +1210,7 @@ Map getParameters(Map parameters) {
*/
Map computeAllGeoIndicators(JdbcDataSource datasource, String zone, String building, String road, String rail, String vegetation,
String water, String impervious, String buildingEstimateTableName,
String sea_land_mask,String urban_areas, String rsuTable,
String sea_land_mask, String urban_areas, String rsuTable,
Map parameters = [:], String prefixName) {
Map inputParameters = getParameters()
if (parameters) {
Expand Down Expand Up @@ -1241,7 +1241,7 @@ Map computeAllGeoIndicators(JdbcDataSource datasource, String zone, String build
water, impervious,
buildingEstimateTableName,
sea_land_mask, urban_areas, rsuTable,
surface_vegetation, surface_hydro,surface_urban_areas,
surface_vegetation, surface_hydro, surface_urban_areas,
snappingTolerance,
buildingHeightModelName, prefixName)
if (!estimHeight) {
Expand Down Expand Up @@ -1330,9 +1330,9 @@ Map computeAllGeoIndicators(JdbcDataSource datasource, String zone, String build
Map spatialUnits = createUnitsOfAnalysis(datasource, zone,
building, road,
rail, vegetation,
water, sea_land_mask, "","",
water, sea_land_mask, "", "",
surface_vegetation,
surface_hydro,surface_urban_areas, snappingTolerance, indicatorUse,
surface_hydro, surface_urban_areas, snappingTolerance, indicatorUse,
prefixName)
if (!spatialUnits) {
error "Cannot create the spatial units"
Expand Down Expand Up @@ -1368,7 +1368,7 @@ Map estimateBuildingHeight(JdbcDataSource datasource, String zone, String buildi
String road, String rail, String vegetation,
String water, String impervious,
String building_estimate, String sea_land_mask, String urban_areas, String rsu,
double surface_vegetation, double surface_hydro,double surface_urban_areas,
double surface_vegetation, double surface_hydro, double surface_urban_areas,
double snappingTolerance, String buildingHeightModelName, String prefixName = "") {
if (!building_estimate) {
error "To estimate the building height a table that contains the list of building to estimate must be provided"
Expand All @@ -1384,7 +1384,7 @@ Map estimateBuildingHeight(JdbcDataSource datasource, String zone, String buildi
Map spatialUnits = createUnitsOfAnalysis(datasource, zone,
building, road, rail, vegetation,
water, sea_land_mask, urban_areas, rsu,
surface_vegetation, surface_hydro,surface_urban_areas, snappingTolerance, ["UTRF"],
surface_vegetation, surface_hydro, surface_urban_areas, snappingTolerance, ["UTRF"],
prefixName)
if (!spatialUnits) {
error "Cannot create the spatial units"
Expand Down Expand Up @@ -1435,7 +1435,7 @@ Map estimateBuildingHeight(JdbcDataSource datasource, String zone, String buildi
prefixName)

def buildingTableName = "BUILDING_TABLE_WITH_RSU_AND_BLOCK_ID"
int nbBuildingEstimated =0
int nbBuildingEstimated = 0
def buildEstimatedHeight
if (datasource.getTable(gatheredScales).isEmpty()) {
info "No building height to estimate"
Expand All @@ -1449,7 +1449,7 @@ Map estimateBuildingHeight(JdbcDataSource datasource, String zone, String buildi
info "Start estimating the building height"
//Apply RF model
buildEstimatedHeight = Geoindicators.TypologyClassification.applyRandomForestModel(datasource,
gatheredScales, buildingHeightModelName,"id_build", prefixName)
gatheredScales, buildingHeightModelName, "id_build", prefixName)
if (!buildEstimatedHeight) {
error "Cannot apply the building height model $buildingHeightModelName"
return
Expand Down Expand Up @@ -1645,6 +1645,9 @@ String rasterizeIndicators(JdbcDataSource datasource,
def grid_column_identifier = "id_grid"
def indicatorTablesToJoin = [:]
indicatorTablesToJoin.put(grid, grid_column_identifier)

//An array to execute some commands on the final table
def sqlUpdateCommand =[]
/*
* Make aggregation process with previous grid and current rsu lcz
*/
Expand Down Expand Up @@ -1845,7 +1848,7 @@ String rasterizeIndicators(JdbcDataSource datasource,
"building_type_fraction")
if (upperScaleAreaStatistics) {
indicatorTablesToJoin.put(upperScaleAreaStatistics, grid_column_identifier)
tablesToDrop<<upperScaleAreaStatistics
tablesToDrop << upperScaleAreaStatistics
} else {
info "Cannot aggregate the building type at grid scale"
}
Expand All @@ -1868,7 +1871,7 @@ String rasterizeIndicators(JdbcDataSource datasource,
if (freeFacadeDensityExact) {
if (list_indicators_upper.intersect(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO"])) {
indicatorTablesToJoin.put(freeFacadeDensityExact, grid_column_identifier)
tablesToDrop<<freeFacadeDensityExact
tablesToDrop << freeFacadeDensityExact
}
if (list_indicators_upper.contains("FREE_EXTERNAL_FACADE_DENSITY")) {
def buildingSurfDensity = Geoindicators.RsuIndicators.buildingSurfaceDensity(
Expand All @@ -1881,8 +1884,8 @@ String rasterizeIndicators(JdbcDataSource datasource,
indicatorTablesToJoin.put(buildingSurfDensity, grid_column_identifier)
}
}
tablesToDrop<<freeFacadeDensityExact
tablesToDrop<<buildingSurfDensity
tablesToDrop << freeFacadeDensityExact
tablesToDrop << buildingSurfDensity
}
} else {
info "Cannot calculate the exact free external facade density"
Expand Down Expand Up @@ -1932,20 +1935,21 @@ String rasterizeIndicators(JdbcDataSource datasource,
}

if (list_indicators_upper.contains("SEA_LAND_FRACTION") && sea_land_mask) {
// If only one type of surface (land or sea) is in the zone, no need for computational fraction calculation
def sea_land_type_rows = datasource.rows("""SELECT $seaLandTypeField, COUNT(*) AS NB_TYPES
if (datasource.getRowCount(sea_land_mask) > 0) {
// If only one type of surface (land or sea) is in the zone, no need for computational fraction calculation
def sea_land_type_rows = datasource.rows("""SELECT $seaLandTypeField, COUNT(*) AS NB_TYPES
FROM $sea_land_mask
GROUP BY $seaLandTypeField""".toString())
if (sea_land_type_rows[seaLandTypeField].count("sea") == 0) {
datasource """
if (sea_land_type_rows[seaLandTypeField].count("sea") == 0) {
datasource """
DROP TABLE IF EXISTS $seaLandFractionTab;
CREATE TABLE $seaLandFractionTab
AS SELECT $grid_column_identifier, 1 AS LAND_FRACTION
FROM $grid"""
indicatorTablesToJoin.put(seaLandFractionTab, grid_column_identifier)
} else {
// Split the potentially big complex seaLand geometries into smaller triangles in order to makes calculation faster
datasource """
indicatorTablesToJoin.put(seaLandFractionTab, grid_column_identifier)
} else {
// Split the potentially big complex seaLand geometries into smaller triangles in order to makes calculation faster
datasource """
DROP TABLE IF EXISTS $tesselatedSeaLandTab;
CREATE TABLE $tesselatedSeaLandTab(id_tesselate serial, the_geom geometry, $seaLandTypeField VARCHAR)
AS SELECT explod_id, the_geom, $seaLandTypeField
Expand All @@ -1954,21 +1958,26 @@ String rasterizeIndicators(JdbcDataSource datasource,
WHERE ST_DIMENSION(the_geom) = 2 AND ST_ISEMPTY(the_geom) IS NOT TRUE
AND ST_AREA(the_geom)>0)')"""

def upperScaleAreaStatistics = Geoindicators.GenericIndicators.upperScaleAreaStatistics(datasource,
grid, grid_column_identifier,
tesselatedSeaLandTab, seaLandTypeField, seaLandTypeField,
prefixName)
tablesToDrop << tesselatedSeaLandTab
if (upperScaleAreaStatistics) {
// Modify columns name to postfix with "_FRACTION"
datasource """
def upperScaleAreaStatistics = Geoindicators.GenericIndicators.upperScaleAreaStatistics(datasource,
grid, grid_column_identifier,
tesselatedSeaLandTab, seaLandTypeField, seaLandTypeField,
prefixName)
tablesToDrop << tesselatedSeaLandTab
if (upperScaleAreaStatistics) {
// Modify columns name to postfix with "_FRACTION"
datasource """
ALTER TABLE ${upperScaleAreaStatistics} RENAME COLUMN TYPE_LAND TO LAND_FRACTION;
ALTER TABLE ${upperScaleAreaStatistics} RENAME COLUMN TYPE_SEA TO SEA_FRACTION;
ALTER TABLE ${upperScaleAreaStatistics} DROP COLUMN THE_GEOM;"""
indicatorTablesToJoin.put(upperScaleAreaStatistics, grid_column_identifier)
} else {
info "Cannot compute the frontal area index."
indicatorTablesToJoin.put(upperScaleAreaStatistics, grid_column_identifier)
} else {
info "Cannot compute the sea land fractions."
}
}
}else{
//Update the final table
sqlUpdateCommand<<"""ALTER TABLE $grid_indicators_table ADD COLUMN LAND_FRACTION DOUBLE PRECISION DEFAULT 1;
ALTER TABLE $grid_indicators_table ADD COLUMN SEA_FRACTION DOUBLE PRECISION DEFAULT 0;"""
}
}

Expand Down Expand Up @@ -2072,6 +2081,9 @@ String rasterizeIndicators(JdbcDataSource datasource,
UPDATE $grid_indicators_table SET ASPECT_RATIO = case when building_fraction = 1 then null else free_external_facade_density / (1 - building_fraction) end
""".toString())
}
//Execute commands
datasource.execute(sqlUpdateCommand.join(" "))

tablesToDrop << createScalesRelationsGridBl
tablesToDrop << buildingCutted

Expand Down
23 changes: 14 additions & 9 deletions osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -652,7 +652,7 @@ class WorflowOSMTest extends WorkflowAbstractTest {
def location = "Redon"
//def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData("Redon")
// location = nominatim.bbox
//location=[47.4, -4.8, 47.6, -4.6]
location=[47.190062,-1.614389,47.196959,-1.602330]
def osm_parmeters = [
"description" : "Example of configuration file to run the OSM workflow and store the result in a folder",
"geoclimatedb": [
Expand All @@ -674,21 +674,21 @@ class WorflowOSMTest extends WorkflowAbstractTest {
["distance" : 0,
"rsu_indicators" : [

"indicatorUse": ["LCZ"] //, "UTRF", "TEB"]
"indicatorUse": ["LCZ", "UTRF", "TEB"]

]/*,"grid_indicators": [
"x_size": 200,
"y_size": 200,
],"grid_indicators": [
"x_size": 100,
"y_size": 100,
//"rowCol": true,
"indicators": ["BUILDING_FRACTION","BUILDING_HEIGHT", "BUILDING_POP",
//"BUILDING_TYPE_FRACTION",
"BUILDING_TYPE_FRACTION",
"WATER_FRACTION","VEGETATION_FRACTION",
"ROAD_FRACTION", "IMPERVIOUS_FRACTION",
"LCZ_PRIMARY",
//"BUILDING_HEIGHT_WEIGHTED", //"BUILDING_SURFACE_DENSITY", "SEA_LAND_FRACTION",
"ASPECT_RATIO",//"SVF",
"BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", "SEA_LAND_FRACTION",
"ASPECT_RATIO","SVF",
"HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"]
], "worldpop_indicators": true,
]/*, "worldpop_indicators": true,
"road_traffic" : true,
"noise_indicators" : [
"ground_acoustic": true
Expand Down Expand Up @@ -870,7 +870,11 @@ class WorflowOSMTest extends WorkflowAbstractTest {
assertTrue h2gis.firstRow("select count(*) as count from $building where HEIGHT_WALL>0 and HEIGHT_ROOF>0").count > 0
}

@Disabled
@Test
//TODO : A fix on OSM must be done because Golfe de Gascogne is tagged as a river.
//This geometry must be redesigned because according the International Hydrographic Organization
// Golfe de Gascogne must be considered as bay where some main rivers empty into it
void testOneSeaLCZ() {
String directory = folder.absolutePath + File.separator + "test_sea_lcz"
File dirFile = new File(directory)
Expand All @@ -893,6 +897,7 @@ class WorflowOSMTest extends WorkflowAbstractTest {
def tableNames = process.values()
def lcz = tableNames.rsu_lcz[0]
H2GIS h2gis = H2GIS.open("${directory + File.separator}sea_lcz_db;AUTO_SERVER=TRUE")
h2gis.save(lcz, "/tmp/sea.geojson", true)
def lcz_group= h2gis.firstRow("select lcz_primary, count(*) as count from $lcz group by lcz_primary".toString())
assertTrue(lcz_group.size()==2)
assertTrue(lcz_group.lcz_primary==107)
Expand Down

0 comments on commit af880c7

Please sign in to comment.