From 746eb3900573972afe8598b0606cc223b8928b79 Mon Sep 17 00:00:00 2001 From: Guilhem Moulin Date: Mon, 2 Oct 2023 17:50:37 +0200 Subject: Use .UnionCascaded() to compute union. --- gis-observation-map | 18 ++++++++---------- 1 file 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') -- cgit v1.2.3