Convert input geometries to multipolygon via union taking.

1 files changed, 55 insertions, 57 deletions

diff --git a/gis-observation-map b/gis-observation-map
index 8fc48de..a1ce006 100755
--- a/gis-observation-map
+++ b/gis-observation-map
@@ -208,15 +208,8 @@ if projectInstance is not None:
         else:
             projectInstance.addMapLayer(layer, False)
             sourceGroup.addLayer(layer)
-        if i == 0:
-            #layerTreeModel = QgsLayerTreeModel(layerTreeRoot)
-            #layerTreeView = QgsLayerTreeView()
-            #layerTreeView.setModel(layerTreeModel)
-            #layerTreeView.setCurrentLayer(layer)
-            canvas = QgsMapCanvas()
-            canvas.setExtent( layer.extent() )
-
-def gpkg_intersects(path, geometries):
+
+def gpkg_intersects(path, source_geometry):
     ds = gpkg_drv.Open(path, update=0)
     for i in range(ds.GetLayerCount()):
         layer = ds.GetLayerByIndex(i)
@@ -227,14 +220,13 @@ def gpkg_intersects(path, geometries):
         while feature is not None:
             geometry = feature.GetGeometryRef()
             geometry.Transform(transform_srs)
-            for source_geometry in geometries:
-                if geometry.Intersects(source_geometry):
-                    ds = None
-                    return True
+            if geometry.Intersects(source_geometry):
+                ds = None
+                return True
             feature = layer.GetNextFeature()
     return False
 
-def idx_intersects(idx_path, geometries):
+def idx_intersects(idx_path, source_geometry):
     ds = geojson_drv.Open(idx_path, update=0)
     if ds.GetLayerCount() != 1:
         print(f'WARN: {idx_path}: has {ds.GetLayerCount()} != 1 layers', file=sys.stderr)
@@ -250,23 +242,16 @@ def idx_intersects(idx_path, geometries):
     while feature is not None:
         geometry = feature.GetGeometryRef()
         geometry.Transform(transform_srs)
-        for source_geometry in geometries:
-            if geometry.Intersects(source_geometry):
-                loc = feature.GetField(locFieldIdx)
-                locations.append(loc)
-                break
+        if geometry.Intersects(source_geometry):
+            loc = feature.GetField(locFieldIdx)
+            locations.append(loc)
         if union is None:
             union = geometry.Clone()
         else:
             union = union.Union(geometry)
         feature = layer.GetNextFeature()
 
-    contains = True
-    for source_geometry in geometries:
-        if not union.Contains(source_geometry):
-            contains = False
-            break
-
+    contains = union.Contains(source_geometry)
     ds = None
     return locations, contains
 
@@ -279,8 +264,9 @@ def find_qlr(basedir):
         return None
     return qlr[0]
 
-geometries = []
+clipGeometry = ogr.Geometry(ogr.wkbMultiPolygon)
 for pt in args.point:
+    # northing,easting in target SRS
     y, x = pt.split(',', 1)
     x = int(x)
     y = int(y)
@@ -300,7 +286,7 @@ for pt in args.point:
     bbox = ogr.Geometry(ogr.wkbPolygon)
     bbox.AddGeometry(ring)
 
-    geometries.append(bbox)
+    clipGeometry = clipGeometry.Union(bbox)
 
 for geom in args.geometry:
     ds = ogr.Open(geom['path'].as_posix(), update=0)
@@ -314,12 +300,12 @@ for geom in args.geometry:
             geometry = feature.GetGeometryRef()
             if not geometry.IsValid():
                 raise Exception('Invalid geometry')
+            geometryType = geometry.GetGeometryType()
+            geometry = geometry.Clone()
+            geometry.FlattenTo2D()
             geometry.Transform(transform_srs)
 
-            if args.margin is None:
-                geometry = geometry.Clone()
-                geometry.FlattenTo2D()
-            else:
+            if args.margin is not None:
                 minX, maxX, minY, maxY = geometry.GetEnvelope()
                 minX -= args.margin
                 maxX += args.margin
@@ -335,15 +321,27 @@ for geom in args.geometry:
 
                 bbox = ogr.Geometry(ogr.wkbPolygon)
                 bbox.AddGeometry(ring)
-                geometry = bbox
+                clipGeometry = clipGeometry.Union(bbox)
+            elif geometryType == ogr.wkbPoint or geometryType == ogr.wkbMultiPoint:
+                parser.error('Point geometry requires positive --margin')
+            elif geometryType == ogr.wkbLineString or geometryType == ogr.wkbMultiLineString:
+                parser.error('Line geometry requires positive --margin')
+            elif geometryType == ogr.wkbPolygon:
+                clipGeometry = clipGeometry.Union(geometry)
+            elif geometryType == ogr.wkbMultiPolygon:
+                for i in range(geometry.GetGeometryCount()):
+                    geom2 = geometry.GetGeometryRef(i)
+                    clipGeometry = clipGeometry.Union(geom2)
+            else:
+                raise Exception('Unsuported geometry (try --margin=0)')
 
