Some improvement to extract urban sprawl, add new unit tests and fix …

…issue #981

.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.

.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

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
     }

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

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)
+    }
 }

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)
+
+    }
+
 
 }

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)
 

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()}')""")
     }

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,

pom.xml

@@ -40,7 +40,7 @@
         <!-- Dependencies versions -->
         <junit-version>5.9.1</junit-version>
         <jts-version>1.19.0</jts-version>
-        <h2gis-version>2.2.1</h2gis-version>
+        <h2gis-version>2.2.2-SNAPSHOT</h2gis-version>
         <h2-version>2.2.224</h2-version>
         <cts-version>1.7.1</cts-version>
         <orbisdata-version>2.1.1-SNAPSHOT</orbisdata-version>