This PR proposes new indicators to investigate the characteristics of…

… a territory: spatial distribution of objects, measurement of heterogeneity, neighborhood analyses, distribution of objects according to geographical directions, distances between objects...

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
 
@@ -93,4 +95,213 @@ 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 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, 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
+    }
+    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 = 0
+    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) 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 the_geom,
+    ID_GRID, ID_ROW, ID_COL, ${grid_levels_query.join(",")},
+    LCZ_PRIMARY,   
+    (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+1=a.ID_COL) 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())
+
+    def tablesToDrop =[]
+    def tableLevelToJoin = grid_scaling_indices
+    tablesToDrop<<grid_scaling_indices
+
+    //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_PRIMARY_URBAN_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_PRIMARY_COOL_LOD_${i},
+    from $lcz_count_lod as a
+    order by ID_ROW_LOD_${i}, ID_COL_LOD_${i}, weight_lcz desc, count asc;""".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}+1=a.ID_COL_LOD_${i}) 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} 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.the_geom;
+        """.toString())
+        tableLevelToJoin=grid_level_join
+    }
+
+    datasource.save(tableLevelToJoin, "/tmp/tableLevelToJoin.geojson", true)
+
+    for (int i in 1..nb_levels) {
+
+        def grid_lod = postfix("grid_lod")
+        tablesToDrop << grid_lod
+        datasource.execute("""
+    create index on $tableLevelToJoin(ID_ROW_LOD_${i}, ID_COL_LOD_${i});
+    DROP TABLE IF EXIsTS $grid_lod;
+    CREATE TABLE $grid_lod AS 
+    SELECT ST_UNION(ST_ACCUM(THE_GEOM)) AS THE_GEOM
+    from $tableLevelToJoin
+    group by ID_ROW_LOD_${i}, ID_COL_LOD_${i};
+    """.toString())
+
+        datasource.save(grid_lod, "/tmp/grid_lod_${i}.geojson", true)
+    }
+
+    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)
+ */
+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_distance")
+
+    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_geom.id},${distance})".toString()
+                }
+            }
+        }
+    }
+    return outputTableName
+}