From 9e4c449d3b7641f6f4c71ed831767ef5ad34e35b Mon Sep 17 00:00:00 2001
From: Rohan Weeden <reweeden@alaska.edu>
Date: Mon, 6 Oct 2025 12:58:28 -0400
Subject: [PATCH 1/7] Refactor transformations module to split out cartesian

---
 geo_extensions/transformations/__init__.py    | 50 +++++++++++++
 .../cartesian.py}                             | 75 +++++--------------
 geo_extensions/transformations/general.py     | 20 +++++
 3 files changed, 87 insertions(+), 58 deletions(-)
 create mode 100644 geo_extensions/transformations/__init__.py
 rename geo_extensions/{transformations.py => transformations/cartesian.py} (60%)
 create mode 100644 geo_extensions/transformations/general.py

diff --git a/geo_extensions/transformations/__init__.py b/geo_extensions/transformations/__init__.py
new file mode 100644
index 0000000..c25f548
--- /dev/null
+++ b/geo_extensions/transformations/__init__.py
@@ -0,0 +1,50 @@
+"""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
+
+__all__ = (
+    "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: <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.
diff --git a/geo_extensions/transformations/general.py b/geo_extensions/transformations/general.py
new file mode 100644
index 0000000..25411ee
--- /dev/null
+++ b/geo_extensions/transformations/general.py
@@ -0,0 +1,20 @@
+"""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(
+        (x, y)
+        for x, y, *_ in polygon.exterior.coords
+    )

From c87fa6a88836587a15cf72d0774c169f1bb7fae3 Mon Sep 17 00:00:00 2001
From: Rohan Weeden <reweeden@alaska.edu>
Date: Mon, 6 Oct 2025 13:00:01 -0400
Subject: [PATCH 2/7] Add function to densify polygon along geodetic lines

---
 geo_extensions/__init__.py                 |  2 +
 geo_extensions/transformations/__init__.py |  2 +
 geo_extensions/transformations/geodetic.py | 97 ++++++++++++++++++++++
 pyproject.toml                             |  1 +
 tests/test_transformations.py              | 94 +++++++++++++++++++++
 5 files changed, 196 insertions(+)
 create mode 100644 geo_extensions/transformations/geodetic.py

