pyproj icon indicating copy to clipboard operation
pyproj copied to clipboard

DOC: Add CRS.get_geod() to geod docs

Open mwtoews opened this issue 3 years ago • 5 comments

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)

mwtoews avatar Oct 31 '22 01:10 mwtoews

Does CRS.get_geod do what you are looking for?

snowman2 avatar Oct 31 '22 02:10 snowman2

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

mwtoews avatar Oct 31 '22 02:10 mwtoews

Probably work mentioning/linking on https://pyproj4.github.io/pyproj/stable/api/geod.html

I agree. This addition would be welcome.

snowman2 avatar Oct 31 '22 14:10 snowman2

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

snowman2 avatar Oct 31 '22 14:10 snowman2

I'd be happy to take up this task

krishnaglodha avatar Jul 12 '23 12:07 krishnaglodha