feat: normalize capture-area geometry TDE-1359 (#1288)

### Motivation

When an existing capture-area geometry is being modified (polygon being
removed/added), as a result the polygons points can be moved around.
When the overall geometry ends up being the same (no visual change), we
prefer the capture-area geometry to not show a change (difference
between existing and new/re-processed one).

### Modifications

- Use
[shapely.normalize()](https://shapely.readthedocs.io/en/2.0.3/reference/shapely.normalize.html)
to keep the geometry consistent.

What does `shapely.normalize()`?

According to [the official
documentation](https://shapely.readthedocs.io/en/latest/geometry.html#canonical-form):

> - the coordinates of exterior rings follow a clockwise orientation and
interior rings have a counter-clockwise orientation
> - the starting point of rings is lower left
> - elements in collections are ordered by geometry type: by descending
dimension and multi-types first (MultiPolygon, Polygon, MultiLineString,
LineString, MultiPoint, Point). Multiple elements from the same type are
ordered from right to left and from top to bottom.

Because, "_the coordinates of exterior rings follow a clockwise
orientation and interior rings have a counter-clockwise orientation_" is
going against [RFC
7946](https://datatracker.ietf.org/doc/html/rfc7946#section-3.1.6) - "_A
linear ring MUST follow the right-hand rule with respect to the area it
bounds, i.e., exterior rings are counterclockwise, and holes are
clockwise._", we need to apply `orient()` after `normalize()` the
geometry.

### Verification

- Unit tests
- Other tests to come  in #1224

scripts/stac/imagery/capture_area.py

@@ -3,8 +3,7 @@
 from typing import Any, Sequence
 
 from linz_logger import get_log
-from shapely import BufferCapStyle, BufferJoinStyle, to_geojson, union_all, wkt
-from shapely.constructive import make_valid
+from shapely import BufferCapStyle, BufferJoinStyle, make_valid, to_geojson, union_all, wkt
 from shapely.geometry.base import BaseGeometry
 from shapely.ops import orient
 
@@ -58,17 +57,22 @@ def merge_polygons(polygons: Sequence[BaseGeometry], buffer_distance: float) ->
     for poly in polygons:
         # Buffer each polygon to round up to the `buffer_distance`
         buffered_poly = poly.buffer(buffer_distance, cap_style=BufferCapStyle.flat, join_style=BufferJoinStyle.mitre)
-        buffered_polygons.append(buffered_poly)
+        buffered_polygons.append(orient(make_valid(buffered_poly).normalize()))
     union_buffered = union_all(buffered_polygons)
     # Negative buffer back in the polygons
     union_unbuffered = union_buffered.buffer(-buffer_distance, cap_style=BufferCapStyle.flat, join_style=BufferJoinStyle.mitre)
     union_simplified = union_unbuffered.simplify(buffer_distance)
     union_rounded = wkt.loads(wkt.dumps(union_simplified, rounding_precision=8))
-    # Apply right-hand rule winding order (exterior rings should be counter-clockwise) to the geometry
+    # Normalize the geometry in order to store consistent geometry.
+    # Shapely.normalize() reorders the exterior ring based on the lexicographically smallest point (lowest (x, y)),
+    # but also make the coordinates of exterior rings following a clockwise orientation and interior rings having a
+    # counter-clockwise orientation
+    # Because of the above, we need to apply right-hand rule winding order
+    # (exterior rings should be counter-clockwise) to the geometry
     # Ref: https://datatracker.ietf.org/doc/html/rfc7946#section-3.1.6
-    oriented_union_simplified = orient(union_rounded, sign=1.0)
+    oriented_union_simplified = orient(make_valid(union_rounded).normalize(), sign=1.0)
 
-    return make_valid(oriented_union_simplified)
+    return oriented_union_simplified
 
 
 def generate_capture_area(polygons: Sequence[BaseGeometry], gsd: Decimal) -> dict[str, Any]:

scripts/stac/imagery/tests/capture_area_test.py

@@ -15,7 +15,7 @@ def test_merge_polygons() -> None:
     polygons = []
     polygons.append(Polygon([(0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0), (0.0, 1.0)]))
     polygons.append(Polygon([(1.0, 1.0), (2.0, 1.0), (2.0, 0.0), (1.0, 0.0), (1.0, 1.0)]))
-    expected_merged_polygon = Polygon([(1.0, 1.0), (0.0, 1.0), (0.0, 0.0), (2.0, 0.0), (2.0, 1.0), (1.0, 1.0)])
+    expected_merged_polygon = Polygon([(0.0, 0.0), (2.0, 0.0), (2.0, 1.0), (0.0, 1.0), (0.0, 0.0)])
     merged_polygons = merge_polygons(polygons, 0)
 
     print(f"Polygon A: {to_feature(polygons[0])}")
@@ -33,7 +33,7 @@ def test_merge_polygons_with_rounding() -> None:
     polygons.append(Polygon([(0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0), (0.0, 1.0)]))
     # The following polygon is off by 0.1 to the "right" from the previous one
     polygons.append(Polygon([(1.1, 1.0), (2.0, 1.0), (2.0, 0.0), (1.0, 0.0), (1.1, 1.0)]))
-    expected_merged_polygon = Polygon([(2.0, 1.0), (0.0, 1.0), (0.0, 0.0), (2.0, 0.0), (2.0, 1.0)])
+    expected_merged_polygon = Polygon([(0.0, 0.0), (2.0, 0.0), (2.0, 1.0), (0.0, 1.0), (0.0, 0.0)])
     # By giving a buffer distance of 0.1, we want to correct this margin of error and have the two polygons being merged
     merged_polygons = merge_polygons(polygons, 0.1)
 
@@ -55,15 +55,15 @@ def test_merge_polygons_with_rounding_margin_too_big() -> None:
     merged_polygons = merge_polygons(polygons, 0.01)
     expected_merged_polygon = Polygon(
         [
-            (1.0, 1.0),
-            (0.0, 1.0),
             (0.0, 0.0),
             (2.0, 0.0),
             (2.0, 1.0),
             (1.1, 1.0),
             (1.01501863, 0.1501863),
             (1.0, 0.15093536),
             (1.0, 1.0),
+            (0.0, 1.0),
+            (0.0, 0.0),
         ]
     )
 
@@ -214,18 +214,39 @@ def test_capture_area_rounding_decimal_places() -> None:
     )
     capture_area_expected = {
         "geometry": {
-            "type": "Polygon",
             "coordinates": [
                 [
+                    [174.67341848, -37.05127777],
                     [174.67342502, -37.05155033],
                     [174.6734799, -37.05128096],
                     [174.67341848, -37.05127777],
-                    [174.67342502, -37.05155033],
                 ]
             ],
+            "type": "Polygon",
         },
         "type": "Feature",
         "properties": {},
     }
     capture_area = generate_capture_area(polygons, Decimal("1"))
     assert capture_area == capture_area_expected
+
+
+def test_ensure_consistent_geometry() -> None:
+    """Test to ensure removing a polygon from a geometry and putting it back gives the same result."""
+    polygons = []
+    poly_a = Polygon([(0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0), (0.0, 1.0)])
+    poly_b = Polygon([(2.0, 1.0), (2.0, 0.0), (1.0, 0.0), (1.0, 1.0), (2.0, 1.0)])
+    polygons.append(poly_a)
+    polygons.append(poly_b)
+    merged_polygons = merge_polygons(polygons, 0)
+
+    print(f"Polygon A: {to_feature(poly_a)}")
+    print(f"Polygon B: {to_feature(poly_b)}")
+    print(f"Merged: {to_feature(merged_polygons)}")
+
+    merged_polygons_without_poly_b = merged_polygons.difference(poly_b)
+    print(f"GeoJSON merged_polygons_without_poly_b: {to_feature(merged_polygons_without_poly_b)}")
+    merged_polygons_with_poly_b_back = merge_polygons([poly_b, merged_polygons_without_poly_b], 0)
+    print(f"GeoJSON result: {to_feature(merged_polygons_with_poly_b_back)}")
+
+    assert merged_polygons_with_poly_b_back.equals_exact(merged_polygons, tolerance=0.0)

scripts/stac/imagery/tests/collection_test.py

@@ -431,15 +431,15 @@ def test_capture_area_added(fake_collection_metadata: CollectionMetadata, fake_l
 
     with subtests.test():
         assert collection.stac["assets"]["capture_area"]["file:checksum"] in (
-            "1220ba57cd77defc7fa72e140f4faa0846e8905ae443de04aef99bf381d4650c17a0",
+            "122024fd55ea5f9812e80748ee78055893ad7bddbe2b5101dec1cc1b949edc295f51",
             # geos 3.11 - geos 3.12 as yet untested
         )
 
     with subtests.test():
         assert "file:size" in collection.stac["assets"]["capture_area"]
 
     with subtests.test():
-        assert collection.stac["assets"]["capture_area"]["file:size"] in (269,)  # geos 3.11 - geos 3.12 as yet untested
+        assert collection.stac["assets"]["capture_area"]["file:size"] in (239,)  # geos 3.11 - geos 3.12 as yet untested
 
 
 def test_should_make_valid_capture_area() -> None: