diff --git a/pyproject.toml b/pyproject.toml index 930d2964..31725b33 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,7 @@ requires-python = ">=3.8" dependencies = [ "Shapely>=1.8.5.post1", "aiohttp>=3.8.3", + "antimeridian>=0.2.3", "click>=8.1.3", "fsspec>=2021.7", "lxml>=4.9.2", diff --git a/src/stactools/core/utils/antimeridian.py b/src/stactools/core/utils/antimeridian.py index ef4918f5..bcecb9fd 100644 --- a/src/stactools/core/utils/antimeridian.py +++ b/src/stactools/core/utils/antimeridian.py @@ -1,17 +1,22 @@ -"""`Antimeridian `_ utilities.""" +"""`Antimeridian `_ utilities. + +Most of the functionality in this module is implemented in `antimeridian +`_. +""" from __future__ import annotations import math -from dataclasses import dataclass +import warnings from enum import Enum, auto -from typing import List, Optional, Tuple +from typing import Optional +import antimeridian import shapely.affinity import shapely.geometry import shapely.ops from pystac import Item -from shapely.geometry import LineString, MultiPolygon, Polygon +from shapely.geometry import MultiPolygon, Polygon class Strategy(Enum): @@ -39,16 +44,16 @@ def fix_item(item: Item, strategy: Strategy) -> Item: Item: The input item, whether it was modified or not. """ geometry = shapely.geometry.shape(item.geometry) - if isinstance(geometry, Polygon): - multi_polygon = False - elif isinstance(geometry, MultiPolygon): - multi_polygon = True - else: - raise ValueError( - "Can only fix antimeridian issues for Polygons or MultiPolygons, " - f"geometry={geometry}" - ) if strategy == Strategy.NORMALIZE: + if isinstance(geometry, Polygon): + multi_polygon = False + elif isinstance(geometry, MultiPolygon): + multi_polygon = True + else: + raise ValueError( + "Can only fix antimeridian issues for Polygons or MultiPolygons, " + f"geometry={geometry}" + ) if multi_polygon: normalized_geometry = normalize_multipolygon(geometry) else: @@ -58,22 +63,11 @@ def fix_item(item: Item, strategy: Strategy) -> Item: item.geometry = shapely.geometry.mapping(normalized_geometry) item.bbox = bbox elif strategy == Strategy.SPLIT: - if multi_polygon: - split_geometry = split_multipolygon(geometry) - else: - 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) + fixed = shapely.geometry.shape(antimeridian.fix_shape(geometry)) + item.bbox = fixed.bounds + item.geometry = shapely.geometry.mapping(fixed) + else: + raise NotImplementedError(f"Unknown strategy: {strategy}") return item @@ -81,14 +75,9 @@ 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 + .. deprecated:: v0.4.7 + Use the `antimeridian `_ package + instead. Args: polygon (:py:class:`shapely.geometry.Polygon`): The input polygon. @@ -97,44 +86,28 @@ def split(polygon: Polygon) -> Optional[MultiPolygon]: Optional[:py:class:`shapely.geometry.MultiPolygon`]: The output polygons, or None if no split occurred. """ - normalized = normalize(polygon) - if normalized is None: - return None - if normalized.bounds[0] < -180: - longitude = -180 - elif normalized.bounds[2] > 180: - longitude = 180 + warnings.warn( + ( + "`stactools.core.utils.antimeridian.split` is deprecated. " + "To fix polygons that cross the antimeridian, use the antimeridian " + "package (https://pypi.org/project/antimeridian/)." + ), + DeprecationWarning, + ) + fixed = antimeridian.fix_polygon(polygon) + if fixed.geom_type == "MultiPolygon": + return fixed 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: - geom = shapely.geometry.polygon.orient(geom) - 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 split_multipolygon(multi_polygon: MultiPolygon) -> Optional[MultiPolygon]: """Splits multiple WGS84 polygons into a multipolygon across the antimeridian. - If none of the contained polygons 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 + .. deprecated:: v0.4.7 + Use the `antimeridian `_ package + instead. Args: multi_polygon (:py:class:`shapely.geometry.MultiPolygon`): The input @@ -142,17 +115,19 @@ def split_multipolygon(multi_polygon: MultiPolygon) -> Optional[MultiPolygon]: Returns: Optional[:py:class:`shapely.geometry.MultiPolygon`]: - The output polygons, or None if no split occurred. + The output polygons. Will not return None, but the Optional return + type is kept to preserve backwards compatibility until this function + is removed. """ - polygons = [] - for polygon in multi_polygon.geoms: - split_polygon = split(polygon) - if split_polygon: - polygons.extend(split_polygon.geoms) - if polygons: - return MultiPolygon(polygons) - else: - return None + warnings.warn( + ( + "`stactools.core.utils.antimeridian.split_multipolygon` is deprecated. " + "To fix polygons that cross the antimeridian, use the antimeridian " + "package (https://pypi.org/project/antimeridian/)." + ), + DeprecationWarning, + ) + return antimeridian.fix_multi_polygon(multi_polygon) def normalize(polygon: Polygon) -> Optional[Polygon]: @@ -177,12 +152,24 @@ def normalize(polygon: Polygon) -> Optional[Polygon]: Todo: Fix this + .. deprecated:: v0.4.7 + "Normalization" does not conform to the GeoJSON specification, and its + use is discouraged. + Args: polygon (:py:class:`shapely.geometry.Polygon`): The input polygon. Returns: Optional[:py:class:`shapely.geometry.Polygon`]: The normalized polygon. """ + warnings.warn( + ( + "`stactools.core.utils.antimeridian.normalize` is deprecated. " + "Normalization does not conform to the GeoJSON spec and its use is " + "discouraged." + ), + DeprecationWarning, + ) coords = list(polygon.exterior.coords) has_changes = False for index, (start, end) in enumerate(zip(coords, coords[1:])): @@ -218,6 +205,10 @@ def normalize_multipolygon(multi_polygon: MultiPolygon) -> Optional[MultiPolygon Todo: Fix this + .. deprecated:: v0.4.7 + "Normalization" does not conform to the GeoJSON specification, and its + use is discouraged. + Args: multi_polygon (:py:class:`shapely.geometry.MultiPolygon`): The input multi-polygon. @@ -226,6 +217,14 @@ def normalize_multipolygon(multi_polygon: MultiPolygon) -> Optional[MultiPolygon Optional[:py:class:`shapely.geometry.MultiPolygon`]: The normalized multi-polygon. """ + warnings.warn( + ( + "`stactools.core.utils.antimeridian.normalize_multipolygon` is deprecated. " + "Normalization does not conform to the GeoJSON spec and its use is " + "discouraged." + ), + DeprecationWarning, + ) polygons = list() changes_made = False for polygon in multi_polygon.geoms: @@ -248,7 +247,9 @@ def enclose_poles(polygon: Polygon) -> Polygon: the geometry up to the north (or down to the south) pole. This is useful for (e.g.) polar-orbiting satellites who have swaths that go over the poles. - This function will raise a value error if the polygon has any interior rings. + .. deprecated:: v0.4.7 + Use the `antimeridian `_ package + instead. Args: polygon (:py:class:`shapely.geometry.Polygon`): An input polygon. @@ -256,96 +257,21 @@ def enclose_poles(polygon: Polygon) -> Polygon: Returns: :py:class:`shapely.geometry.Polygon`: The same polygon, modified to enclose the poles. - """ - if bool(polygon.interiors): - raise ValueError("cannot enclose poles if there is an interior ring") - coords = list(polygon.exterior.coords) - crossings = list() - - # First pass is to detect all antimeridian crossings, without actually - # modifying the coordinates. This is to protect against the case when there - # are additional crossings in between the pole crossings. - for i, (start, end) in enumerate(zip(coords, coords[1:])): - longitude_delta = end[0] - start[0] - if abs(longitude_delta) > 180: - crossings.append(_Crossing.from_points(i, start, end)) - - # We only want the southernmost and northernmost crossings. - crossings = sorted(crossings, key=lambda c: c.latitude) - if len(crossings) > 2: - crossings = [crossings[0], crossings[-1]] - - # If there are two crossings, we know the southernmost one is around the - # south pole and the northernmost is around the north pole, even if the - # crossing latitude is in the other hemisphere. - # - # If we only have one crossing, we just guess that it's crossing on the - # closer pole. - if len(crossings) == 2: - crossings[0].north_pole = False - crossings[1].north_pole = True - - # We work from the back of the list so we can use the indices. - crossings = sorted(crossings, key=lambda c: c.index, reverse=True) - for crossing in crossings: - coords = ( - coords[0 : crossing.index + 1] - + crossing.enclosure() - + coords[crossing.index + 1 :] - ) - - return Polygon(coords) - - -@dataclass -class _Crossing: - index: int - latitude: float - positive_to_negative: bool - north_pole: bool - - @classmethod - def from_points( - cls, index: int, start: Tuple[float, float], end: Tuple[float, float] - ) -> _Crossing: - latitude_delta = end[1] - start[1] - if start[0] > 0: - split_latitude = round( - start[1] - + (180.0 - start[0]) * latitude_delta / (end[0] + 360.0 - start[0]), - 7, - ) - return cls( - index=index, - latitude=split_latitude, - positive_to_negative=True, - north_pole=split_latitude > 0, - ) - else: - split_latitude = round( - start[1] - + (start[0] + 180.0) * latitude_delta / (start[0] + 360.0 - end[0]), - 7, - ) - return cls( - index=index, - latitude=split_latitude, - positive_to_negative=False, - north_pole=split_latitude > 0, - ) - def enclosure(self) -> List[Tuple[float, float]]: - if self.positive_to_negative: - longitudes = (180.0, -180.0) - else: - longitudes = (-180.0, 180.0) - if self.north_pole: - pole_latitude = 90.0 - else: - pole_latitude = -90.0 - return [ - (longitudes[0], self.latitude), - (longitudes[0], pole_latitude), - (longitudes[1], pole_latitude), - (longitudes[1], self.latitude), - ] + Raises: + ValueError: Raised if the polygon was split. This is to keep the return + type the same until this function is removed. + """ + warnings.warn( + ( + "`stactools.core.utils.antimeridian.enclose_poles` is deprecated. " + "To fix polygons that enclose the poles, use the antimeridian " + "package (https://pypi.org/project/antimeridian/)." + ), + DeprecationWarning, + ) + fixed = antimeridian.fix_polygon(polygon) + if fixed.geom_type == "Polygon": + return fixed + else: + raise ValueError("Input polygon was split and we cannot return a MultiPolygon") diff --git a/tests/core/utils/test_antimeridian.py b/tests/core/utils/test_antimeridian.py index c0c2b87c..05cb30db 100644 --- a/tests/core/utils/test_antimeridian.py +++ b/tests/core/utils/test_antimeridian.py @@ -1,5 +1,6 @@ import datetime +import pytest import shapely.geometry from pystac import Item from shapely.geometry import MultiPolygon, Polygon @@ -8,8 +9,9 @@ def test_antimeridian_split() -> 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) + canonical = Polygon(((170, 40), (-170, 40), (-170, 50), (170, 50), (170, 40))) + with pytest.warns(DeprecationWarning): + split = antimeridian.split(canonical) assert split expected = MultiPolygon( ( @@ -22,13 +24,15 @@ def test_antimeridian_split() -> None: assert actual.equals(expected) doesnt_cross = Polygon(((170, 40), (170, 50), (180, 50), (180, 40), (170, 40))) - split = antimeridian.split(doesnt_cross) + with pytest.warns(DeprecationWarning): + split = antimeridian.split(doesnt_cross) assert split is None canonical_other_way = Polygon( - ((-170, 40), (170, 40), (170, 50), (-170, 50), (-170, 40)), + ((-170, 40), (-170, 50), (170, 50), (170, 40), (-170, 40)), ) - split = antimeridian.split(canonical_other_way) + with pytest.warns(DeprecationWarning): + split = antimeridian.split(canonical_other_way) assert split expected = MultiPolygon( ( @@ -41,44 +45,10 @@ def test_antimeridian_split() -> None: assert actual.equals(expected), f"actual={actual}, expected={expected}" -def test_antimeridian_split_complicated() -> None: - complicated = Polygon( - ((170, 40), (170, 50), (-170, 50), (170, 45), (-170, 40), (170, 40)), - ) - split = antimeridian.split(complicated) - assert split - expected = MultiPolygon( - ( - Polygon( - [ - (180.0, 40.0), - (180.0, 42.5), - (170.0, 45.0), - (170.0, 40.0), - (180.0, 40.0), - ], - ), - Polygon([(-180.0, 42.5), (-180.0, 40.0), (-170.0, 40.0), (-180.0, 42.5)]), - Polygon( - [ - (180.0, 47.5), - (180.0, 50.0), - (170.0, 50.0), - (170.0, 45.0), - (180.0, 47.5), - ], - ), - Polygon([(-180.0, 50.0), (-180.0, 47.5), (-170.0, 50.0), (-180.0, 50.0)]), - ), - ) - for actual, expected in zip(split.geoms, expected.geoms): - assert actual.exterior.is_ccw - assert actual.equals(expected), f"actual={actual}, expected={expected}" - - def test_antimeridian_normalize() -> None: canonical = Polygon(((170, 40), (170, 50), (-170, 50), (-170, 40), (170, 40))) - normalized = antimeridian.normalize(canonical) + with pytest.warns(DeprecationWarning): + normalized = antimeridian.normalize(canonical) assert normalized assert normalized.exterior.is_ccw expected = shapely.geometry.box(170, 40, 190, 50) @@ -87,7 +57,8 @@ def test_antimeridian_normalize() -> None: canonical_other_way = Polygon( ((-170, 40), (170, 40), (170, 50), (-170, 50), (-170, 40)), ) - normalized = antimeridian.normalize(canonical_other_way) + with pytest.warns(DeprecationWarning): + normalized = antimeridian.normalize(canonical_other_way) assert normalized assert normalized.exterior.is_ccw expected = shapely.geometry.box(-170, 40, -190, 50) @@ -96,7 +67,8 @@ def test_antimeridian_normalize() -> None: def test_antimeridian_normalize_westerly() -> None: westerly = Polygon(((170, 40), (170, 50), (-140, 50), (-140, 40), (170, 40))) - normalized = antimeridian.normalize(westerly) + with pytest.warns(DeprecationWarning): + normalized = antimeridian.normalize(westerly) assert normalized assert normalized.exterior.is_ccw expected = shapely.geometry.box(-170, 40, -190, 50) @@ -106,7 +78,8 @@ def test_antimeridian_normalize_westerly() -> None: def test_antimeridian_normalize_easterly() -> None: easterly = Polygon(((-170, 40), (140, 40), (140, 50), (-170, 50), (-170, 40))) - normalized = antimeridian.normalize(easterly) + with pytest.warns(DeprecationWarning): + normalized = antimeridian.normalize(easterly) assert normalized assert normalized.exterior.is_ccw expected = shapely.geometry.box(-170, 40, -190, 50) @@ -115,7 +88,7 @@ def test_antimeridian_normalize_easterly() -> None: def test_item_fix_antimeridian_split() -> None: - canonical = Polygon(((170, 40), (170, 50), (-170, 50), (-170, 40), (170, 40))) + canonical = Polygon(((170, 40), (-170, 40), (-170, 50), (170, 50), (170, 40))) item = Item( "an-id", geometry=shapely.geometry.mapping(canonical), @@ -135,7 +108,7 @@ def test_item_fix_antimeridian_split() -> None: expected.geoms, ): assert actual.equals(expected) - assert fix.bbox == [170, 40, -170, 50] + assert fix.bbox == (-180.0, 40.0, 180.0, 50.0) def test_item_fix_antimeridian_normalize() -> None: @@ -147,7 +120,8 @@ def test_item_fix_antimeridian_normalize() -> None: datetime=datetime.datetime.now(), properties={}, ) - fix = antimeridian.fix_item(item, antimeridian.Strategy.NORMALIZE) + with pytest.warns(DeprecationWarning): + fix = antimeridian.fix_item(item, antimeridian.Strategy.NORMALIZE) expected = shapely.geometry.box(170, 40, 190, 50) assert shapely.geometry.shape(fix.geometry).equals(expected) assert fix.bbox @@ -169,17 +143,19 @@ def test_item_fix_antimeridian_multipolygon_ok() -> None: properties={}, ) antimeridian.fix_item(item, antimeridian.Strategy.SPLIT) - antimeridian.fix_item(item, antimeridian.Strategy.NORMALIZE) + with pytest.warns(DeprecationWarning): + antimeridian.fix_item(item, antimeridian.Strategy.NORMALIZE) def test_antimeridian_multipolygon() -> None: multi_polygon = MultiPolygon( [ - Polygon(((170, 40), (170, 42), (-170, 42), (-170, 40), (170, 40))), - Polygon(((170, 48), (170, 50), (-170, 50), (-170, 48), (170, 48))), + Polygon(((170, 40), (-170, 40), (-170, 42), (170, 42), (170, 40))), + Polygon(((170, 48), (-170, 48), (-170, 50), (170, 50), (170, 48))), ], ) - split = antimeridian.split_multipolygon(multi_polygon) + with pytest.warns(DeprecationWarning): + split = antimeridian.split_multipolygon(multi_polygon) assert split expected = MultiPolygon( ( @@ -193,7 +169,8 @@ def test_antimeridian_multipolygon() -> None: assert actual.exterior.is_ccw assert actual.equals(expected), f"actual={actual}, expected={expected}" - normalized = antimeridian.normalize_multipolygon(multi_polygon) + with pytest.warns(DeprecationWarning): + normalized = antimeridian.normalize_multipolygon(multi_polygon) assert normalized expected = MultiPolygon( ( @@ -204,72 +181,3 @@ def test_antimeridian_multipolygon() -> None: for actual, expected in zip(normalized.geoms, expected.geoms): assert actual.exterior.is_ccw assert actual.equals(expected), f"actual={actual}, expected={expected}" - - -def test_antimeridian_enclose_poles() -> None: - before = Polygon(((170, 40), (-170, 50), (-170, -50), (170, -40), (170, 40))) - after = antimeridian.enclose_poles(before) - assert after == Polygon( - ( - (170, 40), - (180, 45), - (180, 90), - (-180, 90), - (-180, 45), - (-170, 50), - (-170, -50), - (-180, -45), - (-180, -90), - (180, -90), - (180, -45), - (170, -40), - (170, 40), - ) - ) - - -def test_antimeridian_enclose_poles_extra_crossing() -> None: - before = Polygon( - ((170, 40), (-170, 50), (-170, -50), (170, -40), (-175, 0), (170, 40)) - ) - after = antimeridian.enclose_poles(before) - assert after == Polygon( - ( - (170, 40), - (180, 45), - (180, 90), - (-180, 90), - (-180, 45), - (-170, 50), - (-170, -50), - (-180, -45), - (-180, -90), - (180, -90), - (180, -45), - (170, -40), - (-175, 0), - (170, 40), - ) - ) - - -def test_antimeridian_enclose_poles_both_northern_hemisphere() -> None: - before = Polygon(((170, 80), (-170, 80), (-170, 10), (170, 10), (170, 80))) - after = antimeridian.enclose_poles(before) - assert after == Polygon( - ( - (170, 80), - (180, 80), - (180, 90), - (-180, 90), - (-180, 80), - (-170, 80), - (-170, 10), - (-180, 10), - (-180, -90), - (180, -90), - (180, 10), - (170, 10), - (170, 80), - ) - )