+R_A appears to have no effect in Geod() and CRS()

[gwerbin-tive]

Problem description

The +R_A option is not being honored in Pyproj, whereas it seems to work fine using the geod CLI tool.

This problem happens if I use pyproj.Geod() or pyproj.CRS().get_geod().

Code sample

My Python code:

import logging
import os

os.environ['PROJ_DEBUG'] = '3'
logging.basicConfig(level=logging.DEBUG)

import pyproj
# DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(proj.ini): call fopen(/.../lib/python3.10/site-packages/pyproj/proj_dir/share/proj/proj.ini) - succeeded

geod1 = pyproj.Geod('+proj=lonlat +ellps=WGS84')
print(geod1.f)
# 0.0033528106647474805
geod1.fwd(-70, 40, 60, 10_000)
# (-69.89851793919364, 40.04498640681145, -119.93473805953444)

geod2 = pyproj.Geod('+proj=lonlat +ellps=WGS84 +R_A')
print(geod2.f)
# 0.0033528106647474805
geod2.fwd(-70, 40, 60, 10_000)
# (-69.89851793919364, 40.04498640681145, -119.93473805953444)

Equivalent CLI invocation using geod:

export PROJ_DEBUG=3

geod -f '%0.4f' +proj=lonlat +ellps=WGS84 <<< '-70E 40N 60d 10000m'
# pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
# pj_ellipsoid - final:    ellps=WGS84
# 70.0447	40.2273	-119.7864m

geod -f '%0.4f' +proj=lonlat +ellps=WGS84 +R_A <<< '-70E 40N 60d 10000m'
# pj_ellipsoid - final: a=6371007.181 f=1/  0.000, errno=0
# pj_ellipsoid - final:   R_A ellps=WGS84
# 70.0448	40.2282	-119.7855m

Expected Output

I expected the output from geod1 to be different from geod2, and that geod2.f should be 0.

I was also surprised to see that the output from the geod command line tool was different from the output from pyproj.Geod.fwd. Presumably the difference can be ascribed to different options when compiling the CLI tool (installed via MacPorts) and the Python package (installed via the wheels on PyPI).

I further expected to see "trace" logging output comparable to what the geod CLI tool emitted, but I suppose that's a separate issue.

Environment Information

MacOS 12.4 on ARM64.

Versions

pyproj info:
    pyproj: 3.4.0
      PROJ: 9.1.0
  data dir: /.../lib/python3.10/site-packages/pyproj/proj_dir/share/proj
user_data_dir: /.../.local/share/proj
PROJ DATA (recommended version): 1.11
PROJ Database: 1.2
EPSG Database: v10.074 [2022-08-01]
ESRI Database: ArcGIS Pro 3.0 [2022-07-09]
IGNF Database: 3.1.0 [2019-05-24]

System:
    python: 3.10.6 (main, Sep 19 2022, 14:40:49) [Clang 13.1.6 (clang-1316.0.21.2.5)]
executable: /.../bin/python
   machine: macOS-12.4-arm64-arm-64bit

Python deps:
   certifi: 2022.9.14
    Cython: None
setuptools: 65.4.1
       pip: 22.2.2

Installation method

Installed using Pip in a Venv:

python -m venv venv
./venv/bin/pip install -U pyproj

[snowman2]

https://pyproj4.github.io/pyproj/stable/api/geod.html#pyproj.Geod.__init__

The parameters of the ellipsoid may also be set directly using the ‘a’ (semi-major or equatorial axis radius) keyword, and any one of the following keywords: ‘b’ (semi-minor, or polar axis radius), ‘e’ (eccentricity), ‘es’ (eccentricity squared), ‘f’ (flattening), or ‘rf’ (reciprocal flattening).

This behavior is expected for Geod in it's current state.

A PR is welcome adding support for R_A.

[gwerbin-tive]

Thanks, I missed that bit in the documentation.

Do the same limitations apply to both Geod and CRS? I don't see an equivalent section in the docs for CRS (looking at this page: https://pyproj4.github.io/pyproj/stable/api/crs/crs.html).

I found that CRS('+proj=lonlat +ellps=WGS84 +R_A') produces an ellipsoid attribute that is not a sphere, which unsurprisingly produces a non-spherical Geod:

pyproj/pyproj/crs/crs.py

Lines 503 to 516 in 5723f37

	def get_geod(self) -> Optional[Geod]:
	"""
	Returns
	-------
	pyproj.geod.Geod:
	Geod object based on the ellipsoid.
	"""
	if self.ellipsoid is None:
	return None
	return Geod(
	a=self.ellipsoid.semi_major_metre,
	rf=self.ellipsoid.inverse_flattening,
	b=self.ellipsoid.semi_minor_metre,
	)

Is that a completely unrelated issue?

[snowman2]

I believe that the CRS issue is related to: OSGeo/PROJ#3170

[gwerbin-tive]

I found that issue before and I suspected as much, thank you for clarifying.

That said, why is +R_A being ignored in this case? Is this an "upstream" behavior in Proj itself, or something specific happening in Pyproj?

[snowman2]

That said, why is +R_A being ignored in this case?

This comment might be helpful: #1061 (comment)

Is this an "upstream" behavior in Proj itself, or something specific happening in Pyproj?

For CRS, this is likely due to how PROJ handles it.

[snowman2]

https://proj.org/usage/ellipsoids.html#ellipsoid-spherification-parameters

If size and shape are given as ellps=xxx, later shape and size parameters are are taken into account as modifiers for the built-in ellipsoid definition.

While this may seem strange, it is in accordance with historical PROJ behavior. It can e.g. be used to define coordinates on the ellipsoid scaled to unit semimajor axis by specifying +ellps=xxx +a=1

I am thinking this may be worth raising an issue upstream: https://github.com/OSGeo/PROJ/

>>> from pyproj import CRS
>>> cc = CRS("+proj=lonlat +ellps=WGS84")
>>> cc.ellipsoid
ELLIPSOID["WGS 84",6378137,298.257223563,
    LENGTHUNIT["metre",1],
    ID["EPSG",7030]]
>>> ccr = CRS("+proj=lonlat +ellps=WGS84 +R_A")
>>> ccr.ellipsoid
ELLIPSOID["WGS 84",6378137,298.257223563,
    LENGTHUNIT["metre",1],
    ID["EPSG",7030]]

$projinfo '+proj=lonlat +ellps=WGS84 +type=crs'
PROJ.4 string:
+proj=longlat +ellps=WGS84 +no_defs +type=crs

WKT2:2019 string:
GEOGCRS["unknown",
    DATUM["Unknown based on WGS84 ellipsoid",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1],
            ID["EPSG",7030]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
$ projinfo '+proj=lonlat +ellps=WGS84 +R_A +type=crs'
PROJ.4 string:
+proj=lonlat +ellps=WGS84 +R_A +type=crs

WKT2:2019 string:
GEOGCRS["unknown",
    DATUM["Unknown based on WGS84 ellipsoid",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1],
            ID["EPSG",7030]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
    REMARK["PROJ CRS string: +proj=lonlat +ellps=WGS84 +R_A"]]

[gwerbin-tive]

Thanks @snowman2, I'll see about asking on their mailing list or issue tracker.

[gwerbin-tive]

Oops sorry, I meant to leave this open as an enhancement request for R_A in Geod!