From 3550018a2ced3c9ef6baedb0e15ecf06ea344a38 Mon Sep 17 00:00:00 2001 From: Michael Johnson Date: Thu, 2 Jul 2026 13:20:18 +0100 Subject: [PATCH 1/3] fix: restore Polygon cutline schema for Shapely 2 Shapely 2 no longer has `__getitem__` for MultiPolygon, so unary_union results were written whole instead of picking the largest polygon. Commit 75a1e1f1 only fixed one half of that mismatch by changing the schema to MultiPolygon; Fiona 1.9+ still rejects Polygon records against a MultiPolygon schema. Use .geoms to select the largest polygon and declare the GPKG layer as Polygon again. --- opendm/cutline.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/opendm/cutline.py b/opendm/cutline.py index 7f715960f..06d3372b7 100644 --- a/opendm/cutline.py +++ b/opendm/cutline.py @@ -140,9 +140,11 @@ def compute_linestrings(direction): return log.ODM_INFO("Merging polygons") - cutline_polygons = unary_union(polygons) - if not hasattr(cutline_polygons, '__getitem__'): - cutline_polygons = [cutline_polygons] + cutline_union = unary_union(polygons) + if cutline_union.geom_type == 'MultiPolygon': + cutline_polygons = list(cutline_union.geoms) + else: + cutline_polygons = [cutline_union] largest_cutline = cutline_polygons[0] max_area = largest_cutline.area @@ -158,7 +160,7 @@ def compute_linestrings(direction): 'driver': 'GPKG', 'schema': { 'properties': {}, - 'geometry': 'MultiPolygon' + 'geometry': 'Polygon' } } From da6e737ab23633e642aec98a368814aaa1f60757 Mon Sep 17 00:00:00 2001 From: Michael Johnson Date: Thu, 2 Jul 2026 13:49:53 +0100 Subject: [PATCH 2/3] Split out largest_polygon function for testing --- opendm/cutline.py | 27 +++++++++++++---------- tests/test_cutline.py | 51 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+), 12 deletions(-) create mode 100644 tests/test_cutline.py diff --git a/opendm/cutline.py b/opendm/cutline.py index 06d3372b7..553197022 100644 --- a/opendm/cutline.py +++ b/opendm/cutline.py @@ -13,7 +13,7 @@ from skimage.draw import line from skimage.graph import route_through_array import shapely -from shapely.geometry import LineString, mapping, shape +from shapely.geometry import LineString, MultiPolygon, mapping, shape from shapely.ops import polygonize, unary_union def write_raster(data, file): @@ -31,6 +31,19 @@ def write_raster(data, file): with rasterio.open(file, 'w', BIGTIFF="IF_SAFER", **profile) as wout: wout.write(data, 1) +def largest_polygon(polygons): + cutline_union = unary_union(polygons) + if isinstance(cutline_union, MultiPolygon): + parts = list(cutline_union.geoms) + else: + parts = [cutline_union] + + largest = parts[0] + for p in parts[1:]: + if p.area > largest.area: + largest = p + return largest + def compute_cutline(orthophoto_file, crop_area_file, destination, max_concurrency=1, scale=1): if io.file_exists(orthophoto_file) and io.file_exists(crop_area_file): log.ODM_INFO("Computing cutline") @@ -140,18 +153,8 @@ def compute_linestrings(direction): return log.ODM_INFO("Merging polygons") - cutline_union = unary_union(polygons) - if cutline_union.geom_type == 'MultiPolygon': - cutline_polygons = list(cutline_union.geoms) - else: - cutline_polygons = [cutline_union] - - largest_cutline = cutline_polygons[0] + largest_cutline = largest_polygon(polygons) max_area = largest_cutline.area - for p in cutline_polygons: - if p.area > max_area: - max_area = p.area - largest_cutline = p log.ODM_INFO("Largest cutline found: %s m^2" % max_area) diff --git a/tests/test_cutline.py b/tests/test_cutline.py new file mode 100644 index 000000000..c02ef9297 --- /dev/null +++ b/tests/test_cutline.py @@ -0,0 +1,51 @@ +import os +import tempfile +import unittest + +import fiona +from shapely.geometry import Polygon, mapping + +from opendm.cutline import largest_polygon + + +class TestCutline(unittest.TestCase): + def test_largest_polygon_from_single_polygon(self): + p = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) + largest = largest_polygon([p]) + self.assertEqual(largest.geom_type, "Polygon") + self.assertAlmostEqual(largest.area, 1.0) + + def test_largest_polygon_from_multipolygon_union(self): + small = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) + large = Polygon([(0, 0), (3, 0), (3, 2), (0, 2)]) + largest = largest_polygon([small, large]) + self.assertEqual(largest.geom_type, "Polygon") + self.assertAlmostEqual(largest.area, 6.0) + + def test_polygon_writes_to_gpkg_with_polygon_schema(self): + """ Tests that fiona can write our polygon to a GPKG file""" + + geom = largest_polygon([ + Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]), + Polygon([(2, 0), (5, 0), (5, 2), (2, 2)]), + ]) + self.assertEqual(mapping(geom)["type"], "Polygon") + + with tempfile.TemporaryDirectory() as tmp: + path = os.path.join(tmp, "cutline.gpkg") + with fiona.open( + path, + "w", + driver="GPKG", + crs="EPSG:4326", + schema={"geometry": "Polygon", "properties": {}}, + ) as sink: + sink.write({"geometry": mapping(geom), "properties": {}}) + + with fiona.open(path) as src: + feature = next(iter(src)) + self.assertEqual(feature["geometry"]["type"], "Polygon") + + +if __name__ == "__main__": + unittest.main() From 8a6dee7629134826c395a4643fadd11763e329f8 Mon Sep 17 00:00:00 2001 From: Michael Johnson Date: Thu, 2 Jul 2026 14:11:58 +0100 Subject: [PATCH 3/3] Added test comments --- tests/test_cutline.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/tests/test_cutline.py b/tests/test_cutline.py index c02ef9297..90ca215a0 100644 --- a/tests/test_cutline.py +++ b/tests/test_cutline.py @@ -10,17 +10,22 @@ class TestCutline(unittest.TestCase): def test_largest_polygon_from_single_polygon(self): - p = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) + """ Tests that the largest polygon from a single polygon is the polygon itself """ + + p = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) # 1x1 square largest = largest_polygon([p]) + self.assertEqual(largest.geom_type, "Polygon") self.assertAlmostEqual(largest.area, 1.0) def test_largest_polygon_from_multipolygon_union(self): - small = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) - large = Polygon([(0, 0), (3, 0), (3, 2), (0, 2)]) + """ Tests that the largest polygon from a multipolygon union is the largest polygon """ + + small = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) # 1x1 square (area = 1) + large = Polygon([(0, 0), (3, 0), (3, 2), (0, 2)]) # 3x2 rectangle (area = 6) largest = largest_polygon([small, large]) self.assertEqual(largest.geom_type, "Polygon") - self.assertAlmostEqual(largest.area, 6.0) + self.assertAlmostEqual(largest.area, 6.0) # matches large polygon def test_polygon_writes_to_gpkg_with_polygon_schema(self): """ Tests that fiona can write our polygon to a GPKG file""" @@ -31,6 +36,8 @@ def test_polygon_writes_to_gpkg_with_polygon_schema(self): ]) self.assertEqual(mapping(geom)["type"], "Polygon") + # Write out the polygon to quickly verify that Fiona can still process + # our largest polygon. with tempfile.TemporaryDirectory() as tmp: path = os.path.join(tmp, "cutline.gpkg") with fiona.open(