New indicators to investigate LCZ classification

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

@@ -39,7 +39,6 @@ abstract class Geoindicators extends AbstractScript {
     public static PopulationIndicators = new PopulationIndicators()
     public static GridIndicators = new GridIndicators()
     public static NoiseIndicators = new NoiseIndicators()
-
     public static WorkflowUtilities = new WorkflowUtilities()
 
     //Cache the XStream models

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

@@ -971,7 +971,7 @@ String upperScaleAreaStatistics(JdbcDataSource datasource, String upperTableName
                                         st_force2d(st_makevalid(b.$upperGeometryColumn)))) AS area
                               FROM $lowerTableName a, $upperTableName b
                               WHERE a.$lowerGeometryColumn && b.$upperGeometryColumn AND 
-                              ST_INTERSECTS(st_force2d(a.$lowerGeometryColumn), st_force2d(b.$upperGeometryColumn));
+                              ST_INTERSECTS(a.$lowerGeometryColumn, b.$upperGeometryColumn);
                               """
     datasource.execute(spatialJoin.toString())
     datasource """CREATE INDEX ON $spatialJoinTable ($lowerColumnName);

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

@@ -20,6 +20,8 @@
 package org.orbisgis.geoclimate.geoindicators
 
 import groovy.transform.BaseScript
+import org.locationtech.jts.geom.Geometry
+import org.locationtech.jts.operation.distance.IndexedFacetDistance
 import org.orbisgis.data.jdbc.JdbcDataSource
 import org.orbisgis.geoclimate.Geoindicators
 
