Add new exact indicators at grid scale

geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy

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

geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy

@@ -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,
@@ -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) {
@@ -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"
         }
@@ -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(
@@ -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"
@@ -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."
@@ -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

geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy

@@ -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

geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy

@@ -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