Merge pull request #2044 from MJohnson459/fix-multipolygon

fix: restore Polygon cutline schema for Shapely 2

Author: MJohnson459

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,16 +153,8 @@ 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]
-
-        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)
 
@@ -158,7 +163,7 @@ def compute_linestrings(direction):
             'driver': 'GPKG',
             'schema': {
                 'properties': {},
-                'geometry': 'MultiPolygon'
+                'geometry': 'Polygon'
             }
         }
 

tests/test_cutline.py

@@ -0,0 +1,58 @@
+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):
+        """ 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):
+        """ 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) # matches large polygon
+
+    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")
+
+        # 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(
+                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()