Add antimeridian helpers to utils

CHANGELOG.md

@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `stac lint` using [stac-check](https://github.com/philvarner/stac-check) ([#254](https://github.com/stac-utils/stactools/pull/254))
 - `stac info` now has option `-s` or `--skip_items` to skip counting and printing catalog item info ([#260](https://github.com/stac-utils/stactools/pull/260))
 - updated `livehtml` target for documentation build to work with latest `sphinx-autobuild` ([#261](https://github.com/stac-utils/stactools/pull/261))
+- Antimeridian helpers ([#259](https://github.com/stac-utils/stactools/pull/259))
 
 ## [0.2.6] - 2022-02-15
 

src/stactools/core/utils/antimeridian.py

@@ -0,0 +1,131 @@
+from enum import Enum, auto
+import math
+from typing import Optional
+
+import shapely.affinity
+import shapely.geometry
+import shapely.ops
+from pystac import Item
+from shapely.geometry import Polygon, MultiPolygon, LineString
+
+
+class Strategy(Enum):
+    """Strategy for handling antimeridian-crossing polygons.
+
+        - SPLIT: split the polygon at the antimeridian
+        - NORMALIZE: keep the polygon as one, but extend it to > 180 or < -180.
+    """
+    SPLIT = auto()
+    NORMALIZE = auto()
+
+
+def fix_item(item: Item, strategy: Strategy) -> None:
+    """Modifies an item in-place to deal with antimeridian issues.
+
+    If the item's geometry is not a `Polygon`, raises a `ValueError`.
+
+    Args:
+        item (pystac.Item): The item to be modified.
+    """
+    geometry = shapely.geometry.shape(item.geometry)
+    if not isinstance(geometry, Polygon):
+        raise ValueError(
+            f"Can only fix antimeridian issues for Polygons, geometry={geometry}"
+        )
+    if strategy == Strategy.NORMALIZE:
+        geometry = normalize(geometry)
+        bbox = geometry.bounds
+        item.geometry = shapely.geometry.mapping(geometry)
+        item.bbox = bbox
+    elif strategy == Strategy.SPLIT:
+        split_geometry = split(geometry)
+        if split_geometry:
+            xmin = 180
+            xmax = -180
+            for geom in split_geometry.geoms:
+                if geom.bounds[0] > xmax:
+                    xmax = geom.bounds[0]
+                if geom.bounds[2] < xmin:
+                    xmin = geom.bounds[2]
+            bounds = split_geometry.bounds
+            # https://datatracker.ietf.org/doc/html/rfc7946#section-5.2
+            item.bbox = [xmax, bounds[1], xmin, bounds[3]]
+            item.geometry = shapely.geometry.mapping(split_geometry)
+
+
+def split(polygon: Polygon) -> Optional[MultiPolygon]:
+    """Splits a single WGS84 polygon into a multipolygon across the antimeridian.
+
+    If the polygon does not cross the antimeridian, returns None. Only handles
+    exterior rings (can't handle interior).
+
+    NOTE: Will not work on polygons that enclose the north or south poles.
+    TODO: Fix this
+
+    Args:
+        polygon (shapely.geometry.Polygon): The input polygon.
+
+    Returns:
+        Optional[shapely.geometry.MultiPolygon]: The output polygons, or None if
+            no split occurred.
+    """
+    normalized = normalize(polygon)
+    if normalized.bounds[0] < -180:
+        longitude = -180
+    elif normalized.bounds[2] > 180:
+        longitude = 180
+    else:
+        return None
+    splitter = LineString(((longitude, -90), (longitude, 90)))
+    split = shapely.ops.split(normalized, splitter)
+    if len(split.geoms) < 2:
+        return None
+    geoms = list()
+    for geom in split.geoms:
+        bounds = geom.bounds
+        if bounds[0] < -180:
+            geoms.append(shapely.affinity.translate(geom, xoff=360))
+        elif bounds[2] > 180:
+            geoms.append(shapely.affinity.translate(geom, xoff=-360))
+        else:
+            geoms.append(geom)
+    return MultiPolygon(geoms)
+
+
+def normalize(polygon: Polygon) -> Polygon:
+    """'Normalizes' a WGS84 lat/lon polygon.
+
+    This converts the polygon's x coordinates to all be the same sign, even if
+    the polygon crosses the antimeridian. E.g.:
+
+    ```
+    canonical = Polygon(((170, 40), (170, 50), (-170, 50), (-170, 40), (170, 40)))
+    normalized = stactools.core.utils.antimeridian.normalize(canonical)
+    assert normalized.equals(shapely.geometry.box(170, 40, 190, 50))
+    ```
+
+    Inspired by
+    https://towardsdatascience.com/around-the-world-in-80-lines-crossing-the-antimeridian-with-python-and-shapely-c87c9b6e1513.
+
+    NOTE: Will not work on polygons that enclose the north or south poles.
+    TODO: Fix this
+
+    Args:
+        polygon (shapely.geometry.Polygon): The input polygon.
+
+    Returns:
+        shapely.geometry.Polygon: The normalized polygon.
+    """
+    coords = list(polygon.exterior.coords)
+    for index, (start, end) in enumerate(zip(coords, coords[1:])):
+        delta = end[0] - start[0]
+        if abs(delta) > 180:
+            coords[index + 1] = (end[0] - math.copysign(360, delta), end[1])
+    polygon = Polygon(coords)
+    centroid = polygon.centroid
+    if centroid.x > 180:
+        return shapely.affinity.translate(polygon, xoff=-360)
+    elif centroid.x < -180:
+        return shapely.affinity.translate(polygon, xoff=+360)
+    else:
+        return polygon

tests/core/test_utils.py

@@ -1,10 +1,16 @@
 from contextlib import contextmanager
+import datetime
 import os.path
 from tempfile import TemporaryDirectory
 import unittest
+from pystac import Item
+import pytest
 
 import rasterio
+import shapely.geometry
+from shapely.geometry import Polygon, MultiPolygon
 
+from stactools.core.utils import antimeridian
 from stactools.core.utils.convert import cogify
 from tests import test_data
 
@@ -32,3 +38,126 @@ def test_profile(self):
             with rasterio.open(outfile) as dataset:
                 self.assertEqual(dataset.compression,
                                  rasterio.enums.Compression.lzw)
+
+    def test_antimeridian_split(self) -> None:
+        # From https://datatracker.ietf.org/doc/html/rfc7946#section-3.1.9
+        canonical = Polygon(
+            ((170, 40), (170, 50), (-170, 50), (-170, 40), (170, 40)))
+        split = antimeridian.split(canonical)
+        assert split
+        expected = MultiPolygon((
+            shapely.geometry.box(170, 40, 180, 50),
+            shapely.geometry.box(-180, 40, -170, 50),
+        ))
+        for actual, expected in zip(split.geoms, expected.geoms):
+            assert actual.equals(expected)
+
+        doesnt_cross = Polygon(
+            ((170, 40), (170, 50), (180, 50), (180, 40), (170, 40)))
+        split = antimeridian.split(doesnt_cross)
+        assert split is None
+
+        canonical_other_way = Polygon(
+            ((-170, 40), (170, 40), (170, 50), (-170, 50), (-170, 40)))
+        split = antimeridian.split(canonical_other_way)
+        assert split
+        expected = MultiPolygon((
+            shapely.geometry.box(-180, 40, -170, 50),
+            shapely.geometry.box(170, 40, 180, 50),
+        ))
+        for actual, expected in zip(split.geoms, expected.geoms):
+            assert actual.equals(
+                expected), f"actual={actual}, expected={expected}"
+
+    def test_antimeridian_split_complicated(self) -> None:
+        complicated = Polygon(((170, 40), (170, 50), (-170, 50), (170, 45),
+                               (-170, 40), (170, 40)))
+        split = antimeridian.split(complicated)
+        assert split
+        expected = MultiPolygon((
+            Polygon(((170, 40), (170, 45), (180, 42.5), (180, 40), (170, 40))),
+            Polygon(((170, 45), (170, 50), (180, 50), (180, 47.5), (170, 45))),
+            Polygon(((-180, 50), (-170, 50), (-180, 47.5), (-180, 50))),
+            Polygon(((-180, 42.5), (-170, 40), (-180, 40), (-180, 42.5))),
+        ))
+        for actual, expected in zip(split.geoms, expected.geoms):
+            assert actual.equals(
+                expected), f"actual={actual}, expected={expected}"
+
+    def test_antimeridian_normalize(self) -> None:
+        canonical = Polygon(
+            ((170, 40), (170, 50), (-170, 50), (-170, 40), (170, 40)))
+        normalized = antimeridian.normalize(canonical)
+        expected = shapely.geometry.box(170, 40, 190, 50)
+        assert normalized.equals(
+            expected), f"actual={normalized}, expected={expected}"
+
+        canonical_other_way = Polygon(
+            ((-170, 40), (170, 40), (170, 50), (-170, 50), (-170, 40)))
+        normalized = antimeridian.normalize(canonical_other_way)
+        expected = shapely.geometry.box(-170, 40, -190, 50)
+        assert normalized.equals(
+            expected), f"actual={normalized}, expected={expected}"
+
+    def test_antimeridian_normalize_westerly(self) -> None:
+        westerly = Polygon(
+            ((170, 40), (170, 50), (-140, 50), (-140, 40), (170, 40)))
+        normalized = antimeridian.normalize(westerly)
+        expected = shapely.geometry.box(-190, 40, -140, 50)
+        assert normalized.equals(
+            expected), f"actual={normalized}, expected={expected}"
+
+    def test_antimeridian_normalize_easterly(self) -> None:
+        easterly = Polygon(
+            ((-170, 40), (140, 40), (140, 50), (-170, 50), (-170, 40)))
+        normalized = antimeridian.normalize(easterly)
+        expected = shapely.geometry.box(140, 40, 190, 50)
+        assert normalized.equals(
+            expected), f"actual={normalized}, expected={expected}"
+
+    def test_item_fix_antimeridian_split(self) -> None:
+        canonical = Polygon(
+            ((170, 40), (170, 50), (-170, 50), (-170, 40), (170, 40)))
+        item = Item("an-id",
+                    geometry=shapely.geometry.mapping(canonical),
+                    bbox=canonical.bounds,
+                    datetime=datetime.datetime.now(),
+                    properties={})
+        antimeridian.fix_item(item, antimeridian.Strategy.SPLIT)
+        expected = MultiPolygon((
+            shapely.geometry.box(170, 40, 180, 50),
+            shapely.geometry.box(-180, 40, -170, 50),
+        ))
+        for actual, expected in zip(
+                shapely.geometry.shape(item.geometry).geoms, expected.geoms):
+            assert actual.equals(expected)
+        assert item.bbox == [170, 40, -170, 50]
+
+    def test_item_fix_antimeridian_normalize(self) -> None:
+        canonical = Polygon(
+            ((170, 40), (170, 50), (-170, 50), (-170, 40), (170, 40)))
+        item = Item("an-id",
+                    geometry=shapely.geometry.mapping(canonical),
+                    bbox=canonical.bounds,
+                    datetime=datetime.datetime.now(),
+                    properties={})
+        antimeridian.fix_item(item, antimeridian.Strategy.NORMALIZE)
+        expected = shapely.geometry.box(170, 40, 190, 50)
+        assert shapely.geometry.shape(item.geometry).equals(expected)
+        assert item.bbox
+        assert list(item.bbox) == [170.0, 40.0, 190.0, 50.0]
+
+    def test_item_fix_antimeridian_multipolygon_failure(self) -> None:
+        split = MultiPolygon((
+            shapely.geometry.box(170, 40, 180, 50),
+            shapely.geometry.box(-180, 40, -170, 50),
+        ))
+        item = Item("an-id",
+                    geometry=shapely.geometry.mapping(split),
+                    bbox=split.bounds,
+                    datetime=datetime.datetime.now(),
+                    properties={})
+        with pytest.raises(ValueError):
+            antimeridian.fix_item(item, antimeridian.Strategy.SPLIT)
+        with pytest.raises(ValueError):
+            antimeridian.fix_item(item, antimeridian.Strategy.NORMALIZE)