From b6e606345d23a5450864155decdd22c0c4183387 Mon Sep 17 00:00:00 2001 From: hksamm Date: Sun, 26 Jul 2026 08:08:51 +0900 Subject: [PATCH 1/2] DOC: describe geometry_area_perimeter size limit concretely The note said only that the function "only works with areas up to half the size of the globe" and that "certain large polygons may return negative values", which does not say what "large" or "certain" mean. Both are the same rule: the area is reduced modulo the total area of the ellipsoid into (-A/2, A/2]. A region covering more than half the ellipsoid therefore comes back as the negated area of its complement -- a polygon enclosing 99% of the globe returns about -1%, not +99%. So "large" means precisely "more than half the ellipsoid" (about 2.55e14 m^2 on WGS84), and shape, proximity to a pole and crossing the antimeridian do not matter on their own. Also cross-reference the existing traversal-direction warning, since that is the other, size-independent reason a result can be negative. Adds a doctest for the wraparound and for recovering the enclosed area. Co-Authored-By: Claude Opus 5 --- docs/history.rst | 1 + pyproj/geod.py | 40 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/docs/history.rst b/docs/history.rst index 48cb1a023..145d64681 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -14,6 +14,7 @@ Latest - ENH: Added :class:`pyproj.enums.GridAvailabilityUse` enum for TransformerGroup ``grid_check`` kwarg - ENH: Added always_xy kwarg to :meth:`pyproj.transformer.Transformer.from_pipeline` - DOC: Update fiona CRS compatibility for fiona 1.9+ (issue #1360) +- DOC: Describe the size limit and negative results of :meth:`pyproj.Geod.geometry_area_perimeter` concretely (issue #1592) - BUG: Clear CONTEXT_THREAD_KEY when destroying PJ_CONTEXT (pull #1541) - BUG: Default skew angle for CF grid mapping oblique mercator to 90 (issue #1506) - BLD: update build-time dependencies diff --git a/pyproj/geod.py b/pyproj/geod.py index 4c3d7b335..c72bd3b41 100644 --- a/pyproj/geod.py +++ b/pyproj/geod.py @@ -1031,15 +1031,24 @@ def geometry_area_perimeter( .. note:: lats should be in the range [-90 deg, 90 deg]. - .. note:: | There are a few limitations : - | - only works with areas up to half the size of the globe ; - | - certain large polygons may return negative values. + .. note:: The area is reduced modulo the total area of the ellipsoid, + into the range ``(-A/2, A/2]``, where ``A`` is that total area + (about ``5.10e14`` m\\ :sup:`2` on WGS84, so ``A/2`` is about + ``2.55e14`` m\\ :sup:`2`). A region covering more than half the + ellipsoid is therefore reported as the negated area of its + complement: a polygon enclosing 99% of the globe returns about + -1% of the total area, not +99%. Size is the only reason a + correctly oriented polygon comes back negative, and "large" + means precisely "more than half the ellipsoid" -- shape, + proximity to a pole and crossing the antimeridian do not + matter on their own. .. warning:: The area returned is signed with counter-clockwise (CCW) traversal being treated as positive. For polygons, holes should use the opposite traversal to the exterior (if the exterior is CCW, the holes/interiors should be CW). You can use `shapely.ops.orient` to - modify the orientation. + modify the orientation. This is the other reason a result can be + negative, and unlike the size limit above it applies at any size. If it is a Polygon, it will return the area and exterior perimeter. It will subtract the area of the interior holes. @@ -1063,6 +1072,29 @@ def geometry_area_perimeter( >>> f"{poly_area:.0f} {poly_perimeter:.0f}" '944373881400 3979008' + A polygon enclosing more than half the ellipsoid wraps around to a + negative value, even though its vertices are ordered counter-clockwise: + + >>> nearly_global = Polygon( + ... LineString([ + ... Point(-179, -89), Point(179, -89), Point(179, 89), Point(-179, 89) + ... ]) + ... ) + >>> big_area, _ = geod.geometry_area_perimeter(nearly_global) + >>> big_area < 0 + True + + Its magnitude is the area of the thin strip left outside the polygon, so + adding the total area of the ellipsoid recovers the enclosed area: + + >>> import math + >>> b = geod.a * (1 - geod.f) + >>> ecc2 = 1 - (b * b) / (geod.a * geod.a) + >>> ecc = math.sqrt(ecc2) + >>> total = 2 * math.pi * geod.a**2 * (1 + (1 - ecc2) / ecc * math.atanh(ecc)) + >>> f"{100 * (big_area + total) / total:.1f}%" + '99.4%' + Parameters ---------- From 0db4de48cdcc6556d4b049fc0216bc264f6eb1fb Mon Sep 17 00:00:00 2001 From: hksamm Date: Mon, 27 Jul 2026 17:18:58 +0900 Subject: [PATCH 2/2] Fix doctest failure count for the added shapely examples test_doctest_wrapper asserts an exact number of expected doctest failures, which is 6 when shapely is missing because the examples in Geod.geometry_length and Geod.geometry_area_perimeter cannot build their geometries. The examples added here are shapely-based too, so the count rose to 10 and the assertion failed on every environment without shapely -- the Docker python 3.12, 3.13 and 3.14 jobs. The 3.11 job has shapely, expects 0 failures, and passed, which confirms the examples themselves are correct. Drop the ellipsoid-area arithmetic, which was a lot of docstring for one point, and keep the sign demonstration. That takes the addition from four failures to two, so the missing-shapely count becomes 8. Verified against a source build with PROJ 9.8.1: with shapely absent the wrapper passes at 8, with shapely present it passes at 0, and unmodified main still passes at 6. Co-Authored-By: Claude Opus 5 --- pyproj/geod.py | 13 ++----------- test/test_doctest_wrapper.py | 5 +++-- 2 files changed, 5 insertions(+), 13 deletions(-) diff --git a/pyproj/geod.py b/pyproj/geod.py index c72bd3b41..caf188ed1 100644 --- a/pyproj/geod.py +++ b/pyproj/geod.py @@ -1080,20 +1080,11 @@ def geometry_area_perimeter( ... Point(-179, -89), Point(179, -89), Point(179, 89), Point(-179, 89) ... ]) ... ) - >>> big_area, _ = geod.geometry_area_perimeter(nearly_global) - >>> big_area < 0 + >>> geod.geometry_area_perimeter(nearly_global)[0] < 0 True Its magnitude is the area of the thin strip left outside the polygon, so - adding the total area of the ellipsoid recovers the enclosed area: - - >>> import math - >>> b = geod.a * (1 - geod.f) - >>> ecc2 = 1 - (b * b) / (geod.a * geod.a) - >>> ecc = math.sqrt(ecc2) - >>> total = 2 * math.pi * geod.a**2 * (1 + (1 - ecc2) / ecc * math.atanh(ecc)) - >>> f"{100 * (big_area + total) / total:.1f}%" - '99.4%' + adding the total area of the ellipsoid recovers the enclosed area. Parameters diff --git a/test/test_doctest_wrapper.py b/test/test_doctest_wrapper.py index ce2a92cfb..e6a43fbdf 100644 --- a/test/test_doctest_wrapper.py +++ b/test/test_doctest_wrapper.py @@ -29,8 +29,9 @@ def test_doctests(): try: import shapely.geometry # noqa: F401 pylint: disable=unused-import except (ImportError, OSError): - # missing shapely - expected_failure_count = 6 + # missing shapely; these are the examples in Geod.geometry_length and + # Geod.geometry_area_perimeter that need a shapely geometry + expected_failure_count = 8 # if the below line fails, doctests have failed assert failure_count == expected_failure_count, (