Skip to content

Commit

Permalink
Merge pull request #937 from ebocher/master_sprawl
Browse files Browse the repository at this point in the history 
New indicators to investigate LCZ classification
  • Loading branch information
@ebocher
ebocher authored Mar 14, 2024
2 parents a2958d4 + 09ac654 commit 7a51ba3
Show file tree
Hide file tree
Showing 7 changed files with 658 additions and 20 deletions.
1 change: 0 additions & 1 deletion geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -971,7 +971,7 @@ String upperScaleAreaStatistics(JdbcDataSource datasource, String upperTableName
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);
Expand Down
230 changes: 227 additions & 3 deletions geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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 to select only one LCZ type
* 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
}
Loading

0 comments on commit 7a51ba3

Please sign in to comment.