pyproj
pyproj copied to clipboard
DOC: Add CRS.get_geod() to geod docs
It would be handy to create a Geod object from an Ellipsoid object. E.g., this is how I would expect:
from pyproj import CRS, Geod
crs = CRS.from_epsg(4326)
geod = Geod.from_ellipsoid(crs.ellipsoid)
or is that too convoluted? Is there a simpler way? A few current working methods:
Geod(ellps='WGS84')
# Geod(ellps='WGS84')
Geod(a=crs.ellipsoid.semi_major_metre, b=crs.ellipsoid.semi_minor_metre)
# Geod('+a=6378137 +f=0.0033528106647475126')
geod = Geod(a=crs.ellipsoid.semi_major_metre, f=1/crs.ellipsoid.inverse_flattening)
# Geod(ellps='WGS84')
and a bonus curiosity, probably a floating-point precision error
assert Geod(ellps='WGS84') == Geod(a=crs.ellipsoid.semi_major_metre, f=1/crs.ellipsoid.inverse_flattening)
assert Geod(ellps='WGS84') != Geod(a=crs.ellipsoid.semi_major_metre, b=crs.ellipsoid.semi_minor_metre)
Does CRS.get_geod do what you are looking for?
It does! Probably work mentioning/linking on https://pyproj4.github.io/pyproj/stable/api/geod.html
Also, it exhibits the same oddity:
crs = CRS.from_epsg(4326)
crs.get_geod()
# Geod('+a=6378137 +f=0.0033528106647475126')
assert crs.get_geod() != Geod(ellps='WGS84')
assert crs.get_geod() == Geod(a=crs.ellipsoid.semi_major_metre, b=crs.ellipsoid.semi_minor_metre)
I'd have to take a look under the hood, but perhaps the constructor is using b instead of 1/f, as the expected result would be Geod(ellps='WGS84').
Update: here it is: https://github.com/pyproj4/pyproj/blob/2a62e36ac5209c8c739587875b4b03959d27f58b/pyproj/crs/crs.py#L512-L516
Probably work mentioning/linking on https://pyproj4.github.io/pyproj/stable/api/geod.html
I agree. This addition would be welcome.
It does look like a precision difference:
>>> from pyproj import pj_ellps
>>> pj_ellps["WGS84"]
{'a': 6378137.0, 'rf': 298.257223563, 'description': 'WGS 84'}
>>> from pyproj import CRS
>>> cc = CRS("WGS84")
>>> geod = cc.get_geod()
>>> geod.a
6378137.0
>>> geod.b
6356752.314245179
>>> geod.f
0.0033528106647475126
>>> from pyproj import Geod
>>> wg = Geod(ellps="WGS84")
>>> wg.a
6378137.0
>>> wg.b
6356752.314245179
>>> wg.f
0.0033528106647474805
I'd be happy to take up this task