@@ -47,8 +49,8 @@ String gridPopulation(JdbcDataSource datasource, String gridTable, String popula
     def outputTableName = postfix BASE_NAME
 
     //Indexing table
-    datasource.createSpatialIndex(gridTable,"the_geom")
-    datasource.createSpatialIndex(populationTable,"the_geom")
+    datasource.createSpatialIndex(gridTable, "the_geom")
+    datasource.createSpatialIndex(populationTable, "the_geom")
     def popColumns = []
     def sum_popColumns = []
     if (populationColumns) {
@@ -93,4 +95,226 @@ String gridPopulation(JdbcDataSource datasource, String gridTable, String popula
             drop table if exists $gridTable_pop,$gridTable_pop_sum, $gridTable_area_sum ;""".toString())
 
     return outputTableName
-}
+}
+
+
+/**
+ * Create a multi-scale grid and aggregate the LCZ_PRIMARY indicators for each level of the grid.
+ * For each level, the adjacent cells are preserved as well as the number of urban and natural cells.
+ * To distinguish between LCZ cells with the same count per level, a weight is used corresponding
+ * to their potential impact on heat
+ *
+ * @param datasource connection to the database
+ * @param grid_indicators a grid that contains for each cell the LCZ_PRIMARY
+ * @param  id_grid grid cell column identifier
+ * @param nb_levels number of aggregate levels. Default is 1
+ *
+ * @return a the initial grid with all aggregated values by levels and the indexes (row, col) for each levels
+ * @author Erwan Bocher (CNRS)
+ */
+String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, String id_grid, int nb_levels = 1) {
+    if (!grid_indicators) {
+        error("No grid_indicators table to aggregate the LCZ values")
+        return
+    }
+    if (nb_levels <= 0 || nb_levels >= 10) {
+        error("The number of levels to aggregate the LCZ values must be between 1 and 10")
+        return
+    }
+
+    def gridColumns = datasource.getColumnNames(grid_indicators)
+
+    if(gridColumns.intersect(["LCZ_PRIMARY", "ID_ROW", "ID_COLUMN", id_grid]).size()==0){
+        error("The grid indicators table must contain the columns LCZ_PRIMARY, ID_ROW, $id_grid")
+        return
+    }
+
+    datasource.execute("""create index on $grid_indicators(id_row,id_col)""".toString())
+    ///First build the index levels for each cell of the input grid
+    def grid_scaling_indices = postfix("grid_scaling_indices")
+    def grid_levels_query = []
+    int grid_offset = 3
+    int offsetCol = 1
+    for (int i in 1..nb_levels) {
+        int level = Math.pow(grid_offset, i)
+        grid_levels_query << " (CAST (ABS(ID_ROW-1)/${level}  AS INT)+1) AS ID_ROW_LOD_$i," +
+                "(CAST (ABS(ID_COL-1)/${level} AS INT)+$offsetCol-1) AS ID_COL_LOD_$i"
+        offsetCol++
+    }
+
+    //Compute the indices for each levels and find the 8 adjacent cells
+    datasource.execute("""DROP TABLE IF EXISTS $grid_scaling_indices;
+    CREATE TABLE $grid_scaling_indices as SELECT *, ${grid_levels_query.join(",")},
+    (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL=a.ID_COL) AS LCZ_PRIMARY_N,
+    (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL=a.ID_COL+1) AS LCZ_PRIMARY_NE,
+    (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW AND ID_COL=a.ID_COL+1) AS LCZ_PRIMARY_E,
+    (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL+1) AS LCZ_PRIMARY_SE,
+    (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL) AS LCZ_PRIMARY_S,
+    (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_SW,
+    (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_W,
+    (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_NW
+    FROM 
+    $grid_indicators as a;  """.toString())
+    //Add LCZ_WARM count at this first level
+    def lcz_warm_first_level = postfix("lcz_warm_first_level")
+    datasource.execute("""DROP TABLE IF EXISTS $lcz_warm_first_level;
+    CREATE TABLE $lcz_warm_first_level as SELECT *,
+    (CASE WHEN LCZ_PRIMARY_N in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END +
+    CASE WHEN  LCZ_PRIMARY_NE in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END +
+    CASE WHEN  LCZ_PRIMARY_E in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END +
+    CASE WHEN  LCZ_PRIMARY_SE in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END +
+    CASE WHEN  LCZ_PRIMARY_S in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END +
+    CASE WHEN  LCZ_PRIMARY_SW in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END +
+    CASE WHEN  LCZ_PRIMARY_W in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END +
+    CASE WHEN  LCZ_PRIMARY_NW in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END+
+    CASE WHEN  LCZ_PRIMARY in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 ELSE 0 END)  AS LCZ_WARM 
+    FROM $grid_scaling_indices """.toString())
+
+    def tablesToDrop = []
+    def tableLevelToJoin = lcz_warm_first_level
+    tablesToDrop << grid_scaling_indices
+    tablesToDrop << lcz_warm_first_level
+
+    //Process all level of details
+    for (int i in 1..nb_levels) {
+        //Index the level row and col
+        datasource.execute("""
+        CREATE INDEX IF NOT EXISTS ${grid_scaling_indices}_idx ON $grid_scaling_indices(id_row_lod_${i},id_col_lod_${i})""".toString())
+        //First compute the number of cells by level of detail
+        //Use the original grid to aggregate the data
+        //A weight is used to select the LCZ value when the mode returns more than one possibility
+        def lcz_count_lod = postfix("lcz_count_lod")
+        tablesToDrop << lcz_count_lod
+        datasource.execute(""" 
+        DROP TABLE IF EXISTS $lcz_count_lod;
+        CREATE TABLE $lcz_count_lod as
+        SELECT  COUNT(*) FILTER (WHERE LCZ_PRIMARY IS NOT NULL) AS COUNT, LCZ_PRIMARY,
+        ID_ROW_LOD_${i}, ID_COL_LOD_${i},
+        CASE WHEN LCZ_PRIMARY=105 THEN 11
+        WHEN LCZ_PRIMARY=107 THEN 12
+        WHEN LCZ_PRIMARY=106 THEN 13
+        WHEN LCZ_PRIMARY= 101 THEN 14
+        WHEN LCZ_PRIMARY =102 THEN 15
+        WHEN LCZ_PRIMARY IN (103,104) THEN 16
+        ELSE LCZ_PRIMARY END AS weight_lcz,
+        FROM $grid_scaling_indices 
+        WHERE LCZ_PRIMARY IS NOT NULL 
+        GROUP BY ID_ROW_LOD_${i}, ID_COL_LOD_${i}, LCZ_PRIMARY;""".toString())
+
+        //Select the LCZ values according the maximum number of cells and the weight
+        //Note that we compute the number of cells for urban and cool LCZ
+        def lcz_count_lod_mode = postfix("lcz_count_lod_mode")
+        tablesToDrop << lcz_count_lod_mode
+        datasource.execute(""" 
+        CREATE INDEX ON $lcz_count_lod(ID_ROW_LOD_${i}, ID_COL_LOD_${i});
+        DROP TABLE IF EXISTS  $lcz_count_lod_mode;    
+        CREATE TABLE $lcz_count_lod_mode as
+        select distinct on (ID_ROW_LOD_${i}, ID_COL_LOD_${i}) *,
+        (select sum(count) from  $lcz_count_lod where   
+        LCZ_PRIMARY in (1,2,3,4,5,6,7,8,9,10,105) 
+        and ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i} and ID_COL_LOD_${i}= a.ID_COL_LOD_${i}) AS  LCZ_WARM_LOD_${i},
+        (select sum(count) from  $lcz_count_lod where  
+        LCZ_PRIMARY in (101,102,103,104,106,107) 
+        and ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i} and ID_COL_LOD_${i}= a.ID_COL_LOD_${i}) AS  LCZ_COOL_LOD_${i},
+        from $lcz_count_lod as a
+        order by count desc, ID_ROW_LOD_${i}, ID_COL_LOD_${i}, weight_lcz;""".toString())
+
+        //Find the 8 adjacent cells for the current level
+        def grid_lod_level_final = postfix("grid_lod_level_final")
+        tablesToDrop << grid_lod_level_final
+        datasource.execute("""
+        CREATE INDEX on $lcz_count_lod_mode(ID_ROW_LOD_${i}, ID_COL_LOD_${i});
+        DROP TABLE IF EXISTS $grid_lod_level_final;
+        CREATE TABLE $grid_lod_level_final as select * EXCEPT(LCZ_PRIMARY, COUNT, weight_lcz), LCZ_PRIMARY AS LCZ_PRIMARY_LOD_${i},
+        (SELECT LCZ_PRIMARY FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}+1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}) AS LCZ_PRIMARY_N_LOD_${i},
+        (SELECT LCZ_PRIMARY FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}+1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}+1) AS LCZ_PRIMARY_NE_LOD_${i},
+        (SELECT LCZ_PRIMARY FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i} AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}+1) AS LCZ_PRIMARY_E_LOD_${i},
+        (SELECT LCZ_PRIMARY FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}-1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}+1) AS LCZ_PRIMARY_SE_LOD_${i},
+        (SELECT LCZ_PRIMARY FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}-1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}) AS LCZ_PRIMARY_S_LOD_${i},
+        (SELECT LCZ_PRIMARY FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}-1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}-1) AS LCZ_PRIMARY_SW_LOD_${i},
+        (SELECT LCZ_PRIMARY FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i} AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}-1) AS LCZ_PRIMARY_W_LOD_${i},
+        (SELECT LCZ_PRIMARY FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}+1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}-1) AS LCZ_PRIMARY_NW_LOD_${i},
+
+        (SELECT LCZ_WARM_LOD_${i} FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}+1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}) AS LCZ_WARM_N_LOD_${i},
+        (SELECT LCZ_WARM_LOD_${i} FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}+1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}+1) AS LCZ_WARM_NE_LOD_${i},
+        (SELECT LCZ_WARM_LOD_${i} FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i} AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}+1) AS LCZ_WARM_E_LOD_${i},
+        (SELECT LCZ_WARM_LOD_${i} FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}-1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}+1) AS LCZ_WARM_SE_LOD_${i},
+        (SELECT LCZ_WARM_LOD_${i} FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}-1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}) AS LCZ_WARM_S_LOD_${i},
+        (SELECT LCZ_WARM_LOD_${i} FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}-1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}-1) AS LCZ_WARM_SW_LOD_${i},
+        (SELECT LCZ_WARM_LOD_${i} FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i} AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}-1) AS LCZ_WARM_W_LOD_${i},
+        (SELECT LCZ_WARM_LOD_${i} FROM $lcz_count_lod_mode WHERE ID_ROW_LOD_${i} = a.ID_ROW_LOD_${i}+1 AND ID_COL_LOD_${i}=a.ID_COL_LOD_${i}-1) AS LCZ_WARM_NW_LOD_${i},
+
+         FROM  $lcz_count_lod_mode as a;  """.toString())
+
+        tableLevelToJoin << grid_lod_level_final
+
+        //Join the final grid level with the original grid
+        def grid_level_join = postfix("grid_level_join")
+        datasource.execute("""
+        CREATE INDEX IF NOT EXISTS ${tableLevelToJoin}_idx ON $tableLevelToJoin (ID_ROW_LOD_${i}, ID_COL_LOD_${i});
+        create index on $grid_lod_level_final(ID_ROW_LOD_${i}, ID_COL_LOD_${i});
+        DROP TABLE IF EXISTS $grid_level_join;
+        CREATE TABLE $grid_level_join as 
+        select a.* EXCEPT(ID_ROW_LOD_${i}, ID_COL_LOD_${i}),  
+        b.* from $tableLevelToJoin as a,  $grid_lod_level_final as b 
+        where a.ID_ROW_LOD_${i} = b.ID_ROW_LOD_${i} and a.ID_COL_LOD_${i}= b.ID_COL_LOD_${i}
+        group by a.ID_ROW_LOD_${i}, a.ID_COL_LOD_${i} , a.id_grid;
+        """.toString())
+        tableLevelToJoin = grid_level_join
+    }
+    datasource.dropTable(tablesToDrop)
+    return tableLevelToJoin
+
+}
+
+
+/**
+ * Compute the distance from each grid cell to the edge of a polygon
+ *
+ * @param datasource a connection to the database
+ * @param input_polygons a set of polygons
+ * @param grid a regular grid
+ * @param id_grid name of the unique identifier column for the cells of the grid
+ * @author Erwan Bocher (CNRS)
+ * TODO : convert this method as a function table in H2GIS
+ */
+String gridDistances(JdbcDataSource datasource, String input_polygons, String grid, String id_grid) {
+    if (!input_polygons) {
+        error("The input polygons cannot be null or empty")
+        return
+    }
+    if (!grid) {
+        error("The grid cannot be null or empty")
+        return
+    }
+    if (!id_grid) {
+        error("Please set the column name identifier for the grid cells")
+        return
+    }
+    int epsg = datasource.getSrid(grid)
+    def outputTableName = postfix("grid_distances")
+
+    datasource.execute(""" DROP TABLE IF EXISTS $outputTableName;
+        CREATE TABLE $outputTableName (THE_GEOM GEOMETRY,ID INT, DISTANCE FLOAT);
+        """.toString())
+
+    datasource.createSpatialIndex(input_polygons)
+    datasource.createSpatialIndex(grid)
+
+    datasource.withBatch(100) { stmt ->
+        datasource.eachRow("SELECT the_geom from $input_polygons".toString()) { row ->
+            Geometry geom = row.the_geom
+            if (geom) {
+                IndexedFacetDistance indexedFacetDistance = new IndexedFacetDistance(geom)
+                datasource.eachRow("""SELECT the_geom, ${id_grid} as id from $grid 
+                where ST_GEOMFROMTEXT('${geom}',$epsg)  && the_geom and 
+            st_intersects(ST_GEOMFROMTEXT('${geom}',$epsg) , ST_POINTONSURFACE(the_geom))""".toString()) { cell ->
+                    Geometry cell_geom = cell.the_geom
+                    double distance = indexedFacetDistance.distance(cell_geom.getCentroid())
+                    stmt.addBatch "insert into $outputTableName values(ST_GEOMFROMTEXT('${cell_geom}',$epsg), ${cell.id},${distance})".toString()
+                }
+            }
+        }
+    }
+    return outputTableName
+}