diff --git a/geo_extensions/__init__.py b/geo_extensions/__init__.py
index 2c09489..319cf1a 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,
@@ -20,6 +21,7 @@
 
 __all__ = (
     "default_transformer",
+    "densify_polygon",
     "drop_z_coordinate",
     "polygon_crosses_antimeridian_ccw",
     "polygon_crosses_antimeridian_fixed_size",
diff --git a/geo_extensions/transformations/__init__.py b/geo_extensions/transformations/__init__.py
index c25f548..a3fe138 100644
--- a/geo_extensions/transformations/__init__.py
+++ b/geo_extensions/transformations/__init__.py
@@ -40,8 +40,10 @@
     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",
diff --git a/geo_extensions/transformations/geodetic.py b/geo_extensions/transformations/geodetic.py
new file mode 100644
index 0000000..5a082b1
--- /dev/null
+++ b/geo_extensions/transformations/geodetic.py
@@ -0,0 +1,97 @@
+"""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: <https://wiki.earthdata.nasa.gov/spaces/CMR/pages/50036858/
+CMR+Data+Partner+User+Guide#CMRDataPartnerUserGuide-GeodeticCoordinateSystem>
+
+This module contains helpers to fulfill the geodetic system CMR requirements.
+"""
+
+import itertools
+import math
+from collections.abc import Generator
+from typing import TypeVar, cast
+
+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(max_edge_len_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 max_edge_len_meters: The maximum desired length for any line segment.
+        Must be greater than 0.
+    :returns: a callable transformation using the passed parameters
+    """
+    if max_edge_len_meters <= 0:
+        raise ValueError("'max_edge_len_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, max_edge_len_meters),
+            holes=[
+                _densify_ring(interior.coords, max_edge_len_meters)
+                for interior in polygon.interiors
+            ],
+        )
+
+    return densify
+
+
+def _densify_ring(
+    coords: CoordinateSequence,
+    max_edge_len_meters: float,
+) -> Generator[tuple[float, ...]]:
+    assert max_edge_len_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
+
+        distance = p1.distanceTo(p2)
+        if distance < max_edge_len_meters:
+            continue
+
+        segments_needed = math.ceil(distance / max_edge_len_meters)
+        points_needed = segments_needed - 1
+        fraction = 1 / segments_needed
+
+        # TODO(reweeden): Safeguard for low edge lengths that would generate
+        # 'absurdly large' amounts of points?
+        for i in range(points_needed):
+            p_new = cast(
+                LatLon,
+                p1.intermediateTo(p2, fraction=fraction * (i + 1)),
+            )
+            yield (p_new.lon, p_new.lat)
+
+    yield c2
+
+
+def _shapely_to_pygeodesy(coord: tuple[float, ...]) -> LatLon:
+    return LatLon(coord[1], coord[0])
diff --git a/pyproject.toml b/pyproject.toml
index 5fa2296..b6ba8cb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -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_transformations.py b/tests/test_transformations.py
index 7df1a51..34610b6 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,98 @@ def test_simplify_line():
     ]
 
 
+def test_densify():
+    polygon = Polygon([
+        (50, 70),
+        (50, 80),
+        (0, 80),
+        (0, 70),
+        (50, 70),
+    ])
+
+    assert list(densify_polygon(500000)(polygon)) == [
+        Polygon([
+            (50, 70),
+            (50, 73.33333333333333),
+            (50, 76.66666666666667),
+            (50, 80),
+            (25, 80.92053252671789),
+            (0, 80),
+            (0, 76.66666666666667),
+            (0, 73.33333333333333),
+            (0, 70),
+            (11.94262888488008, 71.29264444333815),
+            (24.999999999999996, 71.7438759890997),
+            (38.05737111511993, 71.29264444333815),
+            (50, 70),
+        ]),
+    ]
+
+
+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(500000)(polygon)) == [
+        Polygon(
+            shell=[
+                (50, 70),
+                (50, 73.33333333333333),
+                (50, 76.66666666666667),
+                (50, 80),
+                (25, 80.92053252671789),
+                (0, 80),
+                (0, 76.66666666666667),
+                (0, 73.33333333333333),
+                (0, 70),
+                (11.94262888488008, 71.29264444333815),
+                (24.999999999999996, 71.7438759890997),
+                (38.05737111511993, 71.29264444333815),
+                (50, 70),
+            ],
+            holes=[
+                [
+                    (45, 72),
+                    (45, 75),
+                    (45, 78),
+                    (25, 78.70451161084236),
+                    (5, 78),
+                    (4.999999999999999, 75),
+                    (5, 72),
+                    (18.1052765936037, 72.90478732278757),
+                    (31.894723406396295, 72.90478732278757),
+                    (45, 72),
+                ],
+            ],
+        ),
+    ]
+
+
+def test_densify_incomplete():
+    assert list(densify_polygon(500000)(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),

From c596762b8e29e16d828ec4d7f577628f92102b56 Mon Sep 17 00:00:00 2001
From: Rohan Weeden <reweeden@alaska.edu>
Date: Mon, 6 Oct 2025 13:00:24 -0400
Subject: [PATCH 3/7] Remove default_transformer

---
 README.md                    | 12 ------------
 geo_extensions/__init__.py   |  7 -------
 tests/test_geo_extensions.py | 13 +++++++++++--
 3 files changed, 11 insertions(+), 21 deletions(-)

diff --git a/README.md b/README.md
index bb78401..ac277e4 100644
--- a/README.md
+++ b/README.md
@@ -47,15 +47,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 319cf1a..9ef0669 100644
--- a/geo_extensions/__init__.py
+++ b/geo_extensions/__init__.py
@@ -13,14 +13,7 @@
 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",
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),

From 3b0237f1ccdb45a056422aa8752223b1af1b871d Mon Sep 17 00:00:00 2001
From: Rohan Weeden <reweeden@alaska.edu>
Date: Mon, 6 Oct 2025 13:00:59 -0400
Subject: [PATCH 4/7] Update documentation

---
 README.md                | 26 +++++++++++++++++++++++++-
 geo_extensions/checks.py | 20 +++++++++++++++++---
 2 files changed, 42 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index ac277e4..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([
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
 

From b25b1a5c7e10dd129d4fd65079a41c835cad8a2b Mon Sep 17 00:00:00 2001
From: Rohan Weeden <reweeden@alaska.edu>
Date: Mon, 6 Oct 2025 13:01:06 -0400
Subject: [PATCH 5/7] Bump version number

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index b6ba8cb..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 <reweeden@alaska.edu>"]
 license = "APACHE-2"

From 4fbc34a11ee539bad70bfa1580642a9482dfcfd7 Mon Sep 17 00:00:00 2001
From: Rohan Weeden <reweeden@alaska.edu>
Date: Fri, 10 Oct 2025 18:11:40 -0400
Subject: [PATCH 6/7] Densify using crossTrackDistanceTo to determine error
 from cartesian interpretation

---
 geo_extensions/transformations/geodetic.py | 54 ++++++++--------
 tests/test_transformations.py              | 72 ++++++++++++----------
 2 files changed, 67 insertions(+), 59 deletions(-)

diff --git a/geo_extensions/transformations/geodetic.py b/geo_extensions/transformations/geodetic.py
index 5a082b1..d2aaf7f 100644
--- a/geo_extensions/transformations/geodetic.py
+++ b/geo_extensions/transformations/geodetic.py
@@ -16,9 +16,8 @@
 """
 
 import itertools
-import math
 from collections.abc import Generator
-from typing import TypeVar, cast
+from typing import TypeVar
 
 from pygeodesy.sphericalTrigonometry import LatLon
 from shapely.coords import CoordinateSequence
@@ -29,27 +28,28 @@
 T = TypeVar("T")
 
 
-def densify_polygon(max_edge_len_meters: float) -> Transformation:
+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 max_edge_len_meters: The maximum desired length for any line segment.
-        Must be greater than 0.
+    :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 max_edge_len_meters <= 0:
-        raise ValueError("'max_edge_len_meters' must be greater than 0")
+    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, max_edge_len_meters),
+            shell=_densify_ring(polygon.exterior.coords, tolerance_meters),
             holes=[
-                _densify_ring(interior.coords, max_edge_len_meters)
+                _densify_ring(interior.coords, tolerance_meters)
                 for interior in polygon.interiors
             ],
         )
