From 25d7f7204852183797a81394a9eb645015fba5ad Mon Sep 17 00:00:00 2001 From: ebocher Date: Fri, 21 Jun 2024 09:48:59 +0200 Subject: [PATCH] Some improvement to extract urban sprawl, add new unit tests and fix issue #981 --- .github/workflows/CI release.yml | 4 +- .github/workflows/CI_snapshot.yml | 4 +- .../geoindicators/SpatialUnits.groovy | 38 ++++++++--------- .../WorkflowGeoIndicators.groovy | 4 +- .../geoindicators/SpatialUnitsTests.groovy | 20 +++++++++ .../WorkflowGeoIndicatorsTest.groovy | 42 +++++++++++++++++++ .../geoclimate/osm/InputDataFormatting.groovy | 19 ++++----- .../geoclimate/osm/WorkflowOSM.groovy | 4 +- .../geoclimate/osm/WorflowOSMTest.groovy | 21 ++++++---- pom.xml | 2 +- 10 files changed, 111 insertions(+), 47 deletions(-) diff --git a/.github/workflows/CI release.yml b/.github/workflows/CI release.yml index f8a2f528f2..c69ba73343 100644 --- a/.github/workflows/CI release.yml +++ b/.github/workflows/CI release.yml @@ -70,8 +70,8 @@ jobs: release:prepare release:perform \ -Dusername=$GITHUB_ACTOR -Dpassword=$GITHUB_TOKEN ${VERSION:+"-DdevelopmentVersion="$VERSION"-SNAPSHOT"} env: - MAVEN_USERNAME: ${{ secrets.OSSRH_USERNAME }} - MAVEN_PASSWORD: ${{ secrets.OSSRH_TOKEN }} + MAVEN_USERNAME: ${{ secrets.MVN_CENTRAL_USERNAME }} + MAVEN_PASSWORD: ${{ secrets.MVN_CENTRAL_PASSWORD }} MAVEN_GPG_PASSPHRASE: ${{ secrets.OSSRH_GPG_PASSWORD }} # Export the last git tag into env. diff --git a/.github/workflows/CI_snapshot.yml b/.github/workflows/CI_snapshot.yml index 48066234a0..9173789bbc 100644 --- a/.github/workflows/CI_snapshot.yml +++ b/.github/workflows/CI_snapshot.yml @@ -53,8 +53,8 @@ jobs: - name: Deploy run: mvn deploy --no-transfer-progress --batch-mode env: - MAVEN_USERNAME: ${{ secrets.OSSRH_USERNAME }} - MAVEN_PASSWORD: ${{ secrets.OSSRH_TOKEN }} + MAVEN_USERNAME: ${{ secrets.MVN_CENTRAL_USERNAME }} + MAVEN_PASSWORD: ${{ secrets.MVN_CENTRAL_PASSWORD }} MAVEN_GPG_PASSPHRASE: ${{ secrets.OSSRH_GPG_PASSWORD }} # Build jar with dependencies 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 62fcb6b422..2a34100405 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -661,26 +661,25 @@ String computeSprawlAreas(JdbcDataSource datasource, String grid_indicators, } else { def tmp_sprawl = postfix("sprawl_tmp") datasource.execute(""" - DROP TABLE IF EXISTS $tmp_sprawl, $outputTableName; - create table $tmp_sprawl as - select CAST((row_number() over()) as Integer) as id, st_removeholes(the_geom) as the_geom from ST_EXPLODE('( - select st_union(st_accum(the_geom)) as the_geom from - $grid_indicators where - LCZ_PRIMARY NOT IN (101, 102,103,104,106, 107))') - where st_area(st_buffer(the_geom, -$distance,2)) > 1""".toString()) + DROP TABLE IF EXISTS $tmp_sprawl, $outputTableName; + create table $tmp_sprawl as + select CAST((row_number() over()) as Integer) as id, st_removeholes(the_geom) as the_geom from ST_EXPLODE('( + select st_union(st_accum(the_geom)) as the_geom from + $grid_indicators where + LCZ_PRIMARY NOT IN (101, 102,103,104,106, 107))') + where st_area(st_buffer(the_geom, -$distance,2)) > 1""".toString()) datasource.execute("""CREATE TABLE $outputTableName as SELECT CAST((row_number() over()) as Integer) as id, - the_geom - FROM - ST_EXPLODE('( - SELECT - st_removeholes(st_buffer(st_union(st_accum(st_buffer(st_removeholes(the_geom),$distance, ''quad_segs=2 endcap=flat + the_geom FROM + ST_EXPLODE('( + SELECT + st_removeholes(st_buffer(st_union(st_accum(st_buffer(st_removeholes(the_geom),$distance, ''quad_segs=2 endcap=flat join=mitre mitre_limit=2''))), -$distance, ''quad_segs=2 endcap=flat join=mitre mitre_limit=2'')) as the_geom - FROM ST_EXPLODE(''$tmp_sprawl'') )') where (the_geom is not null or - st_isempty(the_geom) = false) and st_area(st_buffer(the_geom, -$distance,2)) >${distance * distance}; - DROP TABLE IF EXISTS $tmp_sprawl; - """.toString()) + FROM ST_EXPLODE(''$tmp_sprawl'') )') where (the_geom is not null or + st_isempty(the_geom) = false) and st_area(st_buffer(the_geom, -$distance,2)) >${distance * distance}; + DROP TABLE IF EXISTS $tmp_sprawl; + """) return outputTableName } } @@ -741,7 +740,7 @@ String inversePolygonsLayer(JdbcDataSource datasource, String input_polygons, St * * @author Erwan Bocher (CNRS) */ -String extractCoolAreas(JdbcDataSource datasource, String grid_indicators,String polygons_mask, +String extractCoolAreas(JdbcDataSource datasource, String grid_indicators, String polygons_mask, float distance = 50) throws Exception { if (!grid_indicators || !polygons_mask) { throw new IllegalArgumentException("No grid_indicators table to extract the cool areas layer") @@ -757,8 +756,9 @@ String extractCoolAreas(JdbcDataSource datasource, String grid_indicators,String SELECT ST_UNION(ST_ACCUM(a.THE_GEOM)) AS THE_GEOM FROM $grid_indicators as a, $polygons_mask as b where a.LCZ_PRIMARY in (101, 102, 103,104, 106, 107) and - a.the_geom && b.the_geom and st_intersects(st_pointonsurface(a.the_geom), b.the_geom))') ${distance > 0 ? - " where (the_geom is not null or st_isempty(the_geom) = false) and st_area(st_buffer(the_geom, -$distance,2)) >${distance * distance}" : ""}; + a.the_geom && b.the_geom and st_intersects(st_pointonsurface(a.the_geom), b.the_geom))') ${distance > 0 ? + " where (the_geom is not null or st_isempty(the_geom) = false) and st_area(st_buffer(the_geom, -$distance,2)) >${distance * distance}" : "" + }; """.toString()) return outputTableName } diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index ee33d62f91..140cb2a338 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -2187,7 +2187,7 @@ Map sprawlIndicators(JdbcDataSource datasource, String grid_indicators, String i if (list_indicators_upper.intersect(["URBAN_SPRAWL_AREAS", "URBAN_SPRAWL_DISTANCES", "URBAN_SPRAWL_COOL_DISTANCES"]) && grid_indicators) { sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(datasource, grid_indicators, distance ) } - if (sprawl_areas) { + if (sprawl_areas && datasource.getRowCount(sprawl_areas)>0) { //Compute the distances if (list_indicators_upper.contains("URBAN_SPRAWL_DISTANCES")) { String inside_sprawl_areas = Geoindicators.GridIndicators.gridDistances(datasource, sprawl_areas, grid_indicators, id_grid, false) @@ -2209,7 +2209,7 @@ Map sprawlIndicators(JdbcDataSource datasource, String grid_indicators, String i } if (list_indicators_upper.contains("URBAN_SPRAWL_COOL_DISTANCES")) { cool_areas = Geoindicators.SpatialUnits.extractCoolAreas(datasource, grid_indicators, sprawl_areas, (distance / 2) as float) - if (cool_areas) { + if (cool_areas && datasource.getRowCount(cool_areas)>0) { String inverse_cool_areas = Geoindicators.SpatialUnits.inversePolygonsLayer(datasource, sprawl_areas,cool_areas) if (inverse_cool_areas) { tablesToDrop << inverse_cool_areas 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 bf0ba8cee8..6949cd2f36 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy @@ -398,4 +398,24 @@ class SpatialUnitsTests { h2GIS.save(rsu, "/tmp/rsu.fgb", true) } } + + @Test + void sprawlAreasTest5() { + //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, 2500 0, 2500 2500, 0 0))'::GEOMETRY, 250, 250); + --Center cell urban + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 6; + """.toString()) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, "grid", 250/2) + assertEquals(1, h2GIS.firstRow("select count(*) as count from $sprawl_areas".toString()).count) + assertEquals(312500, h2GIS.firstRow("select st_area(the_geom) as area from $sprawl_areas".toString()).area, 0.0001) + } } \ No newline at end of file diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy index 7b42b39edf..c8eea73326 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy @@ -629,5 +629,47 @@ class WorkflowGeoIndicatorsTest { """.toString()) } + @Test + void sprawlIndicators1() { + //Data for test + datasource.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, 2500 0, 2500 2500, 0 0))'::GEOMETRY, 250, 250); + --Center cell urban + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 4; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 5; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 6 AND id_col = 6; + UPDATE grid SET LCZ_PRIMARY= 1 WHERE id_row = 5 AND id_col = 6; + """) + def results = Geoindicators.WorkflowGeoIndicators.sprawlIndicators(datasource, "grid", "id_grid", ["URBAN_SPRAWL_AREAS", "URBAN_SPRAWL_DISTANCES", "URBAN_SPRAWL_COOL_DISTANCES"], + 250/2 ) + def urban_cool_areas = results.urban_cool_areas + def urban_sprawl_areas = results.urban_sprawl_areas + assertEquals(312500, datasource.firstRow("select st_area(the_geom) as area from $urban_sprawl_areas").area, 0.0001) + assertEquals(0, datasource.getRowCount(urban_cool_areas)) + } + + + @Test + void sprawlIndicators2() { + //Data for test only vegetation + datasource.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, 2500 0, 2500 2500, 0 0))'::GEOMETRY, 250, 250); + """.toString()) + def results = Geoindicators.WorkflowGeoIndicators.sprawlIndicators(datasource, "grid", "id_grid", ["URBAN_SPRAWL_AREAS", "URBAN_SPRAWL_DISTANCES", "URBAN_SPRAWL_COOL_DISTANCES"], + 250/2 ) + def urban_cool_areas = results.urban_cool_areas + def urban_sprawl_areas = results.urban_sprawl_areas + assertTrue(datasource.isEmpty(urban_sprawl_areas)) + assertNull(urban_cool_areas) + + } + } diff --git a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy index 92b8b41527..fe8a6d77f3 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy @@ -53,18 +53,18 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone def outputTableName = postfix "INPUT_BUILDING" debug 'Formating building layer' def outputEstimateTableName = "EST_${outputTableName}" - datasource """ + datasource.execute(""" DROP TABLE if exists ${outputEstimateTableName}; CREATE TABLE ${outputEstimateTableName} ( id_build INTEGER, ID_SOURCE VARCHAR) - """.toString() + """) - datasource """ + datasource.execute(""" DROP TABLE if exists ${outputTableName}; CREATE TABLE ${outputTableName} (THE_GEOM GEOMETRY, id_build INTEGER, ID_SOURCE VARCHAR, HEIGHT_WALL FLOAT, HEIGHT_ROOF FLOAT, NB_LEV INTEGER, TYPE VARCHAR, MAIN_USE VARCHAR, ZINDEX INTEGER, ROOF_SHAPE VARCHAR); - """.toString() + """) if (building) { def paramsDefaultFile = this.class.getResourceAsStream("buildingParams.json") def parametersMap = parametersMapping(jsonFilename, paramsDefaultFile) @@ -82,8 +82,10 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone Geometry geomZone = datasource.firstRow("select st_union(st_accum(the_geom)) as the_geom from $zone").the_geom PreparedGeometry pZone = PreparedGeometryFactory.prepare(geomZone) def id_build = 1; - datasource.withBatch(200) { stmt -> + datasource.withBatch(100) { stmt -> datasource.eachRow(queryMapper) { row -> + Geometry geom = row.the_geom + if (pZone.intersects(geom)) { def typeAndUseValues = getTypeAndUse(row, columnNames, mappingTypeAndUse) def use = typeAndUseValues[1] def type = typeAndUseValues[0] @@ -99,8 +101,6 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone def zIndex = getZIndex(row.'layer') String roof_shape = row.'roof:shape' if (formatedHeight.nbLevels > 0 && zIndex >= 0 && type) { - Geometry geom = row.the_geom - if (pZone.intersects(geom)) { def srid = geom.getSRID() for (int i = 0; i < geom.getNumGeometries(); i++) { Geometry subGeom = geom.getGeometryN(i) @@ -136,7 +136,7 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone } } else { def id_build = 1; - datasource.withBatch(200) { stmt -> + datasource.withBatch(100) { stmt -> datasource.eachRow(queryMapper) { row -> def typeAndUseValues = getTypeAndUse(row, columnNames, mappingTypeAndUse) def use = typeAndUseValues[1] @@ -152,7 +152,6 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone def formatedHeight = Geoindicators.WorkflowGeoIndicators.formatHeightsAndNbLevels(heightWall, heightRoof, nbLevels, h_lev_min, type, typeAndLevel) def zIndex = getZIndex(row.'layer') String roof_shape = row.'roof:shape' - if (formatedHeight.nbLevels > 0 && zIndex >= 0 && type) { Geometry geom = row.the_geom def srid = geom.getSRID() @@ -201,7 +200,7 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone datasource.execute("""DROP TABLE IF EXISTS $triangles; CREATE TABLE $triangles as SELECT * FROM st_explode('(SELECT CASE WHEN ST_AREA(THE_GEOM) > 100000 - THEN ST_Tessellate(st_buffer(the_geom, -0.001,2)) ELSE THE_GEOM END AS THE_GEOM, + THEN ST_Tessellate(st_snaptoself(the_geom, -0.001)) ELSE THE_GEOM END AS THE_GEOM, type FROM $urban_areas)')""".toString()) datasource.createSpatialIndex(triangles) diff --git a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy index 0848625790..7d97ddd291 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy @@ -667,12 +667,12 @@ void saveLogZoneTable(JdbcDataSource dataSource, String databaseFolder, String i location VARCHAR, info VARCHAR, version VARCHAR, build_number VARCHAR);""") if (osm_geometry == null) { dataSource.execute("""INSERT INTO $logTableZones - VALUES(null,'$id_zone', '$message', + VALUES(null,'$id_zone', '${message.replace("'","''")replace("")}', '${Geoindicators.version()}', '${Geoindicators.buildNumber()}')""") } else { dataSource.execute("""INSERT INTO $logTableZones - VALUES(st_geomfromtext('${osm_geometry}',4326) ,'$id_zone', '$message', + VALUES(st_geomfromtext('${osm_geometry}',4326) ,'$id_zone', '${message.replace("'","''")}', '${Geoindicators.version()}', '${Geoindicators.buildNumber()}')""") } 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 c8992e317e..f62e40c7a2 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -649,12 +649,15 @@ class WorflowOSMTest extends WorkflowAbstractTest { File dirFile = new File(directory) dirFile.delete() dirFile.mkdir() - def location = "Genève, Suisse" + def location = "Redon" //def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(location) - def grid_size = 100 + def grid_size = 250 //location = nominatim.bbox - location=[46.178404,6.095524,46.222959,6.190109] - //location =[62.2, 28.2, 62.4, 28.4] + //location=[46.178404,6.095524,46.222959,6.190109] + location =[ 48.84017284026897, + 2.3061887733275785, + 48.878115442982086, + 2.36742047202511] def osm_parmeters = [ "description" : "Example of configuration file to run the OSM workflow and store the result in a folder", "geoclimatedb": [ @@ -675,9 +678,9 @@ class WorflowOSMTest extends WorkflowAbstractTest { "parameters" : ["distance" : 0, "rsu_indicators" : [ - "indicatorUse": ["LCZ", "TEB"] //, "UTRF", "TEB"] + "indicatorUse": ["LCZ", "TEB"] //, "UTRF"] - ], "grid_indicators" : [ + ], /*"grid_indicators" : [ "x_size" : grid_size, "y_size" : grid_size, //"rowCol": true, @@ -693,9 +696,9 @@ class WorflowOSMTest extends WorkflowAbstractTest { // "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS", "URBAN_SPRAWL_AREAS", "URBAN_SPRAWL_DISTANCES", - "URBAN_SPRAWL_COOL_DISTANCES"]/*, - "lcz_lod":1*/ - ], "worldpop_indicators": true + "URBAN_SPRAWL_COOL_DISTANCES"]*//*, + "lcz_lod":1*//* + ], "worldpop_indicators": true*/ /* "road_traffic" : true, diff --git a/pom.xml b/pom.xml index 9d2481669b..588e82e3c0 100644 --- a/pom.xml +++ b/pom.xml @@ -40,7 +40,7 @@ 5.9.1 1.19.0 - 2.2.1 + 2.2.2-SNAPSHOT 2.2.224 1.7.1 2.1.1-SNAPSHOT