BoundingBox: add cyclic_antimeridian_split

and improve west-east longitude normalization

related to Open-EO/openeo-geopyspark-driver#1568

Author: soxofaan

openeo_driver/_version.py

@@ -1 +1 @@
-__version__ = "0.139.0a2"
+__version__ = "0.139.0a3"

openeo_driver/util/geometry.py

@@ -427,6 +427,17 @@ class CrsRequired(BoundingBoxException):
     pass
 
 
+def normalize_west_east_longitude(west: float, east: float) -> Tuple[float, float]:
+    """
+    Assuming longitude in degrees (e.g. EPSG:4326):
+    normalize west to range [-180, 180) and east to range (-180, 180]
+    which is useful in bounding box contexts.
+    """
+    west = ((west + 180) % 360) - 180
+    east = 180 - ((180 - east) % 360)
+    return west, east
+
+
 @dataclasses.dataclass(frozen=True)
 class BoundingBox:
     """
@@ -452,6 +463,8 @@ class BoundingBox:
     # TODO: using frozen dataclasss, but with custom __init__ (for CRS normalization) makes things a bit messy.
     #       It might be just easier to implement the "frozen dataclass" behavior manually
 
+    # TODO: do longitude normalization (for EPSG:4326) in constructor instead of on the fly in various places?
+
     west: float
     south: float
     east: float
@@ -548,20 +561,15 @@ def _crs_with_cyclic_x(crs: Union[None, str]) -> bool:
         """
         return crs == "EPSG:4326"
 
-    @staticmethod
-    def _normalize_longitude(x: float) -> float:
-        """Normalize an EPSG:4326 longitude coordinate to the range [-180, 180)"""
-        # TODO: do this normalization in constructor, instead of each time on the fly?
-        return (x + 180) % 360 - 180
-
     def cyclic_antimeridian_crossing(self) -> bool:
         """
         Whether this bounding box uses cyclic longitude coordinates
         and crosses the antimeridian (so that west > east).
         """
