Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DOC: Add CRS.get_geod() to geod docs #1170

Open
mwtoews opened this issue Oct 31, 2022 · 5 comments
Open

DOC: Add CRS.get_geod() to geod docs #1170

mwtoews opened this issue Oct 31, 2022 · 5 comments
Labels
documentation Docs need updating good-first-issue Great issue to get started contributing.

Comments

@mwtoews
Copy link
Contributor

mwtoews commented Oct 31, 2022

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 mwtoews added the proposal Idea for a new feature. label Oct 31, 2022
@snowman2
Copy link
Member

Does CRS.get_geod do what you are looking for?

@mwtoews
Copy link
Contributor Author

mwtoews commented Oct 31, 2022

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:

pyproj/pyproj/crs/crs.py

Lines 512 to 516 in 2a62e36

return Geod(
a=self.ellipsoid.semi_major_metre,
rf=self.ellipsoid.inverse_flattening,
b=self.ellipsoid.semi_minor_metre,
)

@snowman2
Copy link
Member

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

I agree. This addition would be welcome.

@snowman2
Copy link
Member

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 snowman2 added documentation Docs need updating good-first-issue Great issue to get started contributing. and removed proposal Idea for a new feature. labels Nov 2, 2022
@snowman2 snowman2 changed the title ENH: add Geod.from_ellipsoid classmethod DOC: Add CRS.get_geod() to geod docs Apr 18, 2023
@krishnaglodha
Copy link

I'd be happy to take up this task

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
documentation Docs need updating good-first-issue Great issue to get started contributing.
Projects
None yet
Development

No branches or pull requests

3 participants