-            geometries.append(geometry)
             feature = layer.GetNextFeature()
 
-if len(args.geometry) == 0 and projectInstance is not None:
-    r = geometries[0].GetGeometryRef(0)
-    pt0X, pt0Y = r.GetPoint_2D(0)
-    pt1X, pt1Y = r.GetPoint_2D(2)
+if not clipGeometry.IsValid():
+    raise Exception('Invalid geometry')
+
+if projectInstance is not None:
+    pt0X, pt1X, pt0Y, pt1Y = clipGeometry.GetEnvelope()
     canvas = QgsMapCanvas()
     extent = QgsRectangle(QgsPointXY(pt0X, pt0Y), QgsPointXY(pt1X, pt1Y))
     canvas.setExtent(extent)
@@ -372,31 +370,31 @@ def geometricFilterPolygon(polygon):
     }
 
 # https://github.com/biodiversitydata-se/SOS/blob/master/Docs/SearchFilter.md#geographics-filter
-def geometricFilter(geometries):
+def geometricFilter(geometry):
     wgs84 = osr.SpatialReference()
     wgs84.ImportFromEPSG(4326) # WGS84
     wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
     towgs84 = osr.CoordinateTransformation(target_srs, wgs84)
 
+    geometry = geometry.Clone()
+    geometry.Transform(towgs84)
+    geometryType = geometry.GetGeometryType()
+
     myGeometries = []
-    for geom in geometries:
-        geom = geom.Clone()
-        geom.Transform(towgs84)
-        geometryType = geom.GetGeometryType()
-        if geometryType == ogr.wkbPoint or geometryType == ogr.wkbMultiPoint:
-            parser.error('Point geometry requires positive --margin')
-        elif geometryType == ogr.wkbLineString or geometryType == ogr.wkbMultiLineString:
-            parser.error('Line geometry requires positive --margin')
-        elif geometryType == ogr.wkbPolygon:
+    if geometryType == ogr.wkbPoint or geometryType == ogr.wkbMultiPoint:
+        parser.error('Point geometry requires positive --margin')
+    elif geometryType == ogr.wkbLineString or geometryType == ogr.wkbMultiLineString:
+        parser.error('Line geometry requires positive --margin')
+    elif geometryType == ogr.wkbPolygon:
+        gfp = geometricFilterPolygon(geometry)
+        myGeometries.append(gfp)
+    elif geometryType == ogr.wkbMultiPolygon:
+        for i in range(geometry.GetGeometryCount()):
+            geom = geometry.GetGeometryRef(i)
             gfp = geometricFilterPolygon(geom)
             myGeometries.append(gfp)
-        elif geometryType == ogr.wkbMultiPolygon:
-            for i in range(geom.GetGeometryCount()):
-                geom2 = geom.GetGeometryRef(i)
-                gfp = geometricFilterPolygon(geom2)
-                myGeometries.append(gfp)
-        else:
-            raise Exception('Unsuported geometry (try --margin=0)')
+    else:
+        raise Exception('Unsuported geometry (try --margin=0)')
 
     return myGeometries
 
@@ -686,7 +684,7 @@ def getTaxonRedlistCategories(taxonLists, i):
 
 if not args.no_observations:
     geographicsFilter = {
-        'geometries': geometricFilter(geometries),
+        'geometries': geometricFilter(clipGeometry),
         'maxAccuracy': 2500,
         'considerObservationAccuracy': True,
         'considerDisturbanceRadius': False
@@ -834,7 +832,7 @@ if args.topo_basedir is not None and projectInstance is not None:
 
         idx = d.joinpath('index.geojson')
         if idx.is_file():
-            tiles, idxContains[topo] = idx_intersects(idx.as_posix(), geometries)
+            tiles, idxContains[topo] = idx_intersects(idx.as_posix(), clipGeometry)
             if len(tiles) == 0:
                 continue
 
@@ -887,7 +885,7 @@ if args.topo_basedir is not None and projectInstance is not None:
                     print(f'WARN: {t}: lacks {gpkg_name}, skipping', file=sys.stderr)
                     continue
 
-                if not gpkg_intersects(gpkg_path.as_posix(), geometries):
+                if not gpkg_intersects(gpkg_path.as_posix(), clipGeometry):
                     continue
 
                 qlr_paths.append(qlr_path)