Use .UnionCascaded() to compute union.

1 files changed, 8 insertions, 10 deletions

diff --git a/gis-observation-map b/gis-observation-map
index e9cda7f..9c58b8e 100755
--- a/gis-observation-map
+++ b/gis-observation-map
@@ -244,7 +244,7 @@ def idx_intersects(idx_path, source_geometry):
     locFieldIdx = layer.FindFieldIndex('location', True)
     locations = []
 
-    union = None
+    union = ogr.Geometry(ogr.wkbMultiPolygon)
     feature = layer.GetNextFeature()
     while feature is not None:
         geometry = feature.GetGeometryRef()
@@ -252,13 +252,10 @@ def idx_intersects(idx_path, source_geometry):
         if geometry.Intersects(source_geometry):
             loc = feature.GetField(locFieldIdx)
             locations.append(loc)
-        if union is None:
-            union = geometry.Clone()
-        else:
-            union = union.Union(geometry)
+        union.AddGeometry(geometry)
         feature = layer.GetNextFeature()
 
-    contains = union.Contains(source_geometry)
+    contains = union.UnionCascaded().Contains(source_geometry)
     ds = None
     return locations, contains
 
@@ -293,7 +290,7 @@ for pt in args.point:
     bbox = ogr.Geometry(ogr.wkbPolygon)
     bbox.AddGeometry(ring)
 
-    clipGeometry = clipGeometry.Union(bbox)
+    clipGeometry.AddGeometry(bbox)
 
 for geom in args.geometry:
     ds = ogr.Open(geom['path'].as_posix(), update=0)
@@ -331,22 +328,23 @@ for geom in args.geometry:
 
             bbox = ogr.Geometry(ogr.wkbPolygon)
             bbox.AddGeometry(ring)
-            clipGeometry = clipGeometry.Union(bbox)
+            clipGeometry.AddGeometry(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)
+            clipGeometry.AddGeometry(geometry)
         elif geometryType == ogr.wkbMultiPolygon:
             for i in range(geometry.GetGeometryCount()):
                 geom2 = geometry.GetGeometryRef(i)
-                clipGeometry = clipGeometry.Union(geom2)
+                clipGeometry.AddGeometry(geom2)
         else:
             raise Exception('Unsuported geometry (try --margin=0)')
 
         feature = layer.GetNextFeature()
 
+clipGeometry = clipGeometry.UnionCascaded()
 if not clipGeometry.IsValid():
     raise Exception('Invalid geometry')