Merge pull request #1003 from j3r3m1/street_width

Street width

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()

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)
+    }
 }

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,
@@ -1805,7 +1837,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)
         }
@@ -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(" "))
 

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

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)
     }

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)