@@ -59,9 +59,9 @@ def densify(polygon: Polygon) -> TransformationResult:
 
 def _densify_ring(
     coords: CoordinateSequence,
-    max_edge_len_meters: float,
+    tolerance_meters: float,
 ) -> Generator[tuple[float, ...]]:
-    assert max_edge_len_meters > 0
+    assert tolerance_meters > 0
 
     if len(coords) < 2:
         yield from coords
@@ -73,25 +73,27 @@ def _densify_ring(
 
         yield c1
 
-        distance = p1.distanceTo(p2)
-        if distance < max_edge_len_meters:
-            continue
-
-        segments_needed = math.ceil(distance / max_edge_len_meters)
-        points_needed = segments_needed - 1
-        fraction = 1 / segments_needed
-
-        # TODO(reweeden): Safeguard for low edge lengths that would generate
-        # 'absurdly large' amounts of points?
-        for i in range(points_needed):
-            p_new = cast(
-                LatLon,
-                p1.intermediateTo(p2, fraction=fraction * (i + 1)),
-            )
+        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/tests/test_transformations.py b/tests/test_transformations.py
index 34610b6..f7b2e3a 100644
--- a/tests/test_transformations.py
+++ b/tests/test_transformations.py
@@ -58,28 +58,22 @@ def test_simplify_line():
 
 def test_densify():
     polygon = Polygon([
-        (50, 70),
-        (50, 80),
-        (0, 80),
-        (0, 70),
-        (50, 70),
+        (50, 75),
+        (10, 80),
+        (0, 77),
+        (40, 70),
+        (50, 75),
     ])
 
-    assert list(densify_polygon(500000)(polygon)) == [
+    assert list(densify_polygon(50_000)(polygon)) == [
         Polygon([
-            (50, 70),
-            (50, 73.33333333333333),
-            (50, 76.66666666666667),
-            (50, 80),
-            (25, 80.92053252671789),
-            (0, 80),
-            (0, 76.66666666666667),
-            (0, 73.33333333333333),
-            (0, 70),
-            (11.94262888488008, 71.29264444333815),
-            (24.999999999999996, 71.7438759890997),
-            (38.05737111511993, 71.29264444333815),
-            (50, 70),
+            (50, 75),
+            (34.100003241169595, 78.2028289318241),
+            (10, 80),
+            (0, 77),
+            (24.297878219303588, 74.40374356383884),
+            (40, 70),
+            (50, 75),
         ]),
     ]
 
@@ -104,34 +98,25 @@ def test_densify_with_holes():
         ],
     )
 
