Skip to content

Commit

Permalink
Fix free external density indicators when any building share a facade
Browse files Browse the repository at this point in the history
Add rasterizeIndicators test
  • Loading branch information
ebocher committed Oct 4, 2023
1 parent c8e51b4 commit 0efec82
Show file tree
Hide file tree
Showing 4 changed files with 95 additions and 6 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -170,13 +170,24 @@ String freeExternalFacadeDensityExact(JdbcDataSource datasource, String building
// 4. Calculates the free facade density by RSU (subtract 3 and 2 and divide by RSU area)
datasource.createIndex(buildLineRsu,idRsu)
datasource.createIndex(sharedLineRsu,idRsu)
datasource """
if(datasource.getRowCount(sharedLineRsu)>0) {
datasource """
DROP TABLE IF EXISTS $onlyBuildRsu;
CREATE TABLE $onlyBuildRsu
AS SELECT a.$idRsu,
(a.$FACADE_AREA-b.$FACADE_AREA)/a.$RSU_AREA AS FREE_EXTERNAL_FACADE_DENSITY
FROM $buildLineRsu AS a LEFT JOIN $sharedLineRsu AS b
ON a.$idRsu = b.$idRsu""".toString()
ON a.$idRsu = b.$idRsu
""".toString()
}else{
datasource """
DROP TABLE IF EXISTS $onlyBuildRsu;
CREATE TABLE $onlyBuildRsu
AS SELECT $idRsu,
$FACADE_AREA/$RSU_AREA AS FREE_EXTERNAL_FACADE_DENSITY
FROM $buildLineRsu
""".toString()
}

// 5. Join RSU having no buildings and set their value to 0
datasource.createIndex(onlyBuildRsu,idRsu)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1613,7 +1613,7 @@ Map computeGeoclimateIndicators(JdbcDataSource datasource, String zone, String b
* @return
*/
String rasterizeIndicators(JdbcDataSource datasource,
String grid, List list_indicators = [],
String grid, List list_indicators,
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,
Expand Down Expand Up @@ -1822,7 +1822,7 @@ String rasterizeIndicators(JdbcDataSource datasource,

}

if (list_indicators_upper.contains("BUILDING_TYPE") && building) {
if (list_indicators_upper.contains("BUILDING_TYPE_FRACTION") && building) {
if (!buildingCutted) {
buildingCutted = cutBuilding(datasource, grid, building)
if (!buildingCutted) {
Expand All @@ -1838,6 +1838,7 @@ String rasterizeIndicators(JdbcDataSource datasource,
"building_type_fraction")
if (upperScaleAreaStatistics) {
indicatorTablesToJoin.put(upperScaleAreaStatistics, grid_column_identifier)
tablesToDrop<<upperScaleAreaStatistics
} else {
info "Cannot aggregate the building type at grid scale"
}
Expand All @@ -1860,6 +1861,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
}
if (list_indicators_upper.contains("FREE_EXTERNAL_FACADE_DENSITY")) {
def buildingSurfDensity = Geoindicators.RsuIndicators.buildingSurfaceDensity(
Expand All @@ -1872,6 +1874,8 @@ String rasterizeIndicators(JdbcDataSource datasource,
indicatorTablesToJoin.put(buildingSurfDensity, grid_column_identifier)
}
}
tablesToDrop<<freeFacadeDensityExact
tablesToDrop<<buildingSurfDensity
}
} else {
info "Cannot calculate the exact free external facade density"
Expand Down Expand Up @@ -1976,7 +1980,7 @@ String rasterizeIndicators(JdbcDataSource datasource,
def svf_fraction = Geoindicators.RsuIndicators.groundSkyViewFactor(datasource, grid, grid_column_identifier, createScalesRelationsGridBl, 0.008, 100, 60, "grid")

if (svf_fraction) {
datasource """ ALTER TABLE ${svf_fraction} RENAME COLUMN GROUND_SKY_VIEW_FACTOR TO SVF_FRACTION""".toString()
datasource """ ALTER TABLE ${svf_fraction} RENAME COLUMN GROUND_SKY_VIEW_FACTOR TO SVF""".toString()
indicatorTablesToJoin.put(svf_fraction, grid_column_identifier)
} else {
info "Cannot compute the sky view factor."
Expand Down Expand Up @@ -2058,7 +2062,7 @@ String rasterizeIndicators(JdbcDataSource datasource,
//Because all other indicators have been yet computed we just alter the table with the aspect_ratio column
datasource.execute("""
ALTER TABLE $grid_indicators_table ADD COLUMN ASPECT_RATIO DOUBLE PRECISION;
UPDATE $grid_indicators_table SET ASPECT_RATIO = free_external_facade_density / (1 - building_fraction)
UPDATE $grid_indicators_table SET ASPECT_RATIO = case when building_fraction = 1 then free_external_facade_density else free_external_facade_density / (1 - building_fraction) end
""".toString())
}
tablesToDrop << createScalesRelationsGridBl
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 0efec82

Please sign in to comment.