diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy index 52bf6d5267..d45ba526a4 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy +++ b/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 diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy index 2bc991c0fd..88bfed584b 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy @@ -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); diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy index 6604ec7738..2ac325ad9d 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy +++ b/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 -} \ No newline at end of file +} + + +/** + * 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 +} diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy index f73a388d89..b5e29b27e2 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -62,7 +62,7 @@ import static org.h2gis.network.functions.ST_ConnectedComponents.getConnectedCom */ String createTSU(JdbcDataSource datasource, String zone, double area = 1f, String road, String rail, String vegetation, - String water, String sea_land_mask,String urban_areas, + String water, String sea_land_mask, String urban_areas, double surface_vegetation, double surface_hydro, double surface_urban_areas, String prefixName) { def BASE_NAME = "rsu" @@ -214,7 +214,7 @@ String prepareTSUData(JdbcDataSource datasource, String zone, String road, Strin } if (vegetation && datasource.hasTable(vegetation)) { debug "Preparing vegetation..." - if(datasource.getColumnNames(vegetation)) { + if (datasource.getColumnNames(vegetation)) { vegetation_tmp = postfix "vegetation_tmp" def vegetation_graph = postfix "vegetation_graph" def subGraphTableNodes = postfix vegetation_graph, "NODE_CC" @@ -270,7 +270,7 @@ String prepareTSUData(JdbcDataSource datasource, String zone, String road, Strin } if (water && datasource.hasTable(water)) { - if(datasource.getColumnNames(water).size()>0) { + if (datasource.getColumnNames(water).size() > 0) { //Extract water debug "Preparing hydrographic..." hydrographic_tmp = postfix "hydrographic_tmp" @@ -318,8 +318,8 @@ String prepareTSUData(JdbcDataSource datasource, String zone, String road, Strin } if (road && datasource.hasTable(road)) { - debug "Preparing road..." - if(datasource.getColumnNames(road).size()>0) { + debug "Preparing road..." + if (datasource.getColumnNames(road).size() > 0) { queryCreateOutputTable += [road_tmp: "(SELECT ST_ToMultiLine(THE_GEOM) FROM $road where (zindex=0 or crossing in ('bridge', 'crossing')) " + "and type not in ('track','service', 'path', 'cycleway', 'steps'))"] } @@ -327,8 +327,8 @@ String prepareTSUData(JdbcDataSource datasource, String zone, String road, Strin if (rail && datasource.hasTable(rail)) { debug "Preparing rail..." - if(datasource.getColumnNames(rail).size()>0){ - queryCreateOutputTable += [rail_tmp: "(SELECT ST_ToMultiLine(THE_GEOM) FROM $rail where (zindex=0 and usage='main') or (crossing = 'bridge' and usage='main'))"] + if (datasource.getColumnNames(rail).size() > 0) { + queryCreateOutputTable += [rail_tmp: "(SELECT ST_ToMultiLine(THE_GEOM) FROM $rail where (zindex=0 and usage='main') or (crossing = 'bridge' and usage='main'))"] } } @@ -605,4 +605,125 @@ String createGrid(JdbcDataSource datasource, Geometry geometry, double deltaX, d } debug "The table $outputTableName has been created" return outputTableName +} + + +/** + * This method allows to compute a sprawl areas layer from a regular grid that + * contains the fraction area of each LCZ type for each cell. + * + * A sprawl geometry represents a continous areas of urban LCZs (1 to 10 plus 105) + * It is important to note that pixels that do not have at least 1 urban neighbor are not kept. + * @param datasource connexion to the database + * @param grid_indicators a grid that contains the LCZ fractions + * @param distance value to erode (delete) small sprawl areas + * @author Erwan Bocher (CNRS) + */ +String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, float distance = 100) { + //We must compute the grid + if (!grid_indicators) { + error("No grid_indicators table to compute the sprawl areas layer") + return + } + if (distance < 0) { + error("Please set a distance greater or equal than 0") + return + } + if (datasource.getRowCount(grid_indicators) == 0) { + error("No grid cells to compute the sprawl areas layer") + return + } + def gridCols = datasource.getColumnNames(grid_indicators) + def lcz_columns_urban = ["LCZ_PRIMARY", "LCZ_WARM"] + def lcz_columns = gridCols.intersect(lcz_columns_urban) + if (lcz_columns.size()>0) { + def outputTableName = postfix("sprawl_areas") + if (distance == 0) { + datasource.execute("""DROP TABLE IF EXISTS $outputTableName; + create table $outputTableName as + select CAST((row_number() over()) as Integer) as id, st_removeholes(the_geom) as the_geom from ST_EXPLODE('( + select st_union(st_accum(the_geom)) as the_geom from + $grid_indicators where lcz_warm>=2 + and LCZ_PRIMARY NOT IN (101, 102,103,104,106, 107))')""".toString()) + return outputTableName + } else { + def tmp_sprawl = postfix("sprawl_tmp") + datasource.execute(""" + DROP TABLE IF EXISTS $tmp_sprawl, $outputTableName; + create table $tmp_sprawl as + select CAST((row_number() over()) as Integer) as id, st_removeholes(the_geom) as the_geom from ST_EXPLODE('( + select st_union(st_accum(the_geom)) as the_geom from + $grid_indicators where lcz_warm>=2 + and LCZ_PRIMARY NOT IN (101, 102,103,104,106, 107))') + where st_isempty(st_buffer(the_geom, -100)) =false""".toString()) + datasource.execute("""CREATE TABLE $outputTableName as SELECT CAST((row_number() over()) as Integer) as id, + the_geom + FROM + ST_EXPLODE('( + SELECT + st_removeholes(st_buffer(st_union(st_accum(st_buffer(st_removeholes(the_geom),$distance, ''quad_segs=2 endcap=flat + join=mitre mitre_limit=2''))), + -$distance, ''quad_segs=2 endcap=flat join=mitre mitre_limit=2'')) as the_geom + FROM ST_EXPLODE(''$tmp_sprawl'') )') ; + DROP TABLE IF EXISTS $tmp_sprawl; + """.toString()) + return outputTableName + } + } + error("No LCZ_PRIMARY column to compute the sprawl areas layer") + return +} + +/** + * This method is used to compute the difference between an input layer of polygons and the bounding box + * of the input layer. + * @param input_polygons a layer that contains polygons + * @author Erwan Bocher (CNRS) + */ +String inversePolygonsLayer(JdbcDataSource datasource, String input_polygons) { + def outputTableName = postfix("inverse_geometries") + def tmp_extent = postfix("tmp_extent") + datasource.execute("""DROP TABLE IF EXISTS $tmp_extent, $outputTableName; + CREATE TABLE $tmp_extent as SELECT ST_EXTENT(THE_GEOM) as the_geom FROM $input_polygons; + CREATE TABLE $outputTableName as + SELECT CAST((row_number() over()) as Integer) as id, the_geom + FROM + ST_EXPLODE('( + select st_difference(a.the_geom, st_accum(b.the_geom)) as the_geom from $tmp_extent as a, $input_polygons + as b where st_dimension(b.the_geom)=2)'); + DROP TABLE IF EXISTS $tmp_extent; + """.toString()) + return outputTableName +} + + + +/** + * This methods allows to extract the cool area geometries inside polygons + * A cool area is continous geometry defined by vegetation and water fractions. + * + * @author Erwan Bocher (CNRS) + */ +String extractCoolAreas(JdbcDataSource datasource, String grid_indicators,float distance = 100) { + if (!grid_indicators) { + error("No grid_indicators table to extract the cool areas layer") + return + } + def gridCols = datasource.getColumnNames(grid_indicators) + def lcz_columns_urban = ["LCZ_PRIMARY"] + def lcz_columns = gridCols.intersect(lcz_columns_urban) + + if (lcz_columns.size()>0) { + def outputTableName = postfix("cool_areas") + datasource.execute(""" + DROP TABLE IF EXISTS $outputTableName; + CREATE TABLE $outputTableName as SELECT CAST((row_number() over()) as Integer) as id, the_geom FROM ST_EXPLODE('( + SELECT ST_UNION(ST_ACCUM(a.THE_GEOM)) AS THE_GEOM FROM $grid_indicators as a + where + a.LCZ_PRIMARY in (101, 102, 103,104, 106, 107))') ${distance>0?" where st_isempty(st_buffer(the_geom, -$distance)) =false":""}; + """.toString()) + return outputTableName + } + error("No LCZ_PRIMARY column to extract the cool areas") + return } \ No newline at end of file diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy new file mode 100644 index 0000000000..a3ec4cc5c4 --- /dev/null +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy @@ -0,0 +1,125 @@ +package org.orbisgis.geoclimate.geoindicators + +import org.junit.jupiter.api.BeforeAll +import org.junit.jupiter.api.BeforeEach +import org.junit.jupiter.api.Disabled +import org.junit.jupiter.api.Test +import org.junit.jupiter.api.io.TempDir +import org.orbisgis.data.H2GIS +import org.orbisgis.data.POSTGIS +import org.orbisgis.geoclimate.Geoindicators + +import static org.junit.jupiter.api.Assertions.assertEquals +import static org.junit.jupiter.api.Assertions.assertTrue +import static org.orbisgis.data.H2GIS.open + +class GridIndicatorsTests { + + @TempDir + static File folder + + private static H2GIS h2GIS + + @BeforeAll + static void beforeAll() { + h2GIS = open(folder.getAbsolutePath() + File.separator + "gridIndicatorsTests") + } + + @BeforeEach + void beforeEach() { + } + + @Test + void multiscaleLCZGridTest() { + //Data for test + h2GIS.execute(""" + --Grid values + DROP TABLE IF EXISTS grid; + CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid, 104 AS lcz_primary FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + + --First cell lod_0 + UPDATE grid SET lcz_primary= 2 WHERE id_row = 2 AND id_col = 2; + + --Center cell lod_0 + UPDATE grid SET lcz_primary= 102 WHERE id_row = 5 AND id_col = 5; + UPDATE grid SET lcz_primary= 2 WHERE id_row = 6 AND id_col = 4; + UPDATE grid SET lcz_primary= 2 WHERE id_row = 6 AND id_col = 5; + UPDATE grid SET lcz_primary= 2 WHERE id_row = 6 AND id_col = 6; + UPDATE grid SET lcz_primary= 2 WHERE id_row = 5 AND id_col = 6; + + --Last cell lod_0 + UPDATE grid SET lcz_primary= 2 WHERE id_row = 8 AND id_col = 7; + UPDATE grid SET lcz_primary= 2 WHERE id_row = 8 AND id_col = 9; + UPDATE grid SET lcz_primary= 2 WHERE id_row = 7 AND id_col = 7; + UPDATE grid SET lcz_primary= 2 WHERE id_row = 7 AND id_col = 8; + UPDATE grid SET lcz_primary= 2 WHERE id_row = 7 AND id_col = 9; + """.toString()) + String grid_scale = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS, "grid","id_grid", 2) + + def values = h2GIS.firstRow("SELECT * EXCEPT(THE_GEOM) FROM $grid_scale WHERE id_row = 2 AND id_col = 2 ".toString()) + + def expectedValues = [ID_COL:2, ID_ROW:2, ID_GRID:10, LCZ_PRIMARY:2, LCZ_PRIMARY_N:104, LCZ_PRIMARY_NE:104, LCZ_PRIMARY_E:104, LCZ_PRIMARY_SE:104, LCZ_PRIMARY_S:104, LCZ_PRIMARY_SW:104, LCZ_PRIMARY_W:104, LCZ_PRIMARY_NW:104, LCZ_WARM:1, ID_ROW_LOD_1:1, ID_COL_LOD_1:0, LCZ_WARM_LOD_1:1, LCZ_COOL_LOD_1:8, LCZ_PRIMARY_LOD_1:104, LCZ_PRIMARY_N_LOD_1:104, LCZ_PRIMARY_NE_LOD_1:2, LCZ_PRIMARY_E_LOD_1:104, LCZ_PRIMARY_SE_LOD_1:null, LCZ_PRIMARY_S_LOD_1:null, LCZ_PRIMARY_SW_LOD_1:null, LCZ_PRIMARY_W_LOD_1:null, LCZ_PRIMARY_NW_LOD_1:null, LCZ_WARM_N_LOD_1:null, LCZ_WARM_NE_LOD_1:4, LCZ_WARM_E_LOD_1:null, LCZ_WARM_SE_LOD_1:null, LCZ_WARM_S_LOD_1:null, LCZ_WARM_SW_LOD_1:null, LCZ_WARM_W_LOD_1:null, LCZ_WARM_NW_LOD_1:null, ID_ROW_LOD_2:1, ID_COL_LOD_2:1, LCZ_WARM_LOD_2:10, LCZ_COOL_LOD_2:71, LCZ_PRIMARY_LOD_2:104, LCZ_PRIMARY_N_LOD_2:null, LCZ_PRIMARY_NE_LOD_2:null, LCZ_PRIMARY_E_LOD_2:null, LCZ_PRIMARY_SE_LOD_2:null, LCZ_PRIMARY_S_LOD_2:null, LCZ_PRIMARY_SW_LOD_2:null, LCZ_PRIMARY_W_LOD_2:null, LCZ_PRIMARY_NW_LOD_2:null, LCZ_WARM_N_LOD_2:null, LCZ_WARM_NE_LOD_2:null, LCZ_WARM_E_LOD_2:null, LCZ_WARM_SE_LOD_2:null, LCZ_WARM_S_LOD_2:null, LCZ_WARM_SW_LOD_2:null, LCZ_WARM_W_LOD_2:null, LCZ_WARM_NW_LOD_2:null] + + assertTrue(values == expectedValues) + + values = h2GIS.firstRow("SELECT * EXCEPT(THE_GEOM) FROM $grid_scale WHERE id_row = 5 AND id_col = 5 ".toString()) + + expectedValues = [ID_COL:5, ID_ROW:5, ID_GRID:40, LCZ_PRIMARY:102, LCZ_PRIMARY_N:2, LCZ_PRIMARY_NE:2, LCZ_PRIMARY_E:2, LCZ_PRIMARY_SE:104, LCZ_PRIMARY_S:104, LCZ_PRIMARY_SW:104, LCZ_PRIMARY_W:104, LCZ_PRIMARY_NW:2, LCZ_WARM:4, ID_ROW_LOD_1:2, ID_COL_LOD_1:1, LCZ_WARM_LOD_1:4, LCZ_COOL_LOD_1:5, LCZ_PRIMARY_LOD_1:2, LCZ_PRIMARY_N_LOD_1:104, LCZ_PRIMARY_NE_LOD_1:2, LCZ_PRIMARY_E_LOD_1:104, LCZ_PRIMARY_SE_LOD_1:104, LCZ_PRIMARY_S_LOD_1:104, LCZ_PRIMARY_SW_LOD_1:104, LCZ_PRIMARY_W_LOD_1:104, LCZ_PRIMARY_NW_LOD_1:104, LCZ_WARM_N_LOD_1:null, LCZ_WARM_NE_LOD_1:5, LCZ_WARM_E_LOD_1:null, LCZ_WARM_SE_LOD_1:null, LCZ_WARM_S_LOD_1:null, LCZ_WARM_SW_LOD_1:1, LCZ_WARM_W_LOD_1:null, LCZ_WARM_NW_LOD_1:null, ID_ROW_LOD_2:1, ID_COL_LOD_2:1, LCZ_WARM_LOD_2:10, LCZ_COOL_LOD_2:71, LCZ_PRIMARY_LOD_2:104, LCZ_PRIMARY_N_LOD_2:null, LCZ_PRIMARY_NE_LOD_2:null, LCZ_PRIMARY_E_LOD_2:null, LCZ_PRIMARY_SE_LOD_2:null, LCZ_PRIMARY_S_LOD_2:null, LCZ_PRIMARY_SW_LOD_2:null, LCZ_PRIMARY_W_LOD_2:null, LCZ_PRIMARY_NW_LOD_2:null, LCZ_WARM_N_LOD_2:null, LCZ_WARM_NE_LOD_2:null, LCZ_WARM_E_LOD_2:null, LCZ_WARM_SE_LOD_2:null, LCZ_WARM_S_LOD_2:null, LCZ_WARM_SW_LOD_2:null, LCZ_WARM_W_LOD_2:null, LCZ_WARM_NW_LOD_2:null] + + assertTrue(values == expectedValues) + + h2GIS.dropTable("grid") + } + + @Test + @Disabled + //Todo a test that shows how to create the a geom layer for each lod + void multiscaleLCZGridGeomTest() { + String grid_indicators = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Dijon/grid_indicators.geojson", true) + int nb_levels= 3 + String grid_scale = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS, grid_indicators,"id_grid", nb_levels) + for (int i in 1..nb_levels) { + def grid_lod = "grid_lod_$i" + h2GIS.execute(""" + create index on $grid_scale(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, MAX(LCZ_PRIMARY_LOD_${i}) AS LCZ_PRIMARY + from $grid_scale + group by ID_ROW_LOD_${i}, ID_COL_LOD_${i}; + """.toString()) + h2GIS.save(grid_lod, "/tmp/grid_lod_${i}.geojson", true) + h2GIS.dropTable(grid_lod) + } + } + + @Test + void gridDistancesTest1() { + //Data for test + h2GIS.execute(""" + --Grid values + DROP TABLE IF EXISTS grid, polygons; + CREATE TABLE grid AS SELECT * FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + CREATE TABLE polygons AS SELECT 'POLYGON ((4 4, 6 4, 6 6, 4 6, 4 4))'::GEOMETRY AS THE_GEOM ; + """.toString()) + + String grid_distances = Geoindicators.GridIndicators.gridDistances(h2GIS, "polygons","grid", "id") + assertEquals(4, h2GIS.firstRow("select count(*) as count from $grid_distances where distance =0.5".toString()).count) + } + + @Test + void gridDistancesTest2() { + //Data for test + h2GIS.execute(""" + --Grid values + DROP TABLE IF EXISTS grid, polygons; + CREATE TABLE grid AS SELECT * FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + CREATE TABLE polygons AS SELECT 'POLYGON ((2 2, 6 2, 6 6, 2 6, 2 2), (3 5, 5 5, 5 3, 3 3, 3 5))'::GEOMETRY AS THE_GEOM ; + """.toString()) + String grid_distances = Geoindicators.GridIndicators.gridDistances(h2GIS, "polygons","grid", "id") + assertEquals(12, h2GIS.firstRow("select count(*) as count from $grid_distances where distance =0.5".toString()).count) + } +} diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy index 2c1a2a52ee..2478dd81d3 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy @@ -21,9 +21,11 @@ package org.orbisgis.geoclimate.geoindicators import org.junit.jupiter.api.BeforeAll import org.junit.jupiter.api.BeforeEach +import org.junit.jupiter.api.Disabled import org.junit.jupiter.api.Test import org.junit.jupiter.api.condition.EnabledIfSystemProperty import org.junit.jupiter.api.io.TempDir +import org.locationtech.jts.geom.Geometry import org.locationtech.jts.io.WKTReader import org.orbisgis.data.POSTGIS import org.orbisgis.geoclimate.Geoindicators @@ -31,6 +33,8 @@ import static org.junit.jupiter.api.Assertions.assertEquals import static org.junit.jupiter.api.Assertions.assertNotNull import org.orbisgis.data.H2GIS +import static org.junit.jupiter.api.Assertions.assertTrue + class SpatialUnitsTests { @TempDir @@ -202,7 +206,6 @@ class SpatialUnitsTests { @EnabledIfSystemProperty(named = "test.h2gis", matches = "false") @Test void regularGridTestH2GIS() { - def wktReader = new WKTReader() def box = wktReader.read('POLYGON((-180 -80, 180 -80, 180 80, -180 80, -180 -80))') def gridP = Geoindicators.SpatialUnits.createGrid(h2GIS, box, 1, 1) @@ -225,4 +228,156 @@ class SpatialUnitsTests { def countRows = postGIS.firstRow "select count(*) as numberOfRows from $outputTable" assert 100 == countRows.numberOfRows } + + @Test + void sprawlAreasTest1() { + //Data for test + h2GIS.execute(""" + --Grid values + DROP TABLE IF EXISTS grid; + CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid, 104 AS LCZ_PRIMARY FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + --Center cell urban + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 6; + """.toString()) + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"grid","id_grid",1) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 0) + assertEquals(1, h2GIS.firstRow("select count(*) as count from $sprawl_areas".toString()).count) + assertEquals(5, h2GIS.firstRow("select st_area(the_geom) as area from $sprawl_areas".toString()).area, 0.0001) + } + + @Test + void sprawlAreasTest2() { + //Data for test + h2GIS.execute(""" + --Grid values + DROP TABLE IF EXISTS grid; + CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid, 104 AS LCZ_PRIMARY FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + """.toString()) + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"grid","id_grid",1) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 0) + assertEquals(1, h2GIS.firstRow("select count(*) as count from $sprawl_areas".toString()).count) + assertTrue(h2GIS.firstRow("select st_union(st_accum(the_geom)) as the_geom from $sprawl_areas".toString()).the_geom.isEmpty()) + } + + @Test + void sprawlAreasTest3() { + //Data for test + h2GIS.execute(""" + --Grid values + DROP TABLE IF EXISTS grid; + CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid,104 AS LCZ_PRIMARY FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 4 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 4 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 4 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 6; + """.toString()) + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"grid","id_grid",1) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 0) + assertEquals(1, h2GIS.firstRow("select count(*) as count from $sprawl_areas".toString()).count) + assertEquals(9, h2GIS.firstRow("select st_area(the_geom) as area from $sprawl_areas".toString()).area, 0.0001) + } + + @Test + void sprawlAreasTest4() { + //Data for test + h2GIS.execute(""" + --Grid values + DROP TABLE IF EXISTS grid; + CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid, 104 AS LCZ_PRIMARY FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 4 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 4 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 4 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 6; + + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 9 AND id_col = 9; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 1 AND id_col = 1; + """.toString()) + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"grid","id_grid",1) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 0) + h2GIS.save(sprawl_areas, "/tmp/sprawl_areas.geojson", true) + h2GIS.save(grid_scales, "/tmp/grid_scales.geojson", true) + assertEquals(1, h2GIS.firstRow("select count(*) as count from $sprawl_areas".toString()).count) + assertEquals(9, h2GIS.firstRow("select st_area(st_accum(the_geom)) as area from $sprawl_areas".toString()).area, 0.0001) + } + + @Test + void inverseGeometriesTest1() { + //Data for test + h2GIS.execute(""" + DROP TABLE IF EXISTS polygons; + CREATE TABLE polygons AS SELECT 'POLYGON ((160 290, 260 290, 260 190, 160 190, 160 290))'::GEOMETRY as the_geom; + """.toString()) + String inverse = Geoindicators.SpatialUnits.inversePolygonsLayer(h2GIS, "polygons") + assertEquals(1, h2GIS.firstRow("select count(*) as count from $inverse".toString()).count) + assertTrue(h2GIS.firstRow("select st_accum(the_geom) as the_geom from $inverse".toString()).the_geom.isEmpty()) + } + + @Test + void inverseGeometriesTest2() { + def wktReader = new WKTReader() + Geometry expectedGeom = wktReader.read("POLYGON ((160 190, 260 190, 260 290, 320 290, 320 150, 240 150, 240 80, 160 80, 160 190))") + //Data for test + h2GIS.execute(""" + DROP TABLE IF EXISTS polygons; + CREATE TABLE polygons AS SELECT 'MULTIPOLYGON (((160 290, 260 290, 260 190, 160 190, 160 290)), + ((240 150, 320 150, 320 80, 240 80, 240 150)))'::GEOMETRY as the_geom; + """.toString()) + String inverse = Geoindicators.SpatialUnits.inversePolygonsLayer(h2GIS, "polygons") + h2GIS.save(inverse, "/tmp/inverse.geojson", true) + assertEquals(1, h2GIS.firstRow("select count(*) as count from $inverse".toString()).count) + assertTrue(expectedGeom.equals(h2GIS.firstRow("select the_geom from $inverse".toString()).the_geom)) + } + + @Test + void inverseGeometriesTest3() { + def wktReader = new WKTReader() + Geometry expectedGeom = wktReader.read("MULTIPOLYGON (((160 190, 260 190, 260 290, 320 290, 320 150, 240 150, 240 80, 160 80, 160 190)), ((230 265, 230 230, 189 230, 189 265, 230 265)))") + //Data for test + h2GIS.execute(""" + DROP TABLE IF EXISTS polygons; + CREATE TABLE polygons AS SELECT 'MULTIPOLYGON (((160 290, 260 290, 260 190, 160 190, 160 290), + (189 265, 230 265, 230 230, 189 230, 189 265)), ((240 150, 320 150, 320 80, 240 80, 240 150)))'::GEOMETRY as the_geom; + """.toString()) + String inverse = Geoindicators.SpatialUnits.inversePolygonsLayer(h2GIS, "polygons") + assertEquals(2, h2GIS.firstRow("select count(*) as count from $inverse".toString()).count) + assertTrue(expectedGeom.equals(h2GIS.firstRow("select st_accum(the_geom) as the_geom from $inverse".toString()).the_geom)) + } + + @Test + @Disabled + void sprawlAreasTestIntegration() { + //Data for test + String path = "/home/ebocher/Autres/data/geoclimate/uhi_lcz/Angers/" + String data = h2GIS.load("${path}grid_indicators.geojson") + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,data,"id_grid", 1) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 100) + h2GIS.save(sprawl_areas, "/tmp/sprawl_areas_indic.geojson", true) + h2GIS.save(grid_scales, "/tmp/grid_indicators.geojson", true) + String distances = Geoindicators.GridIndicators.gridDistances(h2GIS, sprawl_areas, data, "id_grid") + h2GIS.save(distances, "/tmp/distances.geojson", true) + + //Method to compute the cool areas distances + String cool_areas = Geoindicators.SpatialUnits.extractCoolAreas(h2GIS, grid_scales) + h2GIS.save(cool_areas, "/tmp/cool_areas.geojson", true) + String inverse_cool_areas = Geoindicators.SpatialUnits.inversePolygonsLayer(h2GIS,cool_areas) + h2GIS.save(inverse_cool_areas, "/tmp/inverse_cool_areas.geojson", true) + distances = Geoindicators.GridIndicators.gridDistances(h2GIS, inverse_cool_areas, data, "id_grid") + h2GIS.save(distances, "/tmp/cool_inverse_distances.geojson", true) + } } \ No newline at end of file diff --git a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy index 7d32a78f8b..4b10af7d9c 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -649,10 +649,11 @@ class WorflowOSMTest extends WorkflowAbstractTest { File dirFile = new File(directory) dirFile.delete() dirFile.mkdir() - def location = "Redon" - //def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData("Redon") - // location = nominatim.bbox - location=[47.190062,-1.614389,47.196959,-1.602330] + def location = "Dijon" + def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(location) + def grid_size = 100 + location = nominatim.bbox + //location=[47.4, -4.8, 47.6, -4.6] def osm_parmeters = [ "description" : "Example of configuration file to run the OSM workflow and store the result in a folder", "geoclimatedb": [ @@ -677,8 +678,8 @@ class WorflowOSMTest extends WorkflowAbstractTest { "indicatorUse": ["LCZ", "UTRF", "TEB"] ],"grid_indicators": [ - "x_size": 100, - "y_size": 100, + "x_size": grid_size, + "y_size": grid_size, //"rowCol": true, "indicators": ["BUILDING_FRACTION","BUILDING_HEIGHT", "BUILDING_POP", "BUILDING_TYPE_FRACTION", @@ -689,13 +690,26 @@ class WorflowOSMTest extends WorkflowAbstractTest { "ASPECT_RATIO","SVF", "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"] ]/*, "worldpop_indicators": true, + "road_traffic" : true, "noise_indicators" : [ "ground_acoustic": true ]*/ ] ] - OSM.workflow(osm_parmeters) + Map results = OSM.workflow(osm_parmeters) + if(results) { + H2GIS h2gis = H2GIS.open("${directory + File.separator}geoclimate_test_integration;AUTO_SERVER=TRUE") + def tableNames = results.values() + def gridTable = tableNames.grid_indicators[0] + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2gis, gridTable,grid_size, grid_size*grid_size) + def folder_save =location in Collection ? location.join("_") : location + def path = directory + File.separator + "osm_$folder_save" + File.separator + path = "/tmp/" + h2gis.save(sprawl_areas, path + "sprawl_areas.fgb", true) + + h2gis.save(tableNames.rsu_lcz[0], path + "rsu_lcz.fgb", true) + } } @Disabled