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?