From 67ea5760539ac5bc3e3b37d50ea4f74ccbd5e19e Mon Sep 17 00:00:00 2001 From: HugoFara Date: Sat, 25 Jul 2026 22:54:51 +0200 Subject: [PATCH] BUG: Honour PROJ spherification parameters in Geod (issue #1157) Geod does not hand its init string to PROJ; it parses the ellipsoid out in Python. That parser discarded every parameter without a value: # We can ignore safely any parameter that doesn't have a value if kvpair.find("=") == -1: continue which is where +R_A and the rest of the spherification family went. The result was silently wrong rather than unsupported: Geod('+ellps=WGS84 +R_A') returned ellipsoidal answers to a question asked about a sphere, with no error and no warning, while PROJ's own geod CLI honoured the flag. Collect the spherification parameters before valueless parameters are discarded, and apply PROJ's formulas from src/ell_set.cpp, including the mutual exclusion order and the collapse to a sphere afterwards. All eight now agree with PROJ to the millimetre; the +R_A radius also matches the published WGS84 authalic sphere radius of 6371007.1809 m. Unlike PROJ, the latitude of +R_lat_a/+R_lat_g is accepted only as decimal degrees, not DMS. CRS('+ellps=WGS84 +R_A').get_geod() is unchanged: +R_A has no representation in a CRS, so PROJ drops it there too. --- docs/history.rst | 1 + pyproj/geod.py | 125 ++++++++++++++++++++++++++++++++++++++++++++-- test/test_geod.py | 60 ++++++++++++++++++++++ 3 files changed, 183 insertions(+), 3 deletions(-) diff --git a/docs/history.rst b/docs/history.rst index 48cb1a023..94ad12c84 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -16,6 +16,7 @@ Latest - DOC: Update fiona CRS compatibility for fiona 1.9+ (issue #1360) - 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) +- BUG: Honour the ``+R_A``, ``+R_V``, ``+R_a``, ``+R_g``, ``+R_h``, ``+R_lat_a``, ``+R_lat_g`` and ``+R_C`` spherification parameters in :class:`pyproj.Geod` (issue #1157) - BLD: update build-time dependencies 3.7.2 diff --git a/pyproj/geod.py b/pyproj/geod.py index 4c3d7b335..551461dd3 100644 --- a/pyproj/geod.py +++ b/pyproj/geod.py @@ -56,6 +56,101 @@ def _params_from_ellps_map(ellps: str) -> tuple[float, float, float, float, bool return semi_major_axis, semi_minor_axis, flattening, eccentricity_squared, sphere +# Spherification parameters, in the order PROJ resolves them. They are mutually +# exclusive: PROJ applies the first one present in this order and ignores the +# rest. See ellps_spherification in PROJ src/ell_set.cpp. +_SPHERIFICATION_KEYS = ("R_A", "R_V", "R_a", "R_g", "R_h", "R_lat_a", "R_lat_g", "R_C") + +# Series coefficients for the equal-area and equal-volume spheres, +# from PROJ src/ell_set.cpp. +_SIXTH = 1 / 6.0 +_RA4 = 17 / 360.0 +_RA6 = 67 / 3024.0 +_RV4 = 5 / 72.0 +_RV6 = 55 / 1296.0 + + +def _sphere_radius_from_spherification( + requested: dict[str, str], + semi_major_axis: float, + semi_minor_axis: float, + eccentricity_squared: float, + lat_0: float, +) -> float: + """ + Radius of the sphere requested by a PROJ spherification parameter. + + Parameter + --------- + requested: dict[str, str] + The spherification parameters found in the input, mapped to their + (possibly empty) values. + semi_major_axis: float + The semi-major axis of the ellipsoid to spherify. + semi_minor_axis: float + The semi-minor axis of the ellipsoid to spherify. + eccentricity_squared: float + The eccentricity squared of the ellipsoid to spherify. + lat_0: float + Latitude of origin in degrees, used only by ``+R_C``. + + Returns + ------- + float + + """ + key = next(key for key in _SPHERIFICATION_KEYS if key in requested) + if key == "R_A": + # sphere with the same surface area as the ellipsoid + return semi_major_axis * ( + 1.0 + - eccentricity_squared + * (_SIXTH + eccentricity_squared * (_RA4 + eccentricity_squared * _RA6)) + ) + if key == "R_V": + # sphere with the same volume as the ellipsoid + return semi_major_axis * ( + 1.0 + - eccentricity_squared + * (_SIXTH + eccentricity_squared * (_RV4 + eccentricity_squared * _RV6)) + ) + if key == "R_a": + return (semi_major_axis + semi_minor_axis) / 2 + if key == "R_g": + return math.sqrt(semi_major_axis * semi_minor_axis) + if key == "R_h": + if semi_major_axis + semi_minor_axis == 0: + raise GeodError("Invalid ellipsoid for +R_h: a + b is zero.") + return ( + 2 * semi_major_axis * semi_minor_axis / (semi_major_axis + semi_minor_axis) + ) + # The remaining variants are evaluated at a latitude. PROJ parses that + # latitude with proj_dmstor and so also accepts DMS strings; only decimal + # degrees are supported here. + raw_latitude = lat_0 if key == "R_C" else requested[key] + try: + latitude = float(raw_latitude) + except (TypeError, ValueError) as error: + raise GeodError( + f"Invalid latitude for +{key}: {raw_latitude!r}. " + "Only decimal degrees are supported." + ) from error + if math.fabs(latitude) > 90: + raise GeodError(f"Invalid latitude for +{key}: |{latitude}| should be <= 90.") + sin_latitude = math.sin(math.radians(latitude)) + factor = 1 - eccentricity_squared * sin_latitude**2 + if factor == 0: + raise GeodError(f"Invalid eccentricity for +{key}.") + if key == "R_lat_a": + # arithmetic mean of the ellipsoid radii at the given latitude + return semi_major_axis * ( + (1.0 - eccentricity_squared + factor) / (2 * factor * math.sqrt(factor)) + ) + # geometric mean of the ellipsoid radii at the given latitude (+R_lat_g), + # or the radius of the conformal sphere at lat_0 (+R_C) + return semi_major_axis * math.sqrt(1 - eccentricity_squared) / factor + + def _params_from_kwargs(kwargs: dict) -> tuple[float, float, float, float]: """ Build Geodesic parameters from input kwargs: @@ -198,20 +293,30 @@ def __init__(self, initstring: str | None = None, **kwargs) -> None: # if initparams is a proj-type init string, # convert to dict. ellpsd: dict[str, str | float] = {} + spherification: dict[str, str] = {} if initstring is not None: for kvpair in initstring.split(): + key, _, val = kvpair.partition("=") + key = key.lstrip("+") + # The spherification parameters are the only ones that are + # meaningful without a value, so they have to be collected + # before valueless parameters are discarded. + if key in _SPHERIFICATION_KEYS: + spherification[key] = val + continue # Actually only +a and +b are needed # We can ignore safely any parameter that doesn't have a value - if kvpair.find("=") == -1: + if not val: continue - key, val = kvpair.split("=") - key = key.lstrip("+") if key in ["a", "b", "rf", "f", "es", "e"]: ellpsd[key] = float(val) else: ellpsd[key] = val # merge this dict with kwargs dict. kwargs = dict(list(kwargs.items()) + list(ellpsd.items())) + for key in _SPHERIFICATION_KEYS: + if key in kwargs: + spherification.setdefault(key, kwargs.pop(key)) sphere = False if "ellps" in kwargs: ( @@ -229,6 +334,20 @@ def __init__(self, initstring: str | None = None, **kwargs) -> None: eccentricity_squared, ) = _params_from_kwargs(kwargs) + if spherification: + # PROJ collapses the ellipsoid to the requested sphere and clears + # the remaining ellipsoidal parameters. + semi_major_axis = _sphere_radius_from_spherification( + spherification, + semi_major_axis, + semi_minor_axis, + eccentricity_squared, + float(kwargs.get("lat_0", 0)), + ) + semi_minor_axis = semi_major_axis + flattening = 0.0 + eccentricity_squared = 0.0 + if math.fabs(flattening) < 1.0e-8: sphere = True diff --git a/test/test_geod.py b/test/test_geod.py index 59d51d98d..4264fd9ac 100644 --- a/test/test_geod.py +++ b/test/test_geod.py @@ -8,6 +8,7 @@ from numpy.testing import assert_almost_equal, assert_array_equal from pyproj import Geod +from pyproj.exceptions import GeodError from pyproj.geod import GeodIntermediateFlag, reverse_azimuth try: @@ -953,3 +954,62 @@ def test_geod__reverse_azimuth(radians): xx = xy.T[0] yy = xy.T[1] assert_almost_equal(reverse_azimuth(xx * f, radians=radians), yy * f) + + +# Spherification radii for +ellps=WGS84, read from PROJ's own resolved +# ellipsoid (``PROJ_DEBUG=3 geod +ellps=WGS84 ``, "pj_ellipsoid - final"), +# which is the reference implementation pyproj is expected to agree with. +# The +R_A value is corroborated by the published WGS84 authalic sphere radius +# of 6371007.1809 m, which it matches to 0.2 mm. +@pytest.mark.parametrize( + "spherification, expected_radius", + [ + ("+R_A", 6371007.181), + ("+R_V", 6371000.790), + ("+R_a", 6367444.657), + ("+R_g", 6367435.680), + ("+R_h", 6367426.702), + ("+R_lat_a=45", 6378110.053), + ("+R_lat_g=45", 6378101.030), + ("+R_C", 6356752.314), + ], +) +def test_geod__spherification(spherification, expected_radius): + geod = Geod(f"+proj=lonlat +ellps=WGS84 {spherification}") + assert geod.a == pytest.approx(expected_radius, abs=0.001) + # PROJ collapses the ellipsoid to a sphere, clearing the other parameters. + assert geod.b == pytest.approx(expected_radius, abs=0.001) + assert geod.f == 0 + assert geod.es == 0 + + +def test_geod__spherification__changes_result(): + # https://github.com/pyproj4/pyproj/issues/1157 + # +R_A used to be dropped silently, so these two agreed exactly. + ellipsoidal = Geod("+proj=lonlat +ellps=WGS84") + authalic = Geod("+proj=lonlat +ellps=WGS84 +R_A") + assert authalic.a != ellipsoidal.a + assert authalic.f == 0 + assert authalic.fwd(-70, 40, 60, 10000) != ellipsoidal.fwd(-70, 40, 60, 10000) + + +def test_geod__spherification__kwargs(): + assert Geod(ellps="WGS84", R_A=True).a == pytest.approx(6371007.181, abs=0.001) + + +def test_geod__spherification__is_mutually_exclusive(): + # PROJ applies the first parameter in its own fixed order, not the order + # they appear in the string, and ignores the rest. + assert Geod("+ellps=WGS84 +R_g +R_A").a == pytest.approx(6371007.181, abs=0.001) + assert Geod("+ellps=WGS84 +R_A +R_g").a == pytest.approx(6371007.181, abs=0.001) + + +@pytest.mark.parametrize("spherification", ["+R_lat_a=91", "+R_lat_g=-91"]) +def test_geod__spherification__invalid_latitude(spherification): + with pytest.raises(GeodError, match="should be <= 90"): + Geod(f"+ellps=WGS84 {spherification}") + + +def test_geod__spherification__invalid_latitude_value(): + with pytest.raises(GeodError, match="Only decimal degrees"): + Geod("+ellps=WGS84 +R_lat_a=45d30'")