-    assert list(densify_polygon(500000)(polygon)) == [
+    assert list(densify_polygon(50_000)(polygon)) == [
         Polygon(
             shell=[
                 (50, 70),
-                (50, 73.33333333333333),
-                (50, 76.66666666666667),
                 (50, 80),
-                (25, 80.92053252671789),
+                (24.999999999999996, 80.92053252671789),
                 (0, 80),
-                (0, 76.66666666666667),
-                (0, 73.33333333333333),
                 (0, 70),
-                (11.94262888488008, 71.29264444333815),
-                (24.999999999999996, 71.7438759890997),
-                (38.05737111511993, 71.29264444333815),
+                (25, 71.7438759890997),
                 (50, 70),
             ],
             holes=[
                 [
                     (45, 72),
-                    (45, 75),
                     (45, 78),
                     (25, 78.70451161084236),
                     (5, 78),
-                    (4.999999999999999, 75),
                     (5, 72),
-                    (18.1052765936037, 72.90478732278757),
-                    (31.894723406396295, 72.90478732278757),
+                    (25, 73.02127815507072),
                     (45, 72),
                 ],
             ],
@@ -139,8 +124,29 @@ def test_densify_with_holes():
     ]
 
 
+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(500000)(Polygon())) == [Polygon()]
+    assert list(densify_polygon(50_000)(Polygon())) == [Polygon()]
 
 
 def test_densify_error():

From 7e9482ba689ff45d98b48edab5bb1d6d5e5514dc Mon Sep 17 00:00:00 2001
From: Rohan Weeden <reweeden@alaska.edu>
Date: Thu, 18 Dec 2025 19:49:29 -0500
Subject: [PATCH 7/7] Update `drop_z_coordinate` to work with holey polygons
 too

---
 geo_extensions/transformations/general.py |  7 +++--
 tests/test_transformations.py             | 31 +++++++++++++++++++++++
 2 files changed, 36 insertions(+), 2 deletions(-)

diff --git a/geo_extensions/transformations/general.py b/geo_extensions/transformations/general.py
index 25411ee..e01eb73 100644
--- a/geo_extensions/transformations/general.py
+++ b/geo_extensions/transformations/general.py
@@ -15,6 +15,9 @@ def reverse_polygon(polygon: Polygon) -> TransformationResult:
 def drop_z_coordinate(polygon: Polygon) -> TransformationResult:
     """Drop the third element from each coordinate in the polygon."""
     yield Polygon(
-        (x, y)
-        for x, y, *_ in polygon.exterior.coords
+        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/tests/test_transformations.py b/tests/test_transformations.py
index f7b2e3a..5afc007 100644
--- a/tests/test_transformations.py
+++ b/tests/test_transformations.py
@@ -184,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):