Update module to handle both geodetic and cartesian coordinate systems

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)
-```

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",

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
 

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",
+)

geo_extensions/transformations.py → 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: <https://wiki.earthdata.nasa.gov/spaces/CMR/pages/50036858/
+CMR+Data+Partner+User+Guide#CMRDataPartnerUserGuide-CartesianCoordinateSystem>
+
+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.

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
+        ],
+    )