From eb9de183a6d5061b3d31d9cb93c684ea9ee11ba9 Mon Sep 17 00:00:00 2001 From: ebocher Date: Fri, 2 Feb 2024 14:26:26 +0100 Subject: [PATCH 01/12] Start new LCZIndicators.groovy --- .../orbisgis/geoclimate/Geoindicators.groovy | 2 +- .../geoindicators/LCZIndicators.groovy | 275 ++++++++++++++++++ .../SprawlIndicatorsTests.groovy | 92 ++++++ 3 files changed, 368 insertions(+), 1 deletion(-) create mode 100644 geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy create mode 100644 geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy index 52bf6d5267..bc251cac0e 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy @@ -39,7 +39,7 @@ abstract class Geoindicators extends AbstractScript { public static PopulationIndicators = new PopulationIndicators() public static GridIndicators = new GridIndicators() public static NoiseIndicators = new NoiseIndicators() - + public static SprawlIndicators = new LCZIndicators() public static WorkflowUtilities = new WorkflowUtilities() //Cache the XStream models diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy new file mode 100644 index 0000000000..7f5a747f24 --- /dev/null +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy @@ -0,0 +1,275 @@ +/** + * GeoClimate is a geospatial processing toolbox for environmental and climate studies + * https://github.com/orbisgis/geoclimate. + * + * This code is part of the GeoClimate project. GeoClimate is free software; + * you can redistribute it and/or modify it under the terms of the GNU + * Lesser General Public License as published by the Free Software Foundation; + * version 3.0 of the License. + * + * GeoClimate is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License + * for more details . + * + * + * For more information, please consult: + * https://github.com/orbisgis/geoclimate + * + */ +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 + + +@BaseScript Geoindicators geoindicators + +/** + * This method allows to compute a sprawl areas layer. + * A sprawl geometry represents a continous areas of urban LCZs (1 to 10 plus 105) + * @author Erwan Bocher (CNRS) + */ +String compute_sprawl_areas(JdbcDataSource datasource, String grid_indicators, float fraction = 0.65, float distance=100) { + if (!grid_indicators) { + error("No grid_indicators table to compute the sprawl areas layer") + return + } + if (fraction <= 0) { + error("Please set a fraction greater than 0") + return + } + if (distance < 0) { + error("Please set a fraction greater or equal than 0") + return + } + if (datasource.getRowCount(grid_indicators) == 0) { + error("No geometries to compute the sprawl areas layer") + return + } + def gridCols = datasource.getColumnNames(grid_indicators) + + def lcz_columns_urban = ["LCZ_PRIMARY_1", "LCZ_PRIMARY_2", "LCZ_PRIMARY_3", "LCZ_PRIMARY_4", "LCZ_PRIMARY_5", "LCZ_PRIMARY_6", "LCZ_PRIMARY_7", + "LCZ_PRIMARY_8", "LCZ_PRIMARY_9", "LCZ_PRIMARY_10", "LCZ_PRIMARY_105"] + def lcz_fractions = gridCols.intersect(lcz_columns_urban) + + if (lcz_fractions.size() > 0) { + def outputTableName = postfix("sprawl_areas") + def tmp_sprawl = postfix("sprawl_tmp") + datasource.execute(""" + DROP TABLE IF EXISTS $tmp_sprawl, $outputTableName; + CREATE TABLE $tmp_sprawl as select st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat + join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat + join=mitre mitre_limit=2') as the_geom from + st_explode('( + SELECT ST_UNION(st_removeholes(ST_UNION(ST_ACCUM(the_geom)))) + AS THE_GEOM FROM $grid_indicators where ${lcz_fractions.join("+")} >= $fraction + )'); + 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_buffer(st_union(st_accum(st_buffer(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''))') where st_isempty(st_buffer(the_geom, -$distance))=false; + DROP TABLE IF EXISTS $tmp_sprawl; + """.toString()) + + return outputTableName + } + error("No LCZ_PRIMARY columns to compute the sprawl areas layer") + return +} + + +/** + * + * This method is compute the difference between an input layer of polygons and the bounding box + * of the input layer. + * This layer is used to perform the distances out of the sprawl areas + * @author Erwan Bocher (CNRS) + */ +String inverse_geometries(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)'); + DROP TABLE IF EXISTS $tmp_extent; + """.toString()) + return outputTableName +} + + +/** + * This methods allows to extract the cool area geometries. + * A cool area is continous geometry defined by vegetation and water fractions. + * + * @author Erwan Bocher (CNRS) + */ +String sprawl_cool_areas(JdbcDataSource datasource, String sprawl, String grid_indicators, float fraction = 0.4, + float distance =100 ) { + if (!grid_indicators) { + error("No grid_indicators table to compute the sprawl cool areas layer") + return + } + def gridCols = datasource.getColumnNames(grid_indicators) + + def lcz_columns_urban = ["LCZ_PRIMARY_101", "LCZ_PRIMARY_102", + "LCZ_PRIMARY_103", "LCZ_PRIMARY_104","LCZ_PRIMARY_107"] + def lcz_fractions = gridCols.intersect(lcz_columns_urban) + + if (lcz_fractions.size() > 0) { + datasource.createSpatialIndex(sprawl) + datasource.createSpatialIndex(grid_indicators) + 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, st_removeholes(the_geom) as the_geom FROM ST_EXPLODE('( + SELECT ST_UNION(ST_ACCUM(a.THE_GEOM)) AS THE_GEOM FROM $grid_indicators as a, $sprawl as b where + a.the_geom && b.the_geom and st_intersects(b.the_geom, ST_POINTONSURFACE(a.the_geom)) and + a.${lcz_fractions.join("+ a.")} >= $fraction + )') where st_isempty(st_buffer(st_convexhull(the_geom), -$distance))=false; + """.toString()) + return outputTableName + } + error("No LCZ_PRIMARY columns to compute the cool areas") + return +} + + +/** + * @author Erwan Bocher (CNRS) + */ +String inverse_cool_areas(JdbcDataSource datasource, String inverse_sprawl, String cool_areas) { + + def outputTableName = postfix("inverse_cool_areas") + def tmp_inverse = postfix("tmp_inverse") + datasource.execute(""" + DROP TABLE IF EXISTS $tmp_inverse; + CREATE TABLE $tmp_inverse as + SELECT st_union(st_accum(the_geom)) as the_geom + from (select the_geom from $inverse_sprawl union all select the_geom from $cool_areas) """.toString()) + + 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 $tmp_inverse; + CREATE TABLE $outputTableName as + SELECT CAST((row_number() over()) as Integer) as id, st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat + join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat + join=mitre mitre_limit=2') as the_geom + FROM + ST_EXPLODE('( + select st_difference(a.the_geom, st_accum(b.the_geom)) as the_geom from $tmp_extent as a, $tmp_inverse + as b)'); + DROP TABLE IF EXISTS $tmp_extent; + """.toString()) + + return outputTableName +} + +/** + * @author Erwan Bocher (CNRS) + */ +String grid_distances(JdbcDataSource datasource, String sprawl, String grid_indicators) { + int epsg = datasource.getSrid(grid_indicators) + def outputTableName = postfix("grid_sprawl_distance") + + datasource.execute(""" DROP TABLE IF EXISTS $outputTableName; + CREATE TABLE $outputTableName ( THE_GEOM GEOMETRY,ID INT, DISTANCE FLOAT); + """.toString()) + + datasource.createSpatialIndex(sprawl) + datasource.createSpatialIndex(grid_indicators) + + datasource.withBatch(100) { stmt -> + datasource.eachRow("SELECT the_geom from $sprawl".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_indicators + where ST_GEOMFROMTEXT('${geom}',$epsg) && the_geom and + st_intersects(ST_GEOMFROMTEXT('${geom}',$epsg) , ST_POINTONSURFACE(the_geom))""".toString()) { pointRow -> + Geometry cell = pointRow.the_geom + double distance = indexedFacetDistance.distance(cell.getCentroid()) + stmt.addBatch "insert into $outputTableName values(ST_GEOMFROMTEXT('${cell}',$epsg), ${pointRow.id},${distance})".toString() + } + } + } + } + return outputTableName +} + +/** + * @author Erwan Bocher (CNRS) + */ +String scaling_grid(JdbcDataSource datasource, String rsu_lcz) { + + datasource.execute(""" + DROP TABLE IF EXISTS GRID_100, GRID_300, GRID_900; + CREATE TABLE GRID_100 AS SELECT * FROM ST_MakeGrid('$rsu_lcz', 100,100); + CREATE TABLE GRID_300 AS SELECT * FROM ST_MakeGrid('$rsu_lcz', 300,300); + CREATE TABLE GRID_900 AS SELECT * FROM ST_MakeGrid('$rsu_lcz', 900,900); + """.toString()) + + datasource.save("GRID_100", "/tmp/grid_100.geojson", true) + datasource.save("GRID_300", "/tmp/grid_300.geojson", true) + datasource.save("GRID_900", "/tmp/grid_900.geojson", true) + + def values =[:] + int row_index= 0 + datasource.eachRow("select * from grid_100 order by id_row"){row -> + if(row_index<=3){ + row_index++ + } + else { + row_index=0 + } + } + return null + +} + + +/** + * @author Erwan Bocher (CNRS) + */ +String buffering_cool_areas(JdbcDataSource datasource, String rsu_lcz, String grid_indicators) { + + float distance = 10 + def outputTable = postfix("buffering_cool_areas") + def merging_lcz = postfix("merging_lcz") + datasource.execute(""" +DROP TABLE IF EXISTS $merging_lcz; +CREATE TABLE $merging_lcz as +select CAST((row_number() over()) as Integer) as id, the_geom from ST_EXPLODE('( +select st_union(st_accum(st_buffer(THE_GEOM, 0.1))) as the_geom FROM $rsu_lcz where LCZ_PRIMARY IN (101, 102, 103,104, 107))') +""".toString()) + + distance = 100 + + datasource.execute(""" + DROP TABLE IF EXISTS $outputTable; + CREATE TABLE $outputTable as +select st_buffer(the_geom, -$distance,'quad_segs=2 endcap=flat + join=mitre mitre_limit=2') as the_geom from ST_EXPLODE('( + select st_union(st_accum(st_buffer(the_geom, $distance, ''quad_segs=2 endcap=flat join=mitre mitre_limit=2''))) as the_geom from $merging_lcz)') as foo""".toString()) + + def sprawl_layer = postfix("sprawl_layer") + datasource.execute(""" + DROP TABLE IF EXISTS $sprawl_layer; + CREATE TABLE $sprawl_layer as + + """.toString()) + + return outputTable +} \ No newline at end of file diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy new file mode 100644 index 0000000000..da6aa3c52d --- /dev/null +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy @@ -0,0 +1,92 @@ +/** + * GeoClimate is a geospatial processing toolbox for environmental and climate studies + * https://github.com/orbisgis/geoclimate. + * + * This code is part of the GeoClimate project. GeoClimate is free software; + * you can redistribute it and/or modify it under the terms of the GNU + * Lesser General Public License as published by the Free Software Foundation; + * version 3.0 of the License. + * + * GeoClimate is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License + * for more details . + * + * + * For more information, please consult: + * https://github.com/orbisgis/geoclimate + * + */ +package org.orbisgis.geoclimate.geoindicators + +import org.junit.jupiter.api.BeforeAll +import org.junit.jupiter.api.Test +import org.junit.jupiter.api.io.TempDir +import org.orbisgis.data.H2GIS +import org.orbisgis.geoclimate.Geoindicators + +import static org.orbisgis.data.H2GIS.open + +class SprawlIndicatorsTests { + + @TempDir + static File folder + private static H2GIS h2GIS + + @BeforeAll + static void beforeAll() { + h2GIS = open(folder.getAbsolutePath() + File.separator + "sprawlindicators") + } + + @Test + void debug_tests() { + + /*String rsu_lcz = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Angers/rsu_lcz.geojson", true) + def buffering_cool_areas = Geoindicators.SprawlIndicators.buffering_cool_areas(h2GIS, rsu_lcz, "") + + h2GIS.save(buffering_cool_areas, "/tmp/buffering_cool_areas.geojson", true)*/ + + + + String grid_indicators = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Lorient/grid_indicators.geojson", true) + + String scale_grid = Geoindicators.SprawlIndicators.scaling_grid(h2GIS, grid_indicators) + + // h2GIS.save(scale_grid, "/tmp/scale_grid.csv", true) + + /* + def sprawlLayer = Geoindicators.SprawlIndicators.compute_sprawl_areas(h2GIS, grid_indicators) + + h2GIS.save(sprawlLayer, "/tmp/sprawl_areas.geojson", true) + + String grid_distances = Geoindicators.SprawlIndicators.grid_distances(h2GIS, sprawlLayer, grid_indicators) + + h2GIS.save(grid_distances, "/tmp/sprawl_areas_distances.geojson", true) + + + String sprawlLayerInverse = Geoindicators.SprawlIndicators.inverse_geometries(h2GIS, sprawlLayer) + + h2GIS.save(sprawlLayerInverse, "/tmp/spraw_areas_inverse.geojson", true) + + + String grid_distances_out = Geoindicators.SprawlIndicators.grid_distances(h2GIS, sprawlLayerInverse, grid_indicators) + + h2GIS.save(grid_distances_out, "/tmp/sprawl_areas_distances_out.geojson", true) + + def cool_areasLayer = Geoindicators.SprawlIndicators.sprawl_cool_areas(h2GIS, sprawlLayer, grid_indicators) + + h2GIS.save(cool_areasLayer, "/tmp/cool_areas.geojson", true) + + String cool_areas_inverse = Geoindicators.SprawlIndicators.inverse_cool_areas(h2GIS,sprawlLayerInverse, cool_areasLayer) + + h2GIS.save(cool_areas_inverse, "/tmp/cool_areas_inverse.geojson", true) + + + String cool_areas_distances = Geoindicators.SprawlIndicators.grid_distances(h2GIS, cool_areas_inverse, grid_indicators) + + h2GIS.save(cool_areas_distances, "/tmp/cool_areas_distances.geojson", true)*/ + + } + + +} From 3728ad8725795b7bfe7de1283c66188b889aa617 Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 6 Feb 2024 14:05:08 +0100 Subject: [PATCH 02/12] Prepare new indicators --- .../geoindicators/LCZIndicators.groovy | 68 +++++++++++++------ .../SprawlIndicatorsTests.groovy | 21 ++++-- 2 files changed, 65 insertions(+), 24 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy index 7f5a747f24..fabec78ab8 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy @@ -212,29 +212,59 @@ String grid_distances(JdbcDataSource datasource, String sprawl, String grid_indi /** * @author Erwan Bocher (CNRS) */ -String scaling_grid(JdbcDataSource datasource, String rsu_lcz) { +String scaling_grid(JdbcDataSource datasource, String grid_indicators) { - datasource.execute(""" - DROP TABLE IF EXISTS GRID_100, GRID_300, GRID_900; - CREATE TABLE GRID_100 AS SELECT * FROM ST_MakeGrid('$rsu_lcz', 100,100); - CREATE TABLE GRID_300 AS SELECT * FROM ST_MakeGrid('$rsu_lcz', 300,300); - CREATE TABLE GRID_900 AS SELECT * FROM ST_MakeGrid('$rsu_lcz', 900,900); + if (!grid_indicators) { + error("No grid_indicators table to aggregate the LCZ fraction") + return + } + def gridCols = datasource.getColumnNames(grid_indicators) + + def lcz_columns_cool_areas = ["LCZ_PRIMARY_101", "LCZ_PRIMARY_102", + "LCZ_PRIMARY_103", "LCZ_PRIMARY_104","LCZ_PRIMARY_107"] + def lcz_fractions = gridCols.intersect(lcz_columns_cool_areas) + + def grid_scaling_indices = postfix("grid_scaling_indices") + def grid_fraction_lod_1 = postfix("grid_fraction_lod_1") + def grid_fraction_lod_2 = postfix("grid_fraction_lod_2") + datasource.execute("""DROP TABLE IF EXISTS $grid_scaling_indices,$grid_fraction_lod_1,$grid_fraction_lod_2; + CREATE TABLE $grid_scaling_indices as SELECT the_geom, + ID_GRID, ID_COL as ID_COL_LOD_0, ID_ROW as ID_ROW_LOD_0, (CAST (ID_COL/3 AS INT)+1) AS ID_COL_LOD_1, + (CAST (ID_ROW/3 AS INT)+1) AS ID_ROW_LOD_1, + (CAST (ID_COL/9 AS INT)+1) AS ID_COL_LOD_2, + (CAST (ID_ROW/9 AS INT)+1) AS ID_ROW_LOD_2, + ${lcz_fractions.join(",")}, ${lcz_fractions.join("+")} AS SUM_LCZ_COOL FROM $grid_indicators; + CREATE TABLE $grid_fraction_lod_1 as SELECT ID_COL_LOD_1, ID_ROW_LOD_1, SUM(SUM_LCZ_COOL) AS LOD_1_SUM_LCZ_COOL, + ST_UNION(ST_ACCUM(the_geom)) as the_geom, count(*) as count + FROM $grid_scaling_indices GROUP BY ID_COL_LOD_1, ID_ROW_LOD_1; + CREATE TABLE $grid_fraction_lod_2 as SELECT ID_COL_LOD_2, ID_ROW_LOD_2, SUM(SUM_LCZ_COOL) AS LOD_2_SUM_LCZ_COOL, + ST_UNION(ST_ACCUM(the_geom)) as the_geom,count(*) as count + FROM $grid_scaling_indices GROUP BY ID_COL_LOD_2, ID_ROW_LOD_2; """.toString()) + datasource.createIndex(grid_scaling_indices, "ID_COL_LOD_1") + datasource.createIndex(grid_scaling_indices, "ID_ROW_LOD_1") + + datasource.save(grid_fraction_lod_1, "/tmp/grid_lod1.geojson", true) + datasource.save(grid_fraction_lod_2, "/tmp/grid_lod2.geojson", true) + + def grid_fractions = postfix("grid_fractions") + + //"neighbors" + + /*int level =1 + //Collect neighbors + + datasource.execute( """ +DROP TABLE IF EXISTS neighbors; +create table neighbors as + SELECT the_geom, LOD_${level}_SUM_LCZ_COOL as LOD_${level}_LCZ_COOL_2 from $grid_fraction_lod_1 + where ID_COL_LOD_$level=ID_COL_LOD_$level AND ID_ROW_LOD_$level=ID_ROW_LOD_$level-1""".toString()) + + + datasource.save("neighbors", "/tmp/neighbors.geojson", true)*/ + - datasource.save("GRID_100", "/tmp/grid_100.geojson", true) - datasource.save("GRID_300", "/tmp/grid_300.geojson", true) - datasource.save("GRID_900", "/tmp/grid_900.geojson", true) - def values =[:] - int row_index= 0 - datasource.eachRow("select * from grid_100 order by id_row"){row -> - if(row_index<=3){ - row_index++ - } - else { - row_index=0 - } - } return null } diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy index da6aa3c52d..42643a0b8f 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy @@ -47,14 +47,25 @@ class SprawlIndicatorsTests { h2GIS.save(buffering_cool_areas, "/tmp/buffering_cool_areas.geojson", true)*/ + String grid_indicators = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Dijon/grid_indicators.geojson", true) - String grid_indicators = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Lorient/grid_indicators.geojson", true) - - String scale_grid = Geoindicators.SprawlIndicators.scaling_grid(h2GIS, grid_indicators) + //String scale_grid = Geoindicators.SprawlIndicators.scaling_grid(h2GIS, grid_indicators) // h2GIS.save(scale_grid, "/tmp/scale_grid.csv", true) - /* + //GRID SCALES + + h2GIS.execute("""DROP TABLE IF EXISTS lod_1,lod_2, lod_3; + CREATE TABLE lod_1 as select * from ST_MakeGrid('$grid_indicators', 300,300); + CREATE TABLE lod_2 as select * from ST_MakeGrid('$grid_indicators', 900,900); + CREATE TABLE lod_3 as select * from ST_MakeGrid('$grid_indicators', 2700,2700);""".toString()) + + h2GIS.save("lod_1", "/tmp/lod1.geojson", true) + h2GIS.save("lod_2", "/tmp/lod2.geojson", true) + h2GIS.save("lod_3", "/tmp/lod3.geojson", true) + + + def sprawlLayer = Geoindicators.SprawlIndicators.compute_sprawl_areas(h2GIS, grid_indicators) h2GIS.save(sprawlLayer, "/tmp/sprawl_areas.geojson", true) @@ -84,7 +95,7 @@ class SprawlIndicatorsTests { String cool_areas_distances = Geoindicators.SprawlIndicators.grid_distances(h2GIS, cool_areas_inverse, grid_indicators) - h2GIS.save(cool_areas_distances, "/tmp/cool_areas_distances.geojson", true)*/ + h2GIS.save(cool_areas_distances, "/tmp/cool_areas_distances.geojson", true) } From 88d97192680725cb8018ac9a2c489576927c3f10 Mon Sep 17 00:00:00 2001 From: ebocher Date: Fri, 16 Feb 2024 17:09:26 +0100 Subject: [PATCH 03/12] Scaling grid indicators --- .../orbisgis/geoclimate/Geoindicators.groovy | 2 +- .../geoindicators/LCZIndicators.groovy | 174 +++++++++++++----- .../SprawlIndicatorsTests.groovy | 20 +- 3 files changed, 136 insertions(+), 60 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy index bc251cac0e..567cc960c9 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy @@ -39,7 +39,7 @@ abstract class Geoindicators extends AbstractScript { public static PopulationIndicators = new PopulationIndicators() public static GridIndicators = new GridIndicators() public static NoiseIndicators = new NoiseIndicators() - public static SprawlIndicators = new LCZIndicators() + public static LCZIndicators = new LCZIndicators() public static WorkflowUtilities = new WorkflowUtilities() //Cache the XStream models diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy index fabec78ab8..b415c95213 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy @@ -32,8 +32,9 @@ import org.orbisgis.geoclimate.Geoindicators * This method allows to compute a sprawl areas layer. * A sprawl geometry represents a continous areas of urban LCZs (1 to 10 plus 105) * @author Erwan Bocher (CNRS) + * TODO: SPATIAL UNITS */ -String compute_sprawl_areas(JdbcDataSource datasource, String grid_indicators, float fraction = 0.65, float distance=100) { +String compute_sprawl_areas(JdbcDataSource datasource, String grid_indicators, float fraction = 0.65, float distance = 100) { if (!grid_indicators) { error("No grid_indicators table to compute the sprawl areas layer") return @@ -92,6 +93,7 @@ String compute_sprawl_areas(JdbcDataSource datasource, String grid_indicators, f * of the input layer. * This layer is used to perform the distances out of the sprawl areas * @author Erwan Bocher (CNRS) + * TODO : SPATIAL UNITS */ String inverse_geometries(JdbcDataSource datasource, String input_polygons) { def outputTableName = postfix("inverse_geometries") @@ -115,9 +117,10 @@ String inverse_geometries(JdbcDataSource datasource, String input_polygons) { * A cool area is continous geometry defined by vegetation and water fractions. * * @author Erwan Bocher (CNRS) + * TODO: SPATIAL UNITS */ String sprawl_cool_areas(JdbcDataSource datasource, String sprawl, String grid_indicators, float fraction = 0.4, - float distance =100 ) { + float distance = 100) { if (!grid_indicators) { error("No grid_indicators table to compute the sprawl cool areas layer") return @@ -125,7 +128,7 @@ String sprawl_cool_areas(JdbcDataSource datasource, String sprawl, String grid_i def gridCols = datasource.getColumnNames(grid_indicators) def lcz_columns_urban = ["LCZ_PRIMARY_101", "LCZ_PRIMARY_102", - "LCZ_PRIMARY_103", "LCZ_PRIMARY_104","LCZ_PRIMARY_107"] + "LCZ_PRIMARY_103", "LCZ_PRIMARY_104", "LCZ_PRIMARY_107"] def lcz_fractions = gridCols.intersect(lcz_columns_urban) if (lcz_fractions.size() > 0) { @@ -149,6 +152,8 @@ String sprawl_cool_areas(JdbcDataSource datasource, String sprawl, String grid_i /** * @author Erwan Bocher (CNRS) + * + * TODO : spatial units */ String inverse_cool_areas(JdbcDataSource datasource, String inverse_sprawl, String cool_areas) { @@ -179,6 +184,8 @@ String inverse_cool_areas(JdbcDataSource datasource, String inverse_sprawl, Stri /** * @author Erwan Bocher (CNRS) + * + * TODO : GRID INDICATORS */ String grid_distances(JdbcDataSource datasource, String sprawl, String grid_indicators) { int epsg = datasource.getSrid(grid_indicators) @@ -211,59 +218,132 @@ String grid_distances(JdbcDataSource datasource, String sprawl, String grid_indi /** * @author Erwan Bocher (CNRS) + * + * TODO : grid indicators */ -String scaling_grid(JdbcDataSource datasource, String grid_indicators) { +String scaling_grid(JdbcDataSource datasource, String grid_indicators, int nb_levels = 2) { if (!grid_indicators) { error("No grid_indicators table to aggregate the LCZ fraction") return } - def gridCols = datasource.getColumnNames(grid_indicators) - - def lcz_columns_cool_areas = ["LCZ_PRIMARY_101", "LCZ_PRIMARY_102", - "LCZ_PRIMARY_103", "LCZ_PRIMARY_104","LCZ_PRIMARY_107"] - def lcz_fractions = gridCols.intersect(lcz_columns_cool_areas) - + if (nb_levels <= 0 || nb_levels >= 10) { + error("The number of levels to aggregate the LCZ fractions 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_fraction_lod_1 = postfix("grid_fraction_lod_1") - def grid_fraction_lod_2 = postfix("grid_fraction_lod_2") - datasource.execute("""DROP TABLE IF EXISTS $grid_scaling_indices,$grid_fraction_lod_1,$grid_fraction_lod_2; + 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 indice 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_COL as ID_COL_LOD_0, ID_ROW as ID_ROW_LOD_0, (CAST (ID_COL/3 AS INT)+1) AS ID_COL_LOD_1, - (CAST (ID_ROW/3 AS INT)+1) AS ID_ROW_LOD_1, - (CAST (ID_COL/9 AS INT)+1) AS ID_COL_LOD_2, - (CAST (ID_ROW/9 AS INT)+1) AS ID_ROW_LOD_2, - ${lcz_fractions.join(",")}, ${lcz_fractions.join("+")} AS SUM_LCZ_COOL FROM $grid_indicators; - CREATE TABLE $grid_fraction_lod_1 as SELECT ID_COL_LOD_1, ID_ROW_LOD_1, SUM(SUM_LCZ_COOL) AS LOD_1_SUM_LCZ_COOL, - ST_UNION(ST_ACCUM(the_geom)) as the_geom, count(*) as count - FROM $grid_scaling_indices GROUP BY ID_COL_LOD_1, ID_ROW_LOD_1; - CREATE TABLE $grid_fraction_lod_2 as SELECT ID_COL_LOD_2, ID_ROW_LOD_2, SUM(SUM_LCZ_COOL) AS LOD_2_SUM_LCZ_COOL, - ST_UNION(ST_ACCUM(the_geom)) as the_geom,count(*) as count - FROM $grid_scaling_indices GROUP BY ID_COL_LOD_2, ID_ROW_LOD_2; + ID_GRID, ID_ROW as ID_ROW_LOD_0, ID_COL as ID_COL_LOD_0, ${grid_levels_query.join(",")}, + LCZ_PRIMARY AS LCZ_PRIMARY_LOD_0, + (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL=a.ID_COL) AS LCZ_PRIMARY_N_LOD_0, + (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL+1=a.ID_COL) AS LCZ_PRIMARY_NE_LOD_0, + (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW AND ID_COL=a.ID_COL+1) AS LCZ_PRIMARY_E_LOD_0, + (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL+1) AS LCZ_PRIMARY_SE_LOD_0, + (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL) AS LCZ_PRIMARY_S_LOD_0, + (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_SW_LOD_0, + (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_W_LOD_0, + (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_NW_LOD_0 FROM + $grid_indicators as a; """.toString()) + + datasource.save(grid_scaling_indices, "/tmp/grid_lod_0.geojson", true) + + def tablesToDrop =[] + def tablesLevelToJoin =[] + tablesToDrop< Date: Thu, 22 Feb 2024 09:49:43 +0100 Subject: [PATCH 04/12] 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/GridIndicators.groovy | 213 ++++++++++++++- .../geoindicators/LCZIndicators.groovy | 251 ------------------ .../geoindicators/SpatialUnits.groovy | 86 ++++++ .../geoindicators/GridIndicatorsTests.groovy | 72 +++++ .../geoindicators/SpatialUnitsTests.groovy | 10 +- .../SprawlIndicatorsTests.groovy | 10 +- 6 files changed, 382 insertions(+), 260 deletions(-) create mode 100644 geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy 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..2f40d58eb8 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 @@ -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 -} \ 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 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< + 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 +} diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy index b415c95213..1d383fbc4b 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy @@ -20,96 +20,12 @@ 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 @BaseScript Geoindicators geoindicators -/** - * This method allows to compute a sprawl areas layer. - * A sprawl geometry represents a continous areas of urban LCZs (1 to 10 plus 105) - * @author Erwan Bocher (CNRS) - * TODO: SPATIAL UNITS - */ -String compute_sprawl_areas(JdbcDataSource datasource, String grid_indicators, float fraction = 0.65, float distance = 100) { - if (!grid_indicators) { - error("No grid_indicators table to compute the sprawl areas layer") - return - } - if (fraction <= 0) { - error("Please set a fraction greater than 0") - return - } - if (distance < 0) { - error("Please set a fraction greater or equal than 0") - return - } - if (datasource.getRowCount(grid_indicators) == 0) { - error("No geometries to compute the sprawl areas layer") - return - } - def gridCols = datasource.getColumnNames(grid_indicators) - - def lcz_columns_urban = ["LCZ_PRIMARY_1", "LCZ_PRIMARY_2", "LCZ_PRIMARY_3", "LCZ_PRIMARY_4", "LCZ_PRIMARY_5", "LCZ_PRIMARY_6", "LCZ_PRIMARY_7", - "LCZ_PRIMARY_8", "LCZ_PRIMARY_9", "LCZ_PRIMARY_10", "LCZ_PRIMARY_105"] - def lcz_fractions = gridCols.intersect(lcz_columns_urban) - - if (lcz_fractions.size() > 0) { - def outputTableName = postfix("sprawl_areas") - def tmp_sprawl = postfix("sprawl_tmp") - datasource.execute(""" - DROP TABLE IF EXISTS $tmp_sprawl, $outputTableName; - CREATE TABLE $tmp_sprawl as select st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat - join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat - join=mitre mitre_limit=2') as the_geom from - st_explode('( - SELECT ST_UNION(st_removeholes(ST_UNION(ST_ACCUM(the_geom)))) - AS THE_GEOM FROM $grid_indicators where ${lcz_fractions.join("+")} >= $fraction - )'); - 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_buffer(st_union(st_accum(st_buffer(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''))') where st_isempty(st_buffer(the_geom, -$distance))=false; - DROP TABLE IF EXISTS $tmp_sprawl; - """.toString()) - - return outputTableName - } - error("No LCZ_PRIMARY columns to compute the sprawl areas layer") - return -} - - -/** - * - * This method is compute the difference between an input layer of polygons and the bounding box - * of the input layer. - * This layer is used to perform the distances out of the sprawl areas - * @author Erwan Bocher (CNRS) - * TODO : SPATIAL UNITS - */ -String inverse_geometries(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)'); - DROP TABLE IF EXISTS $tmp_extent; - """.toString()) - return outputTableName -} /** @@ -182,173 +98,6 @@ String inverse_cool_areas(JdbcDataSource datasource, String inverse_sprawl, Stri return outputTableName } -/** - * @author Erwan Bocher (CNRS) - * - * TODO : GRID INDICATORS - */ -String grid_distances(JdbcDataSource datasource, String sprawl, String grid_indicators) { - int epsg = datasource.getSrid(grid_indicators) - def outputTableName = postfix("grid_sprawl_distance") - - datasource.execute(""" DROP TABLE IF EXISTS $outputTableName; - CREATE TABLE $outputTableName ( THE_GEOM GEOMETRY,ID INT, DISTANCE FLOAT); - """.toString()) - - datasource.createSpatialIndex(sprawl) - datasource.createSpatialIndex(grid_indicators) - - datasource.withBatch(100) { stmt -> - datasource.eachRow("SELECT the_geom from $sprawl".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_indicators - where ST_GEOMFROMTEXT('${geom}',$epsg) && the_geom and - st_intersects(ST_GEOMFROMTEXT('${geom}',$epsg) , ST_POINTONSURFACE(the_geom))""".toString()) { pointRow -> - Geometry cell = pointRow.the_geom - double distance = indexedFacetDistance.distance(cell.getCentroid()) - stmt.addBatch "insert into $outputTableName values(ST_GEOMFROMTEXT('${cell}',$epsg), ${pointRow.id},${distance})".toString() - } - } - } - } - return outputTableName -} - -/** - * @author Erwan Bocher (CNRS) - * - * TODO : grid indicators - */ -String scaling_grid(JdbcDataSource datasource, String grid_indicators, int nb_levels = 2) { - - if (!grid_indicators) { - error("No grid_indicators table to aggregate the LCZ fraction") - return - } - if (nb_levels <= 0 || nb_levels >= 10) { - error("The number of levels to aggregate the LCZ fractions 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 indice 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 as ID_ROW_LOD_0, ID_COL as ID_COL_LOD_0, ${grid_levels_query.join(",")}, - LCZ_PRIMARY AS LCZ_PRIMARY_LOD_0, - (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL=a.ID_COL) AS LCZ_PRIMARY_N_LOD_0, - (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL+1=a.ID_COL) AS LCZ_PRIMARY_NE_LOD_0, - (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW AND ID_COL=a.ID_COL+1) AS LCZ_PRIMARY_E_LOD_0, - (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL+1) AS LCZ_PRIMARY_SE_LOD_0, - (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL) AS LCZ_PRIMARY_S_LOD_0, - (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW-1 AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_SW_LOD_0, - (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_W_LOD_0, - (SELECT LCZ_PRIMARY FROM $grid_indicators WHERE ID_ROW = a.ID_ROW+1 AND ID_COL=a.ID_COL-1) AS LCZ_PRIMARY_NW_LOD_0 FROM - $grid_indicators as a; """.toString()) - - datasource.save(grid_scaling_indices, "/tmp/grid_lod_0.geojson", true) - - def tablesToDrop =[] - def tablesLevelToJoin =[] - tablesToDrop< 0) { + def outputTableName = postfix("sprawl_areas") + def tmp_sprawl = postfix("sprawl_tmp") + datasource.execute(""" + DROP TABLE IF EXISTS $tmp_sprawl, $outputTableName; + CREATE TABLE $tmp_sprawl as select st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat + join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat + join=mitre mitre_limit=2') as the_geom from + st_explode('( + SELECT ST_UNION(st_removeholes(ST_UNION(ST_ACCUM(the_geom)))) + AS THE_GEOM FROM $grid_indicators where ${lcz_fractions.join("+")} >= $fraction + )'); + 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_buffer(st_union(st_accum(st_buffer(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''))') where st_isempty(st_buffer(the_geom, -$distance))=false; + DROP TABLE IF EXISTS $tmp_sprawl; + """.toString()) + + return outputTableName + } + error("No LCZ fractions columns 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 inverse_geometries(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 } \ 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..fb2e9544d3 --- /dev/null +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy @@ -0,0 +1,72 @@ +package org.orbisgis.geoclimate.geoindicators + +import org.junit.jupiter.api.BeforeAll +import org.junit.jupiter.api.BeforeEach +import org.junit.jupiter.api.Test +import org.junit.jupiter.api.io.TempDir +import org.orbisgis.data.H2GIS +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", 1) + + def values = h2GIS.firstRow("SELECT * EXCEPT(THE_GEOM) FROM $grid_scale WHERE id_row = 2 AND id_col = 2 ".toString()) + + def expectedValues=[ID_GRID:10, ID_ROW:2, ID_COL:2, 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, ID_ROW_LOD_1:1, ID_COL_LOD_1:0, LCZ_PRIMARY_URBAN_LOD_1:1, LCZ_PRIMARY_COOL_LOD_1:8, LCZ_PRIMARY_LOD_1:104, LCZ_PRIMARY_N_LOD_1:104, LCZ_PRIMARY_NE_LOD_1:null, 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] + + assertTrue(values == expectedValues) + + values = h2GIS.firstRow("SELECT * EXCEPT(THE_GEOM) FROM $grid_scale WHERE id_row = 5 AND id_col = 5 ".toString()) + + expectedValues=[ID_GRID:40, ID_ROW:5, ID_COL:5, 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, ID_ROW_LOD_1:2, ID_COL_LOD_1:1, LCZ_PRIMARY_URBAN_LOD_1:4, LCZ_PRIMARY_COOL_LOD_1:5, LCZ_PRIMARY_LOD_1:104, LCZ_PRIMARY_N_LOD_1:104, LCZ_PRIMARY_NE_LOD_1:104, 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] + + assertTrue(values == expectedValues) + + h2GIS.dropTable("grid") + } +} 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..eeb00513ec 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy @@ -202,7 +202,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 +224,13 @@ class SpatialUnitsTests { def countRows = postGIS.firstRow "select count(*) as numberOfRows from $outputTable" assert 100 == countRows.numberOfRows } + + @Test + void sprawlAreasTest() { + String grid_indicators = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Dijon/grid_indicators.geojson", true) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_indicators) + + h2GIS.save(sprawl_areas, "/tmp/sprawl_areas.geojson", true) + + } } \ No newline at end of file diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy index a858625bc9..ceac004de8 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy @@ -52,16 +52,12 @@ class SprawlIndicatorsTests { String grid_indicators = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Dijon/grid_indicators.geojson", true) - String scale_grid = Geoindicators.LCZIndicators.scaling_grid(h2GIS, grid_indicators,1) - // h2GIS.save(scale_grid, "/tmp/scale_grid.csv", true) + def sprawlLayer = Geoindicators.LCZIndicators.compute_sprawl_areas(h2GIS, grid_indicators) - //GRID SCALES - - /*def sprawlLayer = Geoindicators.SprawlIndicators.compute_sprawl_areas(h2GIS, grid_indicators) - - h2GIS.save(sprawlLayer, "/tmp/sprawl_areas.geojson", true) + h2GIS.save(sprawlLayer, "/tmp/sprawl_areas_fractions.geojson", true) + /* String grid_distances = Geoindicators.SprawlIndicators.grid_distances(h2GIS, sprawlLayer, grid_indicators) h2GIS.save(grid_distances, "/tmp/sprawl_areas_distances.geojson", true) From eda745bea3e1322bf60e3716be9801453a0f8677 Mon Sep 17 00:00:00 2001 From: ebocher Date: Fri, 23 Feb 2024 08:21:11 +0100 Subject: [PATCH 05/12] Keep all grid columns --- .../geoclimate/geoindicators/GridIndicators.groovy | 7 +------ .../geoclimate/geoindicators/GridIndicatorsTests.groovy | 1 + 2 files changed, 2 insertions(+), 6 deletions(-) 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 2f40d58eb8..d123487707 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy @@ -134,9 +134,7 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int //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, + 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+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, @@ -267,17 +265,14 @@ String gridDistances(JdbcDataSource datasource, String input_polygons, String gr 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") diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy index fb2e9544d3..7d5d073a3a 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy @@ -5,6 +5,7 @@ import org.junit.jupiter.api.BeforeEach 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 From f960cc5619ed9caa32086026abc46dda1d2f5eda Mon Sep 17 00:00:00 2001 From: ebocher Date: Mon, 4 Mar 2024 17:22:58 +0100 Subject: [PATCH 06/12] Improve multiscale grid --- .../geoindicators/GridIndicators.groovy | 66 ++++++++----------- .../geoindicators/SpatialUnits.groovy | 6 +- .../geoindicators/GridIndicatorsTests.groovy | 35 ++++++++-- .../geoindicators/SpatialUnitsTests.groovy | 17 ++++- 4 files changed, 78 insertions(+), 46 deletions(-) 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 d123487707..7a538febfc 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy @@ -49,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) { @@ -124,11 +124,11 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int def grid_scaling_indices = postfix("grid_scaling_indices") def grid_levels_query = [] int grid_offset = 3 - int offsetCol = 0 + 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) AS ID_COL_LOD_$i" + "(CAST (ABS(ID_COL-1)/${level} AS INT)+$offsetCol-1) AS ID_COL_LOD_$i" offsetCol++ } @@ -136,7 +136,7 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int 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+1=a.ID_COL) AS LCZ_PRIMARY_NE, + (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, @@ -145,9 +145,9 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int (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 tablesToDrop = [] def tableLevelToJoin = grid_scaling_indices - tablesToDrop< Date: Tue, 5 Mar 2024 14:09:31 +0100 Subject: [PATCH 07/12] Add unit tests --- .../geoindicators/SpatialUnits.groovy | 42 ++++--- .../geoindicators/SpatialUnitsTests.groovy | 115 +++++++++++++++++- 2 files changed, 137 insertions(+), 20 deletions(-) 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 fc2f246115..67b6041a8c 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'))"] } } @@ -637,7 +637,6 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, flo error("No grid cells to compute the sprawl areas layer") return } - def gridCols = datasource.getColumnNames(grid_indicators) def lcz_columns_urban = ["LCZ_PRIMARY_1", "LCZ_PRIMARY_2", "LCZ_PRIMARY_3", "LCZ_PRIMARY_4", "LCZ_PRIMARY_5", "LCZ_PRIMARY_6", "LCZ_PRIMARY_7", @@ -646,8 +645,20 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, flo if (lcz_fractions.size() > 0) { def outputTableName = postfix("sprawl_areas") - def tmp_sprawl = postfix("sprawl_tmp") - datasource.execute(""" + 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(st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat + join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat + join=mitre mitre_limit=2')) as the_geom from + st_explode('( + SELECT ST_UNION(st_removeholes(ST_UNION(ST_ACCUM(the_geom)))) + AS THE_GEOM FROM $grid_indicators where ${lcz_fractions.join("+")} >= $fraction + )');""".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 st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat @@ -655,8 +666,9 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, flo st_explode('( SELECT ST_UNION(st_removeholes(ST_UNION(ST_ACCUM(the_geom)))) AS THE_GEOM FROM $grid_indicators where ${lcz_fractions.join("+")} >= $fraction - )'); - CREATE TABLE $outputTableName as SELECT CAST((row_number() over()) as Integer) as id, st_removeholes(the_geom) as the_geom + )');""".toString()) + + datasource.execute("""CREATE TABLE $outputTableName as SELECT CAST((row_number() over()) as Integer) as id, st_removeholes(the_geom) as the_geom FROM ST_EXPLODE('( SELECT @@ -666,8 +678,8 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, flo FROM ST_EXPLODE(''$tmp_sprawl''))') where st_isempty(st_buffer(the_geom, -$distance))=false; DROP TABLE IF EXISTS $tmp_sprawl; """.toString()) - - return outputTableName + return outputTableName + } } error("No LCZ fractions columns to compute the sprawl areas layer") return @@ -679,7 +691,7 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, flo * @param input_polygons a layer that contains polygons * @author Erwan Bocher (CNRS) */ -String inverse_geometries(JdbcDataSource datasource, String input_polygons) { +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; 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 26680673dc..20982dac1a 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy @@ -24,6 +24,7 @@ import org.junit.jupiter.api.BeforeEach 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 +32,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 @@ -226,13 +229,13 @@ class SpatialUnitsTests { } @Test - void sprawlAreasTest() { + 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, 0 AS LCZ_PRIMARY_1, 1 AS LCZ_PRIMARY_104 FROM - ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 100, 100); + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); --Center cell urban UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 5; UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 4; @@ -240,10 +243,112 @@ class SpatialUnitsTests { UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 6; UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 6; """.toString()) - //String grid_indicators = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Dijon/grid_indicators.geojson", true) - String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 10) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 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, 0 AS LCZ_PRIMARY_1, 1 AS LCZ_PRIMARY_104 FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + """.toString()) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 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, 0 AS LCZ_PRIMARY_1, 1 AS LCZ_PRIMARY_104 FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 6; + """.toString()) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 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, 0 AS LCZ_PRIMARY_1, 1 AS LCZ_PRIMARY_104 FROM + ST_MakeGrid('POLYGON((0 0, 9 0, 9 9, 0 0))'::GEOMETRY, 1, 1); + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 6; + + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 9 AND id_col = 9; + UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 1 AND id_col = 1; + """.toString()) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 0) + assertEquals(3, h2GIS.firstRow("select count(*) as count from $sprawl_areas".toString()).count) + assertEquals(11, 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()) + } - h2GIS.save(sprawl_areas, "/tmp/sprawl_areas.geojson", true) + @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)) } } \ No newline at end of file From 8d6e4ca2d463ed4b90efcd42a5cfc175bec511cf Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 5 Mar 2024 15:59:39 +0100 Subject: [PATCH 08/12] Grid distances unit tests --- .../geoindicators/GridIndicators.groovy | 6 ++-- .../geoindicators/GridIndicatorsTests.groovy | 29 ++++++++++++++++++- 2 files changed, 31 insertions(+), 4 deletions(-) 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 7a538febfc..b60caa21be 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy @@ -266,10 +266,10 @@ String gridDistances(JdbcDataSource datasource, String input_polygons, String gr return } int epsg = datasource.getSrid(grid) - def outputTableName = postfix("grid_distance") + def outputTableName = postfix("grid_distances") datasource.execute(""" DROP TABLE IF EXISTS $outputTableName; - CREATE TABLE $outputTableName ( THE_GEOM GEOMETRY,ID INT, DISTANCE FLOAT); + CREATE TABLE $outputTableName (THE_GEOM GEOMETRY,ID INT, DISTANCE FLOAT); """.toString()) datasource.createSpatialIndex(input_polygons) @@ -285,7 +285,7 @@ String gridDistances(JdbcDataSource datasource, String input_polygons, String gr 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() + stmt.addBatch "insert into $outputTableName values(ST_GEOMFROMTEXT('${cell_geom}',$epsg), ${cell.id},${distance})".toString() } } } diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy index 1873ebf0b6..650788542f 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy @@ -79,7 +79,6 @@ class GridIndicatorsTests { 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, nb_levels) - for (int i in 1..nb_levels) { def grid_lod = "grid_lod_$i" h2GIS.execute(""" @@ -93,6 +92,34 @@ class GridIndicatorsTests { 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) } } From 9ad3f1f64214d97ec0c9eaac80163a4c0e76b5de Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 5 Mar 2024 17:16:51 +0100 Subject: [PATCH 09/12] Grid distances unit tests --- .../geoclimate/geoindicators/GridIndicatorsTests.groovy | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy index 650788542f..9085c53897 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy @@ -59,13 +59,13 @@ class GridIndicatorsTests { 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, ID_ROW_LOD_1: 1, ID_COL_LOD_1: 1, 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_WARN_NE_LOD_1: 4, LCZ_WARN_E_LOD_1: null, LCZ_WARN_SE_LOD_1: null, LCZ_WARN_S_LOD_1: null, LCZ_WARN_SW_LOD_1: null, LCZ_WARN_W_LOD_1: null, LCZ_WARN_NW_LOD_1: null] + 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, 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_WARN_NE_LOD_1:4, LCZ_WARN_E_LOD_1:null, LCZ_WARN_SE_LOD_1:null, LCZ_WARN_S_LOD_1:null, LCZ_WARN_SW_LOD_1:null, LCZ_WARN_W_LOD_1:null, LCZ_WARN_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_WARN_NE_LOD_2:null, LCZ_WARN_E_LOD_2:null, LCZ_WARN_SE_LOD_2:null, LCZ_WARN_S_LOD_2:null, LCZ_WARN_SW_LOD_2:null, LCZ_WARN_W_LOD_2:null, LCZ_WARN_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, ID_ROW_LOD_1: 2, ID_COL_LOD_1: 2, 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_WARN_NE_LOD_1: 5, LCZ_WARN_E_LOD_1: null, LCZ_WARN_SE_LOD_1: null, LCZ_WARN_S_LOD_1: null, LCZ_WARN_SW_LOD_1: 1, LCZ_WARN_W_LOD_1: null, LCZ_WARN_NW_LOD_1: null] + 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, 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_WARN_NE_LOD_1:5, LCZ_WARN_E_LOD_1:null, LCZ_WARN_SE_LOD_1:null, LCZ_WARN_S_LOD_1:null, LCZ_WARN_SW_LOD_1:1, LCZ_WARN_W_LOD_1:null, LCZ_WARN_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_WARN_NE_LOD_2:null, LCZ_WARN_E_LOD_2:null, LCZ_WARN_SE_LOD_2:null, LCZ_WARN_S_LOD_2:null, LCZ_WARN_SW_LOD_2:null, LCZ_WARN_W_LOD_2:null, LCZ_WARN_NW_LOD_2:null] assertTrue(values == expectedValues) From 40b14781cee5fee879cf65b9040947ad8d4904d1 Mon Sep 17 00:00:00 2001 From: ebocher Date: Fri, 8 Mar 2024 14:28:10 +0100 Subject: [PATCH 10/12] Improve multiscale grid indicators + sprawl areas --- .../geoindicators/GenericIndicators.groovy | 2 +- .../geoindicators/GridIndicators.groovy | 38 ++++++++- .../geoindicators/SpatialUnits.groovy | 59 ++++++------- .../geoindicators/GridIndicatorsTests.groovy | 4 +- .../geoindicators/SpatialUnitsTests.groovy | 82 +++++++++++-------- 5 files changed, 110 insertions(+), 75 deletions(-) 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 21b9025751..799402b5d4 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 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); 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 b60caa21be..408060db33 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy @@ -142,12 +142,28 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int (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 + (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 = grid_scaling_indices + 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) { @@ -229,14 +245,28 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int 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 + 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.dropTable(tablesToDrop) - + /*def output = postfix("lcz_tiles") + datasource.execute(""" + DROP TABLE IF EXISTS $output; + CREATE TABLE $output as + select *, + SUM(CASE WHEN LCZ_PRIMARY_N in (1,2,3,4,5,6,7,8,9,10,105) + or LCZ_PRIMARY_NE in (1,2,3,4,5,6,7,8,9,10,105) + or LCZ_PRIMARY_E in (1,2,3,4,5,6,7,8,9,10,105) + or LCZ_PRIMARY_SE in (1,2,3,4,5,6,7,8,9,10,105) + or LCZ_PRIMARY_S in (1,2,3,4,5,6,7,8,9,10,105) + or LCZ_PRIMARY_SW in (1,2,3,4,5,6,7,8,9,10,105) + or LCZ_PRIMARY_W in (1,2,3,4,5,6,7,8,9,10,105) + or LCZ_PRIMARY_NW in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 else 0 end) AS LCZ_WARM + FROM $tableLevelToJoin group by the_geom """.toString())*/ return tableLevelToJoin } 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 67b6041a8c..69a0aab597 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -613,24 +613,20 @@ String createGrid(JdbcDataSource datasource, Geometry geometry, double deltaX, d * 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 fraction value to select the cells of the grid to be included in the sprawl areas. * @param distance value to erode (delete) small sprawl areas * @author Erwan Bocher (CNRS) */ -String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, float fraction = 0.65, float distance = 100) { +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 (fraction <= 0) { - error("Please set a fraction greater than 0") - return - } if (distance < 0) { - error("Please set a fraction greater or equal than 0") + error("Please set a distance greater or equal than 0") return } if (datasource.getRowCount(grid_indicators) == 0) { @@ -638,50 +634,43 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, flo return } def gridCols = datasource.getColumnNames(grid_indicators) - - def lcz_columns_urban = ["LCZ_PRIMARY_1", "LCZ_PRIMARY_2", "LCZ_PRIMARY_3", "LCZ_PRIMARY_4", "LCZ_PRIMARY_5", "LCZ_PRIMARY_6", "LCZ_PRIMARY_7", - "LCZ_PRIMARY_8", "LCZ_PRIMARY_9", "LCZ_PRIMARY_10", "LCZ_PRIMARY_105"] - def lcz_fractions = gridCols.intersect(lcz_columns_urban) - - if (lcz_fractions.size() > 0) { + 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(st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat - join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat - join=mitre mitre_limit=2')) as the_geom from - st_explode('( - SELECT ST_UNION(st_removeholes(ST_UNION(ST_ACCUM(the_geom)))) - AS THE_GEOM FROM $grid_indicators where ${lcz_fractions.join("+")} >= $fraction - )');""".toString()) + 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,104,103,107, 106))')""".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 st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat - join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat - join=mitre mitre_limit=2') as the_geom from - st_explode('( - SELECT ST_UNION(st_removeholes(ST_UNION(ST_ACCUM(the_geom)))) - AS THE_GEOM FROM $grid_indicators where ${lcz_fractions.join("+")} >= $fraction - )');""".toString()) - - datasource.execute("""CREATE TABLE $outputTableName as SELECT CAST((row_number() over()) as Integer) as id, st_removeholes(the_geom) as the_geom + 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,104,103,107, 106))') + 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_buffer(st_union(st_accum(st_buffer(the_geom,$distance, ''quad_segs=2 endcap=flat + 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''))') where st_isempty(st_buffer(the_geom, -$distance))=false; + -$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 fractions columns to compute the sprawl areas layer") + error("No LCZ_PRIMARY column to compute the sprawl areas layer") return } diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy index 9085c53897..cc026cf3d7 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy @@ -59,13 +59,13 @@ class GridIndicatorsTests { 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, 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_WARN_NE_LOD_1:4, LCZ_WARN_E_LOD_1:null, LCZ_WARN_SE_LOD_1:null, LCZ_WARN_S_LOD_1:null, LCZ_WARN_SW_LOD_1:null, LCZ_WARN_W_LOD_1:null, LCZ_WARN_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_WARN_NE_LOD_2:null, LCZ_WARN_E_LOD_2:null, LCZ_WARN_SE_LOD_2:null, LCZ_WARN_S_LOD_2:null, LCZ_WARN_SW_LOD_2:null, LCZ_WARN_W_LOD_2:null, LCZ_WARN_NW_LOD_2:null] + 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_WARN_NE_LOD_1:4, LCZ_WARN_E_LOD_1:null, LCZ_WARN_SE_LOD_1:null, LCZ_WARN_S_LOD_1:null, LCZ_WARN_SW_LOD_1:null, LCZ_WARN_W_LOD_1:null, LCZ_WARN_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_WARN_NE_LOD_2:null, LCZ_WARN_E_LOD_2:null, LCZ_WARN_SE_LOD_2:null, LCZ_WARN_S_LOD_2:null, LCZ_WARN_SW_LOD_2:null, LCZ_WARN_W_LOD_2:null, LCZ_WARN_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, 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_WARN_NE_LOD_1:5, LCZ_WARN_E_LOD_1:null, LCZ_WARN_SE_LOD_1:null, LCZ_WARN_S_LOD_1:null, LCZ_WARN_SW_LOD_1:1, LCZ_WARN_W_LOD_1:null, LCZ_WARN_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_WARN_NE_LOD_2:null, LCZ_WARN_E_LOD_2:null, LCZ_WARN_SE_LOD_2:null, LCZ_WARN_S_LOD_2:null, LCZ_WARN_SW_LOD_2:null, LCZ_WARN_W_LOD_2:null, LCZ_WARN_NW_LOD_2:null] + 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_WARN_NE_LOD_1:5, LCZ_WARN_E_LOD_1:null, LCZ_WARN_SE_LOD_1:null, LCZ_WARN_S_LOD_1:null, LCZ_WARN_SW_LOD_1:1, LCZ_WARN_W_LOD_1:null, LCZ_WARN_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_WARN_NE_LOD_2:null, LCZ_WARN_E_LOD_2:null, LCZ_WARN_SE_LOD_2:null, LCZ_WARN_S_LOD_2:null, LCZ_WARN_SW_LOD_2:null, LCZ_WARN_W_LOD_2:null, LCZ_WARN_NW_LOD_2:null] assertTrue(values == expectedValues) 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 20982dac1a..511f76c07b 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy @@ -234,16 +234,17 @@ class SpatialUnitsTests { h2GIS.execute(""" --Grid values DROP TABLE IF EXISTS grid; - CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid, 0 AS LCZ_PRIMARY_1, 1 AS LCZ_PRIMARY_104 FROM + 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= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 5; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 4; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 5; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 6; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 6; + 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 sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 0) + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"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) } @@ -254,10 +255,11 @@ class SpatialUnitsTests { h2GIS.execute(""" --Grid values DROP TABLE IF EXISTS grid; - CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid, 0 AS LCZ_PRIMARY_1, 1 AS LCZ_PRIMARY_104 FROM + 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 sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 0) + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"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()) } @@ -268,18 +270,19 @@ class SpatialUnitsTests { h2GIS.execute(""" --Grid values DROP TABLE IF EXISTS grid; - CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid, 0 AS LCZ_PRIMARY_1, 1 AS LCZ_PRIMARY_104 FROM + 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= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 4; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 5; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 6; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 4; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 6; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 4; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 5; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 6; + 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 sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 0) + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"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) } @@ -290,23 +293,26 @@ class SpatialUnitsTests { h2GIS.execute(""" --Grid values DROP TABLE IF EXISTS grid; - CREATE TABLE grid AS SELECT * EXCEPT(ID), id as id_grid, 0 AS LCZ_PRIMARY_1, 1 AS LCZ_PRIMARY_104 FROM + 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= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 4; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 5; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 4 AND id_col = 6; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 4; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 5 AND id_col = 6; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 4; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 5; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 6 AND id_col = 6; + 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= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 9 AND id_col = 9; - UPDATE grid SET LCZ_PRIMARY_1= 1 , LCZ_PRIMARY_104 =0 WHERE id_row = 1 AND id_col = 1; + 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 sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 0.65, 0) - assertEquals(3, h2GIS.firstRow("select count(*) as count from $sprawl_areas".toString()).count) - assertEquals(11, h2GIS.firstRow("select st_area(st_accum(the_geom)) as area from $sprawl_areas".toString()).area, 0.0001) + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"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 @@ -351,4 +357,14 @@ class SpatialUnitsTests { 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 + void sprawlAreasTestIntegration() { + //Data for test + String path = "/home/ebocher/Autres/data/geoclimate/uhi_lcz/Dijon/" + String data = h2GIS.load("${path}grid_indicators.geojson") + String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,data,1) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 0) + h2GIS.save(sprawl_areas, "/tmp/sprawl_areas_indic.geojson", true) + } } \ No newline at end of file From f4fbeaa5b7e9ead4c335a78628af2d789985710f Mon Sep 17 00:00:00 2001 From: ebocher Date: Thu, 14 Mar 2024 16:40:02 +0100 Subject: [PATCH 11/12] Improve multiscale grid indicators + sprawl areas Clean code --- .../orbisgis/geoclimate/Geoindicators.groovy | 1 - .../geoindicators/GridIndicators.groovy | 124 ++++++++-------- .../geoindicators/LCZIndicators.groovy | 136 ------------------ .../geoindicators/SpatialUnits.groovy | 36 ++++- .../geoindicators/GridIndicatorsTests.groovy | 8 +- .../geoindicators/SpatialUnitsTests.groovy | 29 ++-- .../SprawlIndicatorsTests.groovy | 93 ------------ .../geoclimate/osm/WorflowOSMTest.groovy | 29 ++-- 8 files changed, 140 insertions(+), 316 deletions(-) delete mode 100644 geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy delete mode 100644 geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/Geoindicators.groovy index 567cc960c9..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 LCZIndicators = new LCZIndicators() public static WorkflowUtilities = new WorkflowUtilities() //Cache the XStream models 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 408060db33..0b66997b14 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy @@ -97,6 +97,7 @@ String gridPopulation(JdbcDataSource datasource, String gridTable, String popula 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. @@ -105,12 +106,13 @@ String gridPopulation(JdbcDataSource datasource, String gridTable, String popula * * @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, int nb_levels = 1) { +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 @@ -119,6 +121,14 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int 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") @@ -176,65 +186,65 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int 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()) + 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()) + 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_WARN_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_WARN_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_WARN_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_WARN_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_WARN_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_WARN_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_WARN_NW_LOD_${i}, - - FROM $lcz_count_lod_mode as a; """.toString()) + 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 @@ -248,25 +258,11 @@ String multiscaleLCZGrid(JdbcDataSource datasource, String grid_indicators, int 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; + group by a.ID_ROW_LOD_${i}, a.ID_COL_LOD_${i} , a.id_grid; """.toString()) tableLevelToJoin = grid_level_join } datasource.dropTable(tablesToDrop) - /*def output = postfix("lcz_tiles") - datasource.execute(""" - DROP TABLE IF EXISTS $output; - CREATE TABLE $output as - select *, - SUM(CASE WHEN LCZ_PRIMARY_N in (1,2,3,4,5,6,7,8,9,10,105) - or LCZ_PRIMARY_NE in (1,2,3,4,5,6,7,8,9,10,105) - or LCZ_PRIMARY_E in (1,2,3,4,5,6,7,8,9,10,105) - or LCZ_PRIMARY_SE in (1,2,3,4,5,6,7,8,9,10,105) - or LCZ_PRIMARY_S in (1,2,3,4,5,6,7,8,9,10,105) - or LCZ_PRIMARY_SW in (1,2,3,4,5,6,7,8,9,10,105) - or LCZ_PRIMARY_W in (1,2,3,4,5,6,7,8,9,10,105) - or LCZ_PRIMARY_NW in (1,2,3,4,5,6,7,8,9,10,105) THEN 1 else 0 end) AS LCZ_WARM - FROM $tableLevelToJoin group by the_geom """.toString())*/ return tableLevelToJoin } diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy deleted file mode 100644 index 1d383fbc4b..0000000000 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/LCZIndicators.groovy +++ /dev/null @@ -1,136 +0,0 @@ -/** - * GeoClimate is a geospatial processing toolbox for environmental and climate studies - * https://github.com/orbisgis/geoclimate. - * - * This code is part of the GeoClimate project. GeoClimate is free software; - * you can redistribute it and/or modify it under the terms of the GNU - * Lesser General Public License as published by the Free Software Foundation; - * version 3.0 of the License. - * - * GeoClimate is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License - * for more details . - * - * - * For more information, please consult: - * https://github.com/orbisgis/geoclimate - * - */ -package org.orbisgis.geoclimate.geoindicators - -import groovy.transform.BaseScript -import org.orbisgis.data.jdbc.JdbcDataSource -import org.orbisgis.geoclimate.Geoindicators - - -@BaseScript Geoindicators geoindicators - - - -/** - * This methods allows to extract the cool area geometries. - * A cool area is continous geometry defined by vegetation and water fractions. - * - * @author Erwan Bocher (CNRS) - * TODO: SPATIAL UNITS - */ -String sprawl_cool_areas(JdbcDataSource datasource, String sprawl, String grid_indicators, float fraction = 0.4, - float distance = 100) { - if (!grid_indicators) { - error("No grid_indicators table to compute the sprawl cool areas layer") - return - } - def gridCols = datasource.getColumnNames(grid_indicators) - - def lcz_columns_urban = ["LCZ_PRIMARY_101", "LCZ_PRIMARY_102", - "LCZ_PRIMARY_103", "LCZ_PRIMARY_104", "LCZ_PRIMARY_107"] - def lcz_fractions = gridCols.intersect(lcz_columns_urban) - - if (lcz_fractions.size() > 0) { - datasource.createSpatialIndex(sprawl) - datasource.createSpatialIndex(grid_indicators) - 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, st_removeholes(the_geom) as the_geom FROM ST_EXPLODE('( - SELECT ST_UNION(ST_ACCUM(a.THE_GEOM)) AS THE_GEOM FROM $grid_indicators as a, $sprawl as b where - a.the_geom && b.the_geom and st_intersects(b.the_geom, ST_POINTONSURFACE(a.the_geom)) and - a.${lcz_fractions.join("+ a.")} >= $fraction - )') where st_isempty(st_buffer(st_convexhull(the_geom), -$distance))=false; - """.toString()) - return outputTableName - } - error("No LCZ_PRIMARY columns to compute the cool areas") - return -} - - -/** - * @author Erwan Bocher (CNRS) - * - * TODO : spatial units - */ -String inverse_cool_areas(JdbcDataSource datasource, String inverse_sprawl, String cool_areas) { - - def outputTableName = postfix("inverse_cool_areas") - def tmp_inverse = postfix("tmp_inverse") - datasource.execute(""" - DROP TABLE IF EXISTS $tmp_inverse; - CREATE TABLE $tmp_inverse as - SELECT st_union(st_accum(the_geom)) as the_geom - from (select the_geom from $inverse_sprawl union all select the_geom from $cool_areas) """.toString()) - - 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 $tmp_inverse; - CREATE TABLE $outputTableName as - SELECT CAST((row_number() over()) as Integer) as id, st_buffer(st_buffer(the_geom,-0.1, 'quad_segs=2 endcap=flat - join=mitre mitre_limit=2'),0.1, 'quad_segs=2 endcap=flat - join=mitre mitre_limit=2') as the_geom - FROM - ST_EXPLODE('( - select st_difference(a.the_geom, st_accum(b.the_geom)) as the_geom from $tmp_extent as a, $tmp_inverse - as b)'); - DROP TABLE IF EXISTS $tmp_extent; - """.toString()) - - return outputTableName -} - - -/** - * @author Erwan Bocher (CNRS) - * - * TODO: spatial units - */ -String buffering_cool_areas(JdbcDataSource datasource, String rsu_lcz, String grid_indicators) { - - float distance = 10 - def outputTable = postfix("buffering_cool_areas") - def merging_lcz = postfix("merging_lcz") - datasource.execute(""" -DROP TABLE IF EXISTS $merging_lcz; -CREATE TABLE $merging_lcz as -select CAST((row_number() over()) as Integer) as id, the_geom from ST_EXPLODE('( -select st_union(st_accum(st_buffer(THE_GEOM, 0.1))) as the_geom FROM $rsu_lcz where LCZ_PRIMARY IN (101, 102, 103,104, 107))') -""".toString()) - - distance = 100 - - datasource.execute(""" - DROP TABLE IF EXISTS $outputTable; - CREATE TABLE $outputTable as -select st_buffer(the_geom, -$distance,'quad_segs=2 endcap=flat - join=mitre mitre_limit=2') as the_geom from ST_EXPLODE('( - select st_union(st_accum(st_buffer(the_geom, $distance, ''quad_segs=2 endcap=flat join=mitre mitre_limit=2''))) as the_geom from $merging_lcz)') as foo""".toString()) - - def sprawl_layer = postfix("sprawl_layer") - datasource.execute(""" - DROP TABLE IF EXISTS $sprawl_layer; - CREATE TABLE $sprawl_layer as - - """.toString()) - - return outputTable -} \ No newline at end of file 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 69a0aab597..b5e29b27e2 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -644,7 +644,7 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, flo 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,104,103,107, 106))')""".toString()) + and LCZ_PRIMARY NOT IN (101, 102,103,104,106, 107))')""".toString()) return outputTableName } else { def tmp_sprawl = postfix("sprawl_tmp") @@ -654,7 +654,7 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, flo 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,104,103,107, 106))') + 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 @@ -694,4 +694,36 @@ String inversePolygonsLayer(JdbcDataSource datasource, String input_polygons) { 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 index cc026cf3d7..a3ec4cc5c4 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy @@ -55,17 +55,17 @@ class GridIndicatorsTests { 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", 2) + 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_WARN_NE_LOD_1:4, LCZ_WARN_E_LOD_1:null, LCZ_WARN_SE_LOD_1:null, LCZ_WARN_S_LOD_1:null, LCZ_WARN_SW_LOD_1:null, LCZ_WARN_W_LOD_1:null, LCZ_WARN_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_WARN_NE_LOD_2:null, LCZ_WARN_E_LOD_2:null, LCZ_WARN_SE_LOD_2:null, LCZ_WARN_S_LOD_2:null, LCZ_WARN_SW_LOD_2:null, LCZ_WARN_W_LOD_2:null, LCZ_WARN_NW_LOD_2:null] + 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_WARN_NE_LOD_1:5, LCZ_WARN_E_LOD_1:null, LCZ_WARN_SE_LOD_1:null, LCZ_WARN_S_LOD_1:null, LCZ_WARN_SW_LOD_1:1, LCZ_WARN_W_LOD_1:null, LCZ_WARN_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_WARN_NE_LOD_2:null, LCZ_WARN_E_LOD_2:null, LCZ_WARN_SE_LOD_2:null, LCZ_WARN_S_LOD_2:null, LCZ_WARN_SW_LOD_2:null, LCZ_WARN_W_LOD_2:null, LCZ_WARN_NW_LOD_2:null] + 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) @@ -78,7 +78,7 @@ class GridIndicatorsTests { 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, nb_levels) + 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(""" 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 511f76c07b..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,6 +21,7 @@ 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 @@ -243,7 +244,7 @@ class SpatialUnitsTests { 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",1) + 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) @@ -258,7 +259,7 @@ class SpatialUnitsTests { 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",1) + 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()) @@ -281,7 +282,7 @@ class SpatialUnitsTests { 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",1) + 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) @@ -307,7 +308,7 @@ class SpatialUnitsTests { 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",1) + 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) @@ -335,7 +336,7 @@ class SpatialUnitsTests { 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; + ((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) @@ -359,12 +360,24 @@ class SpatialUnitsTests { } @Test + @Disabled void sprawlAreasTestIntegration() { //Data for test - String path = "/home/ebocher/Autres/data/geoclimate/uhi_lcz/Dijon/" + 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,1) - String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 0) + 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/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy deleted file mode 100644 index ceac004de8..0000000000 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SprawlIndicatorsTests.groovy +++ /dev/null @@ -1,93 +0,0 @@ -/** - * GeoClimate is a geospatial processing toolbox for environmental and climate studies - * https://github.com/orbisgis/geoclimate. - * - * This code is part of the GeoClimate project. GeoClimate is free software; - * you can redistribute it and/or modify it under the terms of the GNU - * Lesser General Public License as published by the Free Software Foundation; - * version 3.0 of the License. - * - * GeoClimate is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License - * for more details . - * - * - * For more information, please consult: - * https://github.com/orbisgis/geoclimate - * - */ -package org.orbisgis.geoclimate.geoindicators - -import org.junit.jupiter.api.BeforeAll -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.geoclimate.Geoindicators - -import static org.orbisgis.data.H2GIS.open - -class SprawlIndicatorsTests { - - @TempDir - static File folder - private static H2GIS h2GIS - - @BeforeAll - static void beforeAll() { - folder = new File("/tmp") - h2GIS = open(folder.getAbsolutePath() + File.separator + "sprawlindicators") - } - - @Disabled - @Test - void debug_tests() { - - /*String rsu_lcz = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Angers/rsu_lcz.geojson", true) - def buffering_cool_areas = Geoindicators.SprawlIndicators.buffering_cool_areas(h2GIS, rsu_lcz, "") - - h2GIS.save(buffering_cool_areas, "/tmp/buffering_cool_areas.geojson", true)*/ - - - String grid_indicators = h2GIS.load("/home/ebocher/Autres/data/geoclimate/uhi_lcz/Dijon/grid_indicators.geojson", true) - - - def sprawlLayer = Geoindicators.LCZIndicators.compute_sprawl_areas(h2GIS, grid_indicators) - - h2GIS.save(sprawlLayer, "/tmp/sprawl_areas_fractions.geojson", true) - - /* - String grid_distances = Geoindicators.SprawlIndicators.grid_distances(h2GIS, sprawlLayer, grid_indicators) - - h2GIS.save(grid_distances, "/tmp/sprawl_areas_distances.geojson", true) - - - String sprawlLayerInverse = Geoindicators.SprawlIndicators.inverse_geometries(h2GIS, sprawlLayer) - - h2GIS.save(sprawlLayerInverse, "/tmp/spraw_areas_inverse.geojson", true) - - - String grid_distances_out = Geoindicators.SprawlIndicators.grid_distances(h2GIS, sprawlLayerInverse, grid_indicators) - - h2GIS.save(grid_distances_out, "/tmp/sprawl_areas_distances_out.geojson", true) - - def cool_areasLayer = Geoindicators.SprawlIndicators.sprawl_cool_areas(h2GIS, sprawlLayer, grid_indicators) - - h2GIS.save(cool_areasLayer, "/tmp/cool_areas.geojson", true) - - String cool_areas_inverse = Geoindicators.SprawlIndicators.inverse_cool_areas(h2GIS,sprawlLayerInverse, cool_areasLayer) - - h2GIS.save(cool_areas_inverse, "/tmp/cool_areas_inverse.geojson", true) - - - String cool_areas_distances = Geoindicators.SprawlIndicators.grid_distances(h2GIS, cool_areas_inverse, grid_indicators) - - h2GIS.save(cool_areas_distances, "/tmp/cool_areas_distances.geojson", true) - */ - - - } - - -} 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 38d0130766..683e5aa116 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -649,9 +649,10 @@ 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 + 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", @@ -676,9 +677,9 @@ class WorflowOSMTest extends WorkflowAbstractTest { "indicatorUse": ["LCZ"] //, "UTRF", "TEB"] - ]/*,"grid_indicators": [ - "x_size": 200, - "y_size": 200, + ],"grid_indicators": [ + "x_size": grid_size, + "y_size": grid_size, //"rowCol": true, "indicators": ["BUILDING_FRACTION","BUILDING_HEIGHT", "BUILDING_POP", //"BUILDING_TYPE_FRACTION", @@ -688,14 +689,26 @@ class WorflowOSMTest extends WorkflowAbstractTest { //"BUILDING_HEIGHT_WEIGHTED", //"BUILDING_SURFACE_DENSITY", "SEA_LAND_FRACTION", "ASPECT_RATIO",//"SVF", "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"] - ], "worldpop_indicators": true, + ], /* "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 From 104b0eec877ef3e14e2c8861b4f18a7f61beb239 Mon Sep 17 00:00:00 2001 From: ebocher Date: Thu, 14 Mar 2024 16:44:29 +0100 Subject: [PATCH 12/12] Doc --- .../orbisgis/geoclimate/geoindicators/GridIndicators.groovy | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 0b66997b14..2ac325ad9d 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicators.groovy @@ -101,12 +101,12 @@ String gridPopulation(JdbcDataSource datasource, String gridTable, String popula /** * 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 + * 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 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