BUG: Honour PROJ spherification parameters in Geod (issue #1157)

[HugoFara]

Geod parses the ellipsoid out of its init string in Python rather than handing the string to PROJ, and that parser dropped every parameter without a value:

# We can ignore safely any parameter that doesn't have a value
if kvpair.find("=") == -1:
    continue

+R_A and the rest of the spherification family went through there. The result was silently wrong rather than unsupported — Geod('+ellps=WGS84 +R_A') returned ellipsoidal answers with no error and no warning, while PROJ's own geod CLI honoured the flag.

This collects those parameters before valueless ones are discarded and applies the formulas from PROJ src/ell_set.cpp (ellps_spherification), including its mutual-exclusion order and the collapse to a sphere afterwards.

Semi-major axis for +ellps=WGS84, against PROJ_DEBUG=3 geod ...:

parameter	before	after	PROJ
+R_A	6378137.000	6371007.181	6371007.181
+R_V	6378137.000	6371000.790	6371000.790
+R_a	6378137.000	6367444.657	6367444.657
+R_g	6378137.000	6367435.680	6367435.680
+R_h	6378137.000	6367426.702	6367426.702
+R_lat_a=45	6378137.000	6378110.053	6378110.053
+R_lat_g=45	6378137.000	6378101.030	6378101.030
+R_C	6378137.000	6356752.314	6356752.314

Geod.fwd now reproduces the geod CLI to 7 decimal places on all of them. The +R_A radius also matches the published WGS84 authalic sphere radius (6371007.1809 m) to 0.2 mm, so the expected values do not rest on PROJ alone.

Two deliberate limits:

• PROJ parses the +R_lat_a/+R_lat_g latitude with proj_dmstor and so accepts DMS; this accepts decimal degrees only and raises GeodError otherwise.
• CRS('+ellps=WGS84 +R_A').get_geod() is unchanged. +R_A has no representation in a CRS, and PROJ drops it there too, so honouring it would diverge from PROJ rather than match it.

Verification: 14 new tests, all failing before the change and passing after. Full suite 1047 passed, mypy clean, pre-commit clean. Built against PROJ 9.8.1.

• Closes +R_A appears to have no effect in Geod() and CRS() #1157
• Tests added
• Fully documented, including history.rst for all changes and api/*.rst for new API

[jjimenezshaw]

Hi @HugoFara
Reading the PR description I thought that it was passing the parameters like +R_A to PROJ. However I see that it is re-implementing all the logic in pyproj. That is not the idea of pyproj I guess.
I think it should use the PROJ functions pj_ell_set or pj_ellipsoid.

[HugoFara]

Hi @jjimenezshaw . I couldn't find a way to reach either function from pyproj, maybe you can help me here.

Both are declared in proj_internal.h, not proj.h, and that header isn't installed. PROJ ships five: HEADERS_LIBPROJ in src/lib_proj.cmake is proj.h, proj_experimental.h, proj_constants.h, proj_symbol_rename.h, geodesic.h (plus apps/projapps_lib.h). I checked the actual proj-devel 9.8.1 package and proj_internal.h isn't in it.

Only pj_ell_set is exported as _Z10pj_ell_setP6pj_ctxP8ARG_listPdS3_

pj_ellipsoid isn't callable from outside libproj at all. pj_ell_set has PROJ_DLL but no extern "C", so what's actually exported is a C++-mangled name rather than stable API (should we build an internal paralist?).

Have I missed a supported route? If pyproj is expected to build against proj_internal.h, I'll happily redo it that way.

I also notived there is a way to let PROJ do the work through the public API only. +proj=cart applies spherification like any other PJ, so the ellipsoid can be read back by transforming two points: (0, 0, 0) gives x = a, and (0, 90, 0) gives z = b:

params	a	b
+ellps=WGS84	6378137.0000	6356752.3142
+ellps=WGS84 +R_A	6371007.1811	6371007.1811
+ellps=WGS84 +R_lat_a=45	6378110.0529	6378110.0529
+ellps=bessel	6377397.1550	6356078.9628
+a=6378137 +rf=298.257223563	6378137.0000	6356752.3142

That covers named ellipsoids, explicit +a/+rf and all eight spherification variants, with PROJ doing every bit of the math. It would remove the formulas in this PR along with the hand-rolled proj-string parser and the pj_ellps lookup in _params_from_ellps_map, so it fixes the class of bug rather than this one instance, and DMS latitudes for +R_lat_a/+R_lat_g, which this PR doesn't support, would come for free.

The trade-off is that it's a much bigger change than #1157 warrants: it moves Geod.__init__ into the Cython layer and touches every construction path, and it's a behaviour change for anyone relying on the current parsing quirks.

Would you prefer I take this PR that way, or keep it narrow and open a separate issue for the delegation? Happy either way.

[jjimenezshaw]

You are right. It is only exposed in the C++ API.

I don't know what would be the best option. Should we expose something in the C API, @rouault @kbevers ? Is there another function that we don't see?

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.95%. Comparing base (e81e5e5) to head (cd59643).
⚠️ Report is 52 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1610      +/-   ##
==========================================
+ Coverage   95.85%   95.95%   +0.09%     
==========================================
  Files          20       20              
  Lines        1880     1926      +46     
==========================================
+ Hits         1802     1848      +46     
  Misses         78       78              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[HugoFara]

Hi, I saw you marked the work as planned @snowman2

I updated the branch, closed the coverage gap and fixed failing tests. Let me know if you'd need anything else.

[snowman2]

The code looks good from an overview and tests standpoint. I agree with @jjimenezshaw about the logic living in PROJ in the long term. However, I don't see a function API for converting the geod string to a geod object.

@HugoFara are you comfortable enough in the PROJ codebase to attempt adding a public method to create a geod object from a geod string?

[HugoFara]

I can try it. I'll require the help of my Claude to explain the specifics of the codebase, thanks for trusting it with me @snowman2 .