From 4c249ca35cdc858ae98146f1435038d108cba39f Mon Sep 17 00:00:00 2001 From: BERNARD Jeremy Date: Wed, 6 Nov 2024 15:12:52 +0100 Subject: [PATCH 1/2] Modify indicator name from building_height_dist to ...distribution --- .../org/orbisgis/geoclimate/bdtopo/WorkflowDebugTest.groovy | 2 +- .../geoclimate/geoindicators/WorkflowGeoIndicators.groovy | 2 +- .../groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) 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 7e727d8cff..6ada845a8c 100644 --- a/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowDebugTest.groovy +++ b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowDebugTest.groovy @@ -138,7 +138,7 @@ class WorkflowDebugTest { void testIntegrationFolderInput() { def input_data = "/home/bernardj/Data/BDT/V3/BDTOPO_3-3_TOUSTHEMES_SHP_LAMB93_R11_2024-03-15" def locations = [ - [6864452.136340265, 651978.2850116278, 6865452.136340265, 652978.2850116278]] + [6857743.1274865065, 653030.9684552355, 6858743.1274865065, 654030.9684552355]] String directory = "/tmp/bdtopo3" File dirFile = new File(directory) dirFile.delete() 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 5295583b10..be48d24ac0 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -1805,7 +1805,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } - if (list_indicators_upper.contains("BUILDING_HEIGHT_DIST") && building) { + if (list_indicators_upper.contains("BUILDING_HEIGHT_DISTRIBUTION") && building) { if (!buildingCutted) { buildingCutted = cutBuilding(datasource, grid, building) } diff --git a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy index fbea159870..fdc30a15bc 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -733,7 +733,7 @@ class WorflowOSMTest extends WorkflowAbstractTest { "BUILDING_TYPE_FRACTION", "WATER_FRACTION", "VEGETATION_FRACTION", "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "FREE_EXTullERNAL_FACADE_DENSITY", "BUILDING_SURFACE_DENSITY", - "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION"] + "BUILDING_HEIGHT_DISTRIBUTION", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION"] def databasePath = '/tmp/geoclimate_chain_dbtest;AUTO_SERVER=TRUE' def h2gis_properties = ["databaseName": databasePath, "user": "sa", "password": ""] def datasource = H2GIS.open(h2gis_properties) From eba16b5d4ec57079da121199b4ca4b43a86ac730 Mon Sep 17 00:00:00 2001 From: BERNARD Jeremy Date: Wed, 6 Nov 2024 16:54:03 +0100 Subject: [PATCH 2/2] Add street_width in grid indic and modify aspect ratio --- .../geoindicators/RsuIndicators.groovy | 41 +++++++++++++++++ .../WorkflowGeoIndicators.groovy | 46 ++++++++++++++----- .../geoindicators/RsuIndicatorsTests.groovy | 21 +++++++++ .../WorkflowGeoIndicatorsTest.groovy | 4 +- 4 files changed, 100 insertions(+), 12 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy index f1acd23dfb..9dd155bd25 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -2471,4 +2471,45 @@ String groundLayer(JdbcDataSource datasource, String zone, String id_zone, error """Cannot compute the unified ground layer""" } return outputTableName +} + +/** + * Process used to compute the average street width needed by models such as TARGET. A simple approach based + * on the street canyons assumption is used for the calculation. The area weighted mean building height is divided + * by the aspect ratio (defined as the sum of facade area within a given RSU area divided by the area of free + * surfaces of the given RSU (not covered by buildings). The "avg_height_roof_area_weighted", + * "rsu_free_external_facade_density" and "rsu_building_density" are used for the calculation. + * + * @param datasource A connexion to a database (H2GIS, PostGIS, ...) where are stored the input table and in which + * the resulting database will be stored + * @param rsuTable The name of the input ITable where are stored the RSU + * @param avgHeightRoofAreaWeightedColumn The name of the column where are stored the area weighted mean building height + * values (within the rsu Table) + * @param aspectRatioColumn The name of the column where are stored the aspect ratio values (within the rsu Table) + * @param prefixName String use as prefix to name the output table + * + * @return A database table name. + * @author Jérémy Bernard + */ +String streetWidth(JdbcDataSource datasource, String rsuTable, String avgHeightRoofAreaWeightedColumn, + String aspectRatioColumn, prefixName) throws Exception { + try { + def COLUMN_ID_RSU = "id_rsu" + def BASE_NAME = "street_width" + debug "Executing RSU street width" + + // The name of the outputTableName is constructed + def outputTableName = prefix prefixName, "rsu_" + BASE_NAME + datasource.execute(""" + DROP TABLE IF EXISTS $outputTableName; + CREATE TABLE $outputTableName AS + SELECT CASE WHEN $aspectRatioColumn = 0 + THEN null + ELSE $avgHeightRoofAreaWeightedColumn / $aspectRatioColumn END + AS $BASE_NAME, $COLUMN_ID_RSU FROM $rsuTable""") + + return outputTableName + } catch (SQLException e) { + throw new SQLException("Cannot compute the street width at RSU scale", e) + } } \ No newline at end of file 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 be48d24ac0..85ef13d3d1 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -1592,6 +1592,8 @@ String rasterizeIndicators(JdbcDataSource datasource, def grid_column_identifier = "id_grid" def indicatorTablesToJoin = [:] indicatorTablesToJoin.put(grid, grid_column_identifier) + // Needed for intermediate calculations... + def tablesToJoinForWidth = [:] //An array to execute some commands on the final table def sqlUpdateCommand = [] @@ -1685,7 +1687,7 @@ String rasterizeIndicators(JdbcDataSource datasource, list_indicators_upper.each { if (it == "BUILDING_FRACTION" || it == "BUILDING_SURFACE_DENSITY" || - it == "ASPECT_RATIO" || it == "FREE_EXTERNAL_FACADE_DENSITY") { + it == "ASPECT_RATIO" || it == "FREE_EXTERNAL_FACADE_DENSITY" || it == "STREET_WIDTH") { columnFractionsList.put(priorities.indexOf("building"), "building") } else if (it == "WATER_FRACTION") { columnFractionsList.put(priorities.indexOf("water"), "water") @@ -1723,6 +1725,7 @@ String rasterizeIndicators(JdbcDataSource datasource, datasource, grid, grid_column_identifier, superpositionsTableGrid, [:], priorities_tmp, prefixName) indicatorTablesToJoin.put(surfaceFractionsProcess, grid_column_identifier) + tablesToJoinForWidth.put(surfaceFractionsProcess, grid_column_identifier) tablesToDrop << superpositionsTableGrid } else { @@ -1753,7 +1756,9 @@ String rasterizeIndicators(JdbcDataSource datasource, weightedBuildingIndicators, prefixName) indicatorTablesToJoin.put(computeWeightedAggregStat, grid_column_identifier) + tablesToJoinForWidth.put(computeWeightedAggregStat, grid_column_identifier) + tablesToDrop << computeWeightedAggregStat } if (list_indicators_upper.contains("BUILDING_TYPE_FRACTION") && building) { @@ -1771,7 +1776,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } - if ((list_indicators_upper.intersect(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO", "BUILDING_SURFACE_DENSITY"]) && building)) { + if (list_indicators_upper.intersect(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO", "BUILDING_SURFACE_DENSITY", "STREET_WIDTH"]) && building) { if (!createScalesRelationsGridBl) { // Create the relations between grid cells and buildings createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin(datasource, @@ -1780,11 +1785,38 @@ String rasterizeIndicators(JdbcDataSource datasource, } def freeFacadeDensityExact = Geoindicators.RsuIndicators.freeExternalFacadeDensityExact(datasource, createScalesRelationsGridBl, grid, grid_column_identifier, prefixName) + tablesToJoinForWidth.put(freeFacadeDensityExact, grid_column_identifier) if (freeFacadeDensityExact) { - if (list_indicators_upper.intersect(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO"])) { + if (list_indicators_upper.intersect(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO", "STREET_WIDTH"])) { indicatorTablesToJoin.put(freeFacadeDensityExact, grid_column_identifier) tablesToDrop << freeFacadeDensityExact } + + def gridForWidth = Geoindicators.DataUtils.joinTables(datasource, tablesToJoinForWidth, "grid_for_width") + tablesToDrop << gridForWidth + + // Compute the aspect_ratio (and street width if needed) + if (list_indicators_upper.intersect(["ASPECT_RATIO", "STREET_WIDTH"]) && building) { + datasource "ALTER TABLE $gridForWidth RENAME COLUMN $grid_column_identifier TO ID_RSU" + def aspectRatio = Geoindicators.RsuIndicators.aspectRatio(datasource, gridForWidth, + "free_external_facade_density", "building_fraction", prefixName) + datasource "ALTER TABLE $aspectRatio RENAME COLUMN ID_RSU TO $grid_column_identifier" + indicatorTablesToJoin.put(aspectRatio, grid_column_identifier) + tablesToDrop << aspectRatio + + // Calculates the street width if needed + if (list_indicators_upper.contains("STREET_WIDTH") && building) { + tablesToJoinForWidth.put(aspectRatio, grid_column_identifier) + gridForWidth = Geoindicators.DataUtils.joinTables(datasource, tablesToJoinForWidth, "grid_for_width") + datasource "ALTER TABLE $gridForWidth RENAME COLUMN $grid_column_identifier TO ID_RSU" + def streetWidth = Geoindicators.RsuIndicators.streetWidth(datasource, gridForWidth, + "avg_height_roof_area_weighted", "aspect_ratio", prefixName) + datasource "ALTER TABLE $streetWidth RENAME COLUMN ID_RSU TO $grid_column_identifier" + indicatorTablesToJoin.put(streetWidth, grid_column_identifier) + tablesToDrop << streetWidth + } + } + if (list_indicators_upper.contains("FREE_EXTERNAL_FACADE_DENSITY")) { def buildingSurfDensity = Geoindicators.RsuIndicators.buildingSurfaceDensity( datasource, freeFacadeDensityExact, @@ -1935,14 +1967,6 @@ String rasterizeIndicators(JdbcDataSource datasource, //Join all indicators at grid scale def joinGrids = Geoindicators.DataUtils.joinTables(datasource, indicatorTablesToJoin, grid_indicators_table) - //Compute the aspect_ratio - if (list_indicators_upper.contains("ASPECT_RATIO") && building) { - //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 = case when building_fraction = 1 then null else free_external_facade_density / (1 - building_fraction) end - """.toString()) - } //Execute commands datasource.execute(sqlUpdateCommand.join(" ")) diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy index 25b01cb47b..7ca07eb672 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy @@ -174,6 +174,27 @@ class RsuIndicatorsTests { assertEquals(null, result["aspect_ratio"]) } + @Test + void streetWidthTest() { + def p = Geoindicators.RsuIndicators.streetWidth(h2GIS, "rsu_test_all_indics_for_lcz", + "GEOM_AVG_HEIGHT_ROOF", "aspect_ratio", "test") + assertNotNull(p) + def concat = 0 + h2GIS.eachRow("SELECT * FROM test_rsu_street_width WHERE id_rsu = 1") { + row -> concat += row.street_width + } + assertEquals(7.5, concat, 0.001) + } + + @Test + void streetWidthTest2() { + def p = Geoindicators.RsuIndicators.streetWidth(h2GIS, "rsu_test_all_indics_for_lcz", + "GEOM_AVG_HEIGHT_ROOF", "aspect_ratio", "test") + assert (p) + def result = h2GIS.firstRow("SELECT street_width FROM test_rsu_street_width WHERE id_rsu = 18") + assertEquals(null, result["street_width"]) + } + @Test void projectedFacadeAreaDistributionTest() { // Only the first 5 first created buildings are selected for the tests diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy index f7c78ea120..0d83364dc8 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy @@ -531,7 +531,8 @@ class WorkflowGeoIndicatorsTest { "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"] + "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS", + "STREET_WIDTH"] grid_indicators = Geoindicators.WorkflowGeoIndicators.rasterizeIndicators(datasource, grid, list_indicators, "building", null, null, null, null, null, null, null, null) @@ -556,6 +557,7 @@ class WorkflowGeoIndicatorsTest { assertTrue(1.5 - rows.EFFECTIVE_TERRAIN_ROUGHNESS_LENGTH < 0.001) assertEquals(8, rows.EFFECTIVE_TERRAIN_ROUGHNESS_CLASS) assertNull(rows.ASPECT_RATIO) + assertNull(rows.STREET_WIDTH) assertTrue(0.5 - rows.SVF < 0.1) assertEquals(1d, rows.TYPE_OFFICE) }