diff --git a/README.md b/README.md index bb78401..c6218f6 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,30 @@ This library consists of some common functions needed to manipulate polygons before posting them to CMR. This includes functions to split polygons along the antimeridian and perform other 'clean up' operations. +## Beware the Coordinate System + +CMR supports two modes for interpreting how the points of a polygon are supposed +to be connected to each other: Geodetic and Cartesian. + +In the geodetic coordinate system, points are connected to each other by the +shortest arc on a spherical approximation of the Earth. This is the closest +system to how a satellite would actually be seeing the area in the real world. +However, it is not typically the way that polygons would be rendered in a web +search tool as these usually use web mercator projection. + +In the cartesian coordinate system, points are connected to each other more or +less along longitude and latitude lines. This means that points will be +connected by straight lines when viewed in a mercator projection, which is +useful for web tools but not necessarily an accurate representation of the data +acquired by the satellite. + +Keep these differences in mind when implementing a polygon transformation +pipeline, and be aware of what mode the collection you are posting to is set up +in. As a general rule, the higher the point density is on a polygon, the less +the difference between the two coordinate system interpretations will be because +each line segment will be shorter. It is therefore a good idea to start your +pipeline with a 'densify' operation, do the desired transformations, and then +simplify at the end. ## Example Usage @@ -35,9 +59,9 @@ def my_custom_transformation(polygon): transformer = Transformer([ - simplify_polygon(0.1), my_custom_transformation, split_polygon_on_antimeridian_ccw, + simplify_polygon(0.1), ]) final_polygons = transformer.transform([ @@ -47,15 +71,3 @@ final_polygons = transformer.transform([ ]) ]) ``` - -The default transformer performs some standard transformations that are usually -needed. Check the definition for what those transformations are. - -```python -from geo_extensions import default_transformer - - -WKT = "MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)), ((15 5, 40 10, 10 20, 5 10, 15 5)))" - -polygons = default_transformer.from_wkt(WKT) -``` diff --git a/geo_extensions/__init__.py b/geo_extensions/__init__.py index 2c09489..9ef0669 100644 --- a/geo_extensions/__init__.py +++ b/geo_extensions/__init__.py @@ -3,6 +3,7 @@ polygon_crosses_antimeridian_fixed_size, ) from geo_extensions.transformations import ( + densify_polygon, drop_z_coordinate, reverse_polygon, simplify_polygon, @@ -12,14 +13,8 @@ from geo_extensions.transformer import Transformer, to_polygons from geo_extensions.types import Transformation, TransformationResult -default_transformer = Transformer([ - simplify_polygon(0.1), - split_polygon_on_antimeridian_ccw, -]) - - __all__ = ( - "default_transformer", + "densify_polygon", "drop_z_coordinate", "polygon_crosses_antimeridian_ccw", "polygon_crosses_antimeridian_fixed_size", diff --git a/geo_extensions/checks.py b/geo_extensions/checks.py index 60bdbad..e461670 100644 --- a/geo_extensions/checks.py +++ b/geo_extensions/checks.py @@ -1,12 +1,25 @@ +"""Geospatial aware polygon tests. + +When polygons are required to be in counter-clockwise order, this means counter- +clockwise in the real world, on the surface of the Earth. Specifically, polygons +MUST NOT be oriented using the shapely `orient` function, as this function +treats polygons as existing on an infinite flat plane and may end up actually +mis-ordering the polygons. Knowing that a polygon is in fact counter-clockwise +ordered on the surface of the Earth makes the shapely `is_ccw` property a very +useful and easy check to determine if the polygon crosses the antimeridian, as +in this case, the polygon will appear to be mis-ordered in the infinite flat +plane space. +""" + from shapely.geometry import Polygon def polygon_crosses_antimeridian_ccw(polygon: Polygon) -> bool: """Checks if the longitude coordinates 'wrap around' the 180/-180 line. - The polygon must be oriented in counter-clockwise order. - - :param polygon: the polygon to check + :param polygon: the polygon to check, must be known to be in counter- + clockwise order. + :returns: true if the polygon crosses the antimeridian """ # Polygons crossing the antimeridian will appear to be mis-ordered or @@ -25,6 +38,7 @@ def polygon_crosses_antimeridian_fixed_size( :param min_lon_extent: the lower bound for the distance between the longitude values of the bounding box enclosing the entire polygon. Must be between (0, 180) exclusive. + :returns: true if the polygon crosses the antimeridian """ assert 0 < min_lon_extent < 180 diff --git a/geo_extensions/transformations/__init__.py b/geo_extensions/transformations/__init__.py new file mode 100644 index 0000000..a3fe138 --- /dev/null +++ b/geo_extensions/transformations/__init__.py @@ -0,0 +1,52 @@ +"""Geospatial helpers to prepare spatial extents for submission to CMR. + +CMR can interpret polygons as being in one of two coordinate systems, Cartesian +or Geodetic. Transformation helpers that must assume they are working in one of +these two coordinate systems all found in their own submodules `cartesian` and +`geodetic`. + +Each coordinate system also has its own set of constraints that polygons must +follow. These constraints are documented in the corresponding module for the +coordinate system. + +There are some additional constraints that depend on the data format being used. +For UMM-G, polygons must + - Be counter clockwise ordered + - Include closure points + +In this module, polygons are expected to follow the UMM-G constraints. + +A table describing the differences between data format requirements can be found +here: +https://wiki.earthdata.nasa.gov/display/CMR/Polygon+Support+in+CMR+Search+and+Ingest+Interfaces + +There are several challenges with representing polygons on a spherical surface, +the primary being that since all straight lines will 'wrap around' the surface, +it becomes impossible to unabmiguously define a polygon using only an ordered +set of points. This is the primary reason for the CMR requirements, as those +additional constraints make it possible to determine exactly which area is +meant by a set of points. Unfortunately, the polygons that we get from the data +provider won't necessarily meet those same requirements and we must use mission +specific knowledge to convert them to an unambiguous set of polygons to be used +by CMR. + +This module aims to assist in that conversion to unambiguous polygons using the +CMR additional requirements. +""" + +from geo_extensions.transformations.cartesian import ( + simplify_polygon, + split_polygon_on_antimeridian_ccw, + split_polygon_on_antimeridian_fixed_size, +) +from geo_extensions.transformations.general import drop_z_coordinate, reverse_polygon +from geo_extensions.transformations.geodetic import densify_polygon + +__all__ = ( + "densify_polygon", + "drop_z_coordinate", + "reverse_polygon", + "simplify_polygon", + "split_polygon_on_antimeridian_ccw", + "split_polygon_on_antimeridian_fixed_size", +) diff --git a/geo_extensions/transformations.py b/geo_extensions/transformations/cartesian.py similarity index 60% rename from geo_extensions/transformations.py rename to geo_extensions/transformations/cartesian.py index 0e1a993..416abd0 100644 --- a/geo_extensions/transformations.py +++ b/geo_extensions/transformations/cartesian.py @@ -1,41 +1,14 @@ -"""Geospatial helpers to prepare spatial extents for submission to CMR. - -CMR has the following constraints on spatial extents: - - The implemented Geodetic model uses the great circle distance to connect - two vertices for constructing a polygon area or line. If there is not - enough density (that is, the number of points) for a set of vertices, - then the line or the polygon area might be misinterpreted or the - metadata might be considered invalid. - - Any single spatial area may cross the International Date Line and/or Poles - - Any single spatial area may not cover more than one half of the earth. - -Taken from: https://wiki.earthdata.nasa.gov/pages/viewpage.action?spaceKey=CMR&title=CMR+Data+Partner+User+Guide - -There are also additional constraints that depend on the data format being used. -For UMM-G polygons must - - Be counter clockwise ordered - - Include closure points - -A table describing the differences between data format requirements can be found here: -https://wiki.earthdata.nasa.gov/display/CMR/Polygon+Support+in+CMR+Search+and+Ingest+Interfaces - -There are several challenges with representing polygons on a spherical surface, -the primary being that since all straight lines will 'wrap around' the surface, -it becomes impossible to unabmiguously define a polygon using only an ordered -set of points. This is the primary reason for the CMR requirements, as those -additional constraints make it possible to determine exactly which area is -meant by a set of points. Unfortunately, the polygons that we get from the data -provider won't necessarily meet those same requirements and we must use mission -specific knowledge to convert them to an unambiguous set of polygons to be used -by CMR. - -This module aims to assist in that conversion to unambiguous polygons using the -CMR additional requirements. Any polygons passed in as arguments or returned -from functions in this module are assumed to be in counter clockwise order as -seen in the spherical space. This makes detecting whether a polygon crosses the -antimeridian (wraps around the edge of the flat coordinate system) very easy -even in the general case, because such a polygon will appear to be clockwise -ordered in the shapely flat space. +"""Geospatial helpers for working in a cartesian coordinate system. + +CMR has the following constraints for cartesian polygons: + - Any single spatial area may not cross the International Date Line (unless + it is a bounding box) or Poles. + - Two vertices will be connected with a straight line. + +Taken from: + +This module contains helpers to fulfill the cartesian system CMR requirements. """ from typing import cast @@ -50,17 +23,15 @@ ) from geo_extensions.types import Transformation, TransformationResult -Point = tuple[float, float] -Bbox = list[Point] - ANTIMERIDIAN = LineString([(180, 90), (180, -90)]) def simplify_polygon(tolerance: float, preserve_topology: bool = True) -> Transformation: - """Create a transformation that calls polygon.simplify. + """CARTESIAN: Create a transformation that calls polygon.simplify. :returns: a callable transformation using the passed parameters """ + def simplify(polygon: Polygon) -> TransformationResult: """Perform a shapely simplify operation on the polygon.""" # NOTE(reweeden): I have been unable to produce a situation where a @@ -76,21 +47,9 @@ def simplify(polygon: Polygon) -> TransformationResult: return simplify -def reverse_polygon(polygon: Polygon) -> TransformationResult: - """Perform a shapely reverse operation on the polygon.""" - yield polygon.reverse() - - -def drop_z_coordinate(polygon: Polygon) -> TransformationResult: - yield Polygon( - (x, y) - for x, y, *_ in polygon.exterior.coords - ) - - def split_polygon_on_antimeridian_ccw(polygon: Polygon) -> TransformationResult: - """Perform adjustment when the polygon crosses the antimeridian and is known - to be wound in counter clockwise order. + """CARTESIAN: Perform adjustment when the polygon crosses the antimeridian + and is known to be wound in counter clockwise order. CMR requires the polygon to be split into two separate polygons to avoid it being interpreted as wrapping the long way around the Earth. @@ -116,8 +75,8 @@ def split_polygon_on_antimeridian_ccw(polygon: Polygon) -> TransformationResult: def split_polygon_on_antimeridian_fixed_size( min_lon_extent: float, ) -> Transformation: - """Perform adjustment when the polygon crosses the antimeridian using a - heuristic to determine if the polygon needs to be split. + """CARTESIAN: Perform adjustment when the polygon crosses the antimeridian + using a heuristic to determine if the polygon needs to be split. CMR requires the polygon to be split into two separate polygons to avoid it being interpreted as wrapping the long way around the Earth. diff --git a/geo_extensions/transformations/general.py b/geo_extensions/transformations/general.py new file mode 100644 index 0000000..e01eb73 --- /dev/null +++ b/geo_extensions/transformations/general.py @@ -0,0 +1,23 @@ +"""Geospatial helpers that work the same regardless of which coordinate system +the polygons are using. +""" + +from shapely.geometry import Polygon + +from geo_extensions.types import TransformationResult + + +def reverse_polygon(polygon: Polygon) -> TransformationResult: + """Perform a shapely reverse operation on the polygon.""" + yield polygon.reverse() + + +def drop_z_coordinate(polygon: Polygon) -> TransformationResult: + """Drop the third element from each coordinate in the polygon.""" + yield Polygon( + shell=((x, y) for x, y, *_ in polygon.exterior.coords), + holes=[ + ((x, y) for x, y, *_ in interior.coords) + for interior in polygon.interiors + ], + ) diff --git a/geo_extensions/transformations/geodetic.py b/geo_extensions/transformations/geodetic.py new file mode 100644 index 0000000..d2aaf7f --- /dev/null +++ b/geo_extensions/transformations/geodetic.py @@ -0,0 +1,99 @@ +"""Geospatial helpers for working in a geodetic coordinate system. + +CMR has the following constraints for geodetic polygons: + - The implemented Geodetic model uses the great circle distance to connect + two vertices for constructing a polygon area or line. If there is not + enough density (that is, the number of points) for a set of vertices, + then the line or the polygon area might be misinterpreted or the + metadata might be considered invalid. + - Any single spatial area may cross the International Date Line and/or Poles. + - Any single spatial area may not cover more than one half of the earth. + +Taken from: + +This module contains helpers to fulfill the geodetic system CMR requirements. +""" + +import itertools +from collections.abc import Generator +from typing import TypeVar + +from pygeodesy.sphericalTrigonometry import LatLon +from shapely.coords import CoordinateSequence +from shapely.geometry import Polygon + +from geo_extensions.types import Transformation, TransformationResult + +T = TypeVar("T") + + +def densify_polygon(tolerance_meters: float) -> Transformation: + """GEODETIC: Create a transformation that increases the point density of a + polygon along great circle arcs between each point. + + In a sense, this is an opposite operation from 'simplify'. + + :param tolerance_meters: The maximum allowable cross track error between + a line segment when interpreted as a cartesian point. Must be greater + than 0. + :returns: a callable transformation using the passed parameters + """ + if tolerance_meters <= 0: + raise ValueError("'tolerance_meters' must be greater than 0") + + def densify(polygon: Polygon) -> TransformationResult: + """Densify the polygon by adding additional points along the great + circle arcs between the existing points. + """ + yield Polygon( + shell=_densify_ring(polygon.exterior.coords, tolerance_meters), + holes=[ + _densify_ring(interior.coords, tolerance_meters) + for interior in polygon.interiors + ], + ) + + return densify + + +def _densify_ring( + coords: CoordinateSequence, + tolerance_meters: float, +) -> Generator[tuple[float, ...]]: + assert tolerance_meters > 0 + + if len(coords) < 2: + yield from coords + return + + for c1, c2 in itertools.pairwise(coords): + p1 = _shapely_to_pygeodesy(c1) + p2 = _shapely_to_pygeodesy(c2) + + yield c1 + + for p_new in _densify_edge(p1, p2, tolerance_meters): + yield (p_new.lon, p_new.lat) + + yield c2 + + +def _densify_edge(p1: LatLon, p2: LatLon, tolerance_meters: float) -> Generator[LatLon]: + # Cartesian midpoint + c_mid_cartesian = ((p1.lon + p2.lon) / 2, (p1.lat + p2.lat) / 2) + p_mid_cartesian = _shapely_to_pygeodesy(c_mid_cartesian) + + error_meters = abs(p_mid_cartesian.crossTrackDistanceTo(p1, p2)) + if error_meters < tolerance_meters: + return + + # Add a point in the middle and recursively densify the resulting edges + p_mid = p1.midpointTo(p2) + yield from _densify_edge(p1, p_mid, tolerance_meters) + yield p_mid + yield from _densify_edge(p_mid, p2, tolerance_meters) + + +def _shapely_to_pygeodesy(coord: tuple[float, ...]) -> LatLon: + return LatLon(coord[1], coord[0]) diff --git a/pyproject.toml b/pyproject.toml index 5fa2296..137bb1a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.core.masonry.api" [tool.poetry] name = "geo-extensions" -version = "0.3.1" +version = "0.4.0" description = "" authors = ["Rohan Weeden "] license = "APACHE-2" @@ -14,6 +14,7 @@ packages = [{include = "geo_extensions"}] [tool.poetry.dependencies] python = "^3.10" +pygeodesy = "^25.9.9" shapely = "^2.0.3" [tool.poetry.group.dev.dependencies] diff --git a/tests/test_geo_extensions.py b/tests/test_geo_extensions.py index 7b31b25..90bca8e 100644 --- a/tests/test_geo_extensions.py +++ b/tests/test_geo_extensions.py @@ -1,12 +1,21 @@ from shapely.geometry import Polygon -from geo_extensions import default_transformer +from geo_extensions import ( + Transformer, + simplify_polygon, + split_polygon_on_antimeridian_ccw, +) def test_default_transformer(data_path): wkt_path = data_path / "OPERA_L2_RTC-S1_T114-243299-IW1_20230722T122534Z_20230818T184635Z_S1A_30_v0.4.wkt" - polygons = default_transformer.from_wkt(wkt_path.read_text()) + transformer = Transformer([ + split_polygon_on_antimeridian_ccw, + simplify_polygon(0.1), + ]) + + polygons = transformer.from_wkt(wkt_path.read_text()) assert polygons == [ Polygon([ (-64.18242781701073, 80.92318071697005), diff --git a/tests/test_transformations.py b/tests/test_transformations.py index 7df1a51..5afc007 100644 --- a/tests/test_transformations.py +++ b/tests/test_transformations.py @@ -1,3 +1,4 @@ +import pytest import shapely.geometry import strategies from hypothesis import HealthCheck, given, settings @@ -5,6 +6,7 @@ from shapely.geometry import Polygon from geo_extensions.transformations import ( + densify_polygon, drop_z_coordinate, simplify_polygon, split_polygon_on_antimeridian_ccw, @@ -54,6 +56,104 @@ def test_simplify_line(): ] +def test_densify(): + polygon = Polygon([ + (50, 75), + (10, 80), + (0, 77), + (40, 70), + (50, 75), + ]) + + assert list(densify_polygon(50_000)(polygon)) == [ + Polygon([ + (50, 75), + (34.100003241169595, 78.2028289318241), + (10, 80), + (0, 77), + (24.297878219303588, 74.40374356383884), + (40, 70), + (50, 75), + ]), + ] + + +def test_densify_with_holes(): + polygon = Polygon( + shell=[ + (50, 70), + (50, 80), + (0, 80), + (0, 70), + (50, 70), + ], + holes=[ + [ + (45, 72), + (45, 78), + (5, 78), + (5, 72), + (45, 72), + ], + ], + ) + + assert list(densify_polygon(50_000)(polygon)) == [ + Polygon( + shell=[ + (50, 70), + (50, 80), + (24.999999999999996, 80.92053252671789), + (0, 80), + (0, 70), + (25, 71.7438759890997), + (50, 70), + ], + holes=[ + [ + (45, 72), + (45, 78), + (25, 78.70451161084236), + (5, 78), + (5, 72), + (25, 73.02127815507072), + (45, 72), + ], + ], + ), + ] + + +def test_densify_idempotent(): + polygon = Polygon([ + (50, 75), + (10, 80), + (0, 77), + (40, 70), + (50, 75), + ]) + + transformation = densify_polygon(50_000) + + densified_polygons = list(transformation(polygon)) + double_densified_polygons = [ + poly + for densified_polygon in densified_polygons + for poly in transformation(densified_polygon) + ] + + assert densified_polygons == double_densified_polygons + + +def test_densify_incomplete(): + assert list(densify_polygon(50_000)(Polygon())) == [Polygon()] + + +def test_densify_error(): + with pytest.raises(ValueError, match="must be greater than 0"): + densify_polygon(0) + + def test_drop_z_coordinate(): polygon = Polygon([ (180, 1, 10), @@ -84,6 +184,37 @@ def test_drop_z_coordinate_noop(): assert list(drop_z_coordinate(polygon)) == [polygon] +def test_drop_z_coordinate_holes(): + polygon = Polygon( + shell=[ + (100, 10, 10), + (100, 0, 10), + (80, 0, 10), + (80, 10, 10), + (100, 10, 10), + ], + holes=[ + [(93, 8, 10), (83, 8, 10), (83, 2, 10), (93, 8, 10)], + [(97, 2, 10), (97, 8, 10), (87, 2, 10), (97, 2, 10)], + ], + ) + assert list(drop_z_coordinate(polygon)) == [ + Polygon( + shell=[ + (100, 10), + (100, 0), + (80, 0), + (80, 10), + (100, 10), + ], + holes=[ + [(93, 8), (83, 8), (83, 2), (93, 8)], + [(97, 2), (97, 8), (87, 2), (97, 2)], + ], + ), + ] + + @given(polygon=strategies.rectangles()) @settings(suppress_health_check=[HealthCheck.filter_too_much]) def test_split_polygon_on_antimeridian_ccw_returns_ccw(polygon):