Merge pull request #989 from ebocher/improve_intersection

Improve OSM multipolygons construction
orbisgis · Jul 5, 2024 · 51c1bae · 51c1bae
2 parents a323350 + 1b85331
commit 51c1bae
Show file tree

Hide file tree

Showing 6 changed files with 5,027 additions and 10 deletions.
diff --git a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy
@@ -297,7 +297,7 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils {
             CREATE TABLE $logTableZones (the_geom GEOMETRY(GEOMETRY, 4326), 
             location VARCHAR, info VARCHAR, version  VARCHAR, build_number VARCHAR);""")
         //Find the geometry of the location
-        Geometry geom = dataSource.firstRow("SELECT st_union(st_accum(THE_GEOM)) as the_geom FROM WHERE commune").the_geom
+        Geometry geom = dataSource.firstRow("SELECT st_union(st_accum(THE_GEOM)) as the_geom FROM commune").the_geom
         if (geom == null) {
             dataSource.execute("""INSERT INTO $logTableZones 
                     VALUES(null,'$location', '$message', 

diff --git a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy
@@ -279,7 +279,6 @@ class InputDataFormattingTest {
     @Test
     //enable it to test data extraction from the overpass api
     void extractCreateFormatGISLayers() {
-
         def directory = "/tmp/geoclimate"
 
         File file = new File(directory)
@@ -288,17 +287,17 @@ class InputDataFormattingTest {
         }
         def h2GIS = H2GIS.open("${file.absolutePath + File.separator}osm_gislayers;AUTO_SERVER=TRUE".toString())
 
-        def zoneToExtract = "Vannes"
+        def zoneToExtract = "Nantes"
 
         //def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(zoneToExtract)
         // zoneToExtract = nominatim.bbox
 
         //zoneToExtract = [62.2, 28.2, 62.4, 28.4]
 
-        //zoneToExtract =[45.784554,4.861279,45.796015,4.883981]
+        zoneToExtract =[47.0619, -1.8145005, 47.394558, -1.2849174]
         Map extractData = OSM.InputDataLoading.extractAndCreateGISLayers(h2GIS, zoneToExtract)
 
-        String formatedPlaceName = zoneToExtract.join("-").trim().split("\\s*(,|\\s)\\s*").join("_");
+        String formatedPlaceName = zoneToExtract.join("_").trim().split("\\s*(,|\\s)\\s*").join("_");
 
         if (!formatedPlaceName) {
             formatedPlaceName = zoneToExtract

diff --git a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy
@@ -649,7 +649,7 @@ class WorflowOSMTest extends WorkflowAbstractTest {
         File dirFile = new File(directory)
         dirFile.delete()
         dirFile.mkdir()
-        def location = "Villeurbanne"
+        def location = "Nantes"
         def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(location)
         def grid_size = 250
         location = nominatim.bbox

diff --git a/osmtools/src/main/groovy/org/orbisgis/geoclimate/osmtools/Transform.groovy b/osmtools/src/main/groovy/org/orbisgis/geoclimate/osmtools/Transform.groovy
@@ -529,7 +529,7 @@ def extractRelationsAsPolygons(JdbcDataSource datasource, String osmTablesPrefix
     datasource """
                 DROP TABLE IF EXISTS $relationsPolygonsInnerExploded;
                 CREATE TABLE $relationsPolygonsInnerExploded AS 
-                    SELECT the_geom AS the_geom, id_relation 
+                    SELECT st_makepolygon(the_geom) AS the_geom, id_relation 
                     FROM st_explode('$relationsPolygonsInner') 
                     WHERE ST_STARTPOINT(the_geom) = ST_ENDPOINT(the_geom)
                     AND ST_NPoints(the_geom)>=4; 
@@ -544,10 +544,11 @@ def extractRelationsAsPolygons(JdbcDataSource datasource, String osmTablesPrefix
                 CREATE INDEX ON $relationsPolygonsInnerExploded(id_relation);       
                 DROP TABLE IF EXISTS $relationsMpHoles;
                 CREATE TABLE $relationsMpHoles AS 
-                    SELECT ST_MAKEPOLYGON(ST_EXTERIORRING(a.the_geom), ST_ACCUM(b.the_geom)) AS the_geom, a.ID_RELATION
+                    SELECT ST_MAKEPOLYGON(ST_EXTERIORRING(a.the_geom), ST_ToMultiLine(ST_ACCUM(b.the_geom))) AS the_geom, a.ID_RELATION
                     FROM $relationsPolygonsOuterExploded AS a 
                     LEFT JOIN $relationsPolygonsInnerExploded AS b 
-                    ON( a.ID_RELATION=b.ID_RELATION)
+                    ON( a.ID_RELATION=b.ID_RELATION and  a.the_geom && b.the_geom 
+                    AND st_intersects(a.the_geom, st_pointonsurface(b.the_geom)))
                     GROUP BY a.the_geom, a.id_relation;
                 CREATE INDEX ON $relationsMpHoles(id_relation);
         """.toString()

diff --git a/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/TransformTest.groovy b/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/TransformTest.groovy
@@ -840,5 +840,16 @@ class TransformTest extends AbstractOSMToolsTest {
         h2GIS.save(lines, "/tmp/building.fgb")
     }
 
-
+    @Test
+    void buildGISLayersTest2() {
+        def prefix = "OSM_WATER"
+        assertTrue OSMTools.Loader.load(ds, prefix,
+                new File(this.class.getResource("water.osm").toURI()).getAbsolutePath())
+        //Create building layer
+        def tags = ["natural":"water"]
+        String outputTableName = OSMTools.Transform.toPolygons(ds, prefix, 4326, tags)
+        def res = ds.firstRow("select count(*) as count, sum(ST_NumInteriorRings(the_geom)) as nb_holes from $outputTableName")
+        assertEquals(7, res.count)
+        assertEquals(3, res.nb_holes)
+    }
 }