-        return self._crs_with_cyclic_x(self.crs) and (
-            self._normalize_longitude(self.west) > self._normalize_longitude(self.east)
-        )
+        if self._crs_with_cyclic_x(self.crs):
+            west, east = normalize_west_east_longitude(west=self.west, east=self.east)
+            return west > east
+        return False
 
     def as_dict(self) -> dict:
         return {
@@ -586,9 +594,9 @@ def as_polygon(self) -> shapely.geometry.Polygon:
         """
         west, east = self.west, self.east
         if self._crs_with_cyclic_x(self.crs):
-            west = self._normalize_longitude(west)
-            east = self._normalize_longitude(east)
+            west, east = normalize_west_east_longitude(west=west, east=east)
             if east < west:
+                # TODO: this assumes "cyclic" implies longitude in degrees (EPSG:4326)
                 east += 360
 
         return shapely.geometry.box(minx=west, miny=self.south, maxx=east, maxy=self.north)
@@ -597,9 +605,9 @@ def as_geometry(self) -> Union[shapely.geometry.Polygon, shapely.geometry.MultiP
         """Get bounding box as a shapely geometry (Polygon or MultiPolygon when crossing antimeridian)"""
         west, east = self.west, self.east
         if self._crs_with_cyclic_x(self.crs):
-            east = self._normalize_longitude(east)
-            west = self._normalize_longitude(west)
+            west, east = normalize_west_east_longitude(west=west, east=east)
             if east < west:
+                # TODO: this assumes "cyclic" implies longitude in degrees (EPSG:4326), and split is at +/-180
                 return shapely.geometry.MultiPolygon(
                     [
                         shapely.geometry.box(west, self.south, 180, self.north),
@@ -625,12 +633,11 @@ def as_geojson(self) -> dict:
     def centroid(self) -> Tuple[float, float]:
         if self._crs_with_cyclic_x(self.crs):
             # Properly handle cyclic longitude coordinates, and antimeridian crossing
-            west = self._normalize_longitude(self.west)
-            east = self._normalize_longitude(self.east)
+            west, east = normalize_west_east_longitude(west=self.west, east=self.east)
             if west <= east:
                 x = 0.5 * (west + east)
             else:
-                x = self._normalize_longitude(0.5 * (west + east + 360))
+                x, _ = normalize_west_east_longitude(west=0.5 * (west + east + 360), east=0)
         else:
             x = 0.5 * (self.west + self.east)
 
@@ -792,3 +799,18 @@ def intersection(self, other: "BoundingBox") -> Union["BoundingBox", None]:
         if west >= east or south >= north:
             return None
         return BoundingBox(west=west, south=south, east=east, north=north, crs=self.crs)
+
+    def cyclic_antimeridian_split(self) -> List["BoundingBox"]:
+        """
+        Split this bounding box into two bounding boxes
+        when it uses cyclic longitude coordinates and crosses the antimeridian.
+        Otherwise, just return this bounding box as is.
+        """
+        if self.cyclic_antimeridian_crossing():
+            # TODO: this assumes "cyclic" implies longitude in degrees (EPSG:4326), and split is at +/-180
+            return [
+                BoundingBox(west=self.west, south=self.south, east=180, north=self.north, crs=self.crs),
+                BoundingBox(west=-180, south=self.south, east=self.east, north=self.north, crs=self.crs),
+            ]
+        else:
+            return [self]

tests/util/test_geometry.py

@@ -24,6 +24,7 @@
     spatial_extent_union,
     validate_geojson_basic,
     epsg_code_or_none,
+    normalize_west_east_longitude,
 )
 
 from ..data import get_path
@@ -865,6 +866,10 @@ def test_centroid(self):
             (170, -160, -175),
             (170 + 360, -160 + 5 * 360, -175),
             (170, -170, -180),
+            (-180, -170, -175),
+            (-190, -180, 175),
+            (170, 180, 175),
+            (180, 190, -175),
         ],
     )
     def test_centroid_epsg4326_wrap_around(self, west, east, expected):
@@ -911,6 +916,18 @@ def test_contains_generic(self):
                 [(-167.1, 51.5), (-167.1 - 360, 51.5), (0, 51), (360, 51.5), (169.9, 52), (169.9 + 360, 52)],
                 [(-189, 51.5), (-171, 51.5), (171, 51.5), (189, 51.5), (0, 50.9), (0, 52.1)],
             ),
+            (
+                # Bounding box with west bound at  -180
+                BoundingBox(-180, 51, -170, 52, crs=4326),
+                [(-180, 51.5), (-175, 51.5), (180, 51.5), (185, 51.5)],
+                [(-180.5, 51.5), (179.5, 51.5)],
+            ),
+            (
+                # Bounding box with east bound at 180
+                BoundingBox(170, 51, 180, 52, crs=4326),
+                [(180, 51.5), (175, 51.5), (-180, 51.5), (-185, 51.5)],
+                [(180.5, 51.5), (-179.5, 51.5)],
+            ),
         ],
     )
     def test_contains_epsg4326_wrap_around(self, bbox, inside, outside):
@@ -1115,6 +1132,60 @@ def test_intersection(self, bbox1, bbox2, expected):
         if bbox1.crs == bbox2.crs:
             assert bbox2.intersection(bbox1) == expected
 
+    @pytest.mark.parametrize(
+        ["bbox", "expected"],
+        [
+            (
+                BoundingBox(1, 2, 3, 4),
+                [BoundingBox(1, 2, 3, 4)],
+            ),
+            (
+                BoundingBox(-170, 2, 170, 4, crs=4326),
+                [BoundingBox(-170, 2, 170, 4, crs=4326)],
+            ),
+            (
+                BoundingBox(175, 2, -170, 4, crs=4326),
+                [
+                    BoundingBox(175, 2, 180, 4, crs=4326),
+                    BoundingBox(-180, 2, -170, 4, crs=4326),
+                ],
+            ),
+            (
+                BoundingBox(175, 2, 180, 4, crs=4326),
+                [BoundingBox(175, 2, 180, 4, crs=4326)],
+            ),
+            (
+                BoundingBox(-180, 2, -170, 4, crs=4326),
+                [BoundingBox(-180, 2, -170, 4, crs=4326)],
+            ),
+        ],
+    )
+    def test_cyclic_antimeridian_split(self, bbox, expected):
+        assert bbox.cyclic_antimeridian_split() == expected
+
+
+def test_normalize_west_east_longitude():
+    assert normalize_west_east_longitude(0, 0) == (0, 0)
+    assert normalize_west_east_longitude(10, 20) == (10, 20)
+    assert normalize_west_east_longitude(10.5, 20.5) == (10.5, 20.5)
+    assert normalize_west_east_longitude(190, 200) == (-170, -160)
+    assert normalize_west_east_longitude(190.5, 200.5) == (-169.5, -159.5)
+
+    assert normalize_west_east_longitude(-10, -20) == (-10, -20)
+    assert normalize_west_east_longitude(-10.5, -20.5) == (-10.5, -20.5)
+    assert normalize_west_east_longitude(-200, -190) == (160, 170)
+    assert normalize_west_east_longitude(-200.5, -190.5) == (159.5, 169.5)
+
+    assert normalize_west_east_longitude(-180, -170) == (-180, -170)
+    assert normalize_west_east_longitude(-190, -180) == (170, 180)
+    assert normalize_west_east_longitude(-180.5, -170.5) == (179.5, -170.5)
+    assert normalize_west_east_longitude(-190.5, -180.5) == (169.5, 179.5)
+
+    assert normalize_west_east_longitude(170, 180) == (170, 180)
+    assert normalize_west_east_longitude(180, 190) == (-180, -170)
+    assert normalize_west_east_longitude(170.5, 180.5) == (170.5, -179.5)
+    assert normalize_west_east_longitude(180.5, 190.5) == (-179.5, -169.5)
+
 
 class TestValidateGeoJSON:
     @staticmethod