I am using pyproj inverse transform to add azimuth and distance "info" to an ordered geodataframe (gdf), but my datasets are in different locations around the world. I need to use a local UTM EPSG to get accurate azimuths and distances (e.g., discussion here; this is common knowledge).

For a given EPSG, how can I systematically extract the g = pyproj.Geod(ellps='X') info from the CRS in the right format for X?

Below is my best attempt using myellipsoid = CRS.from_user_input(myepsg).ellipsoid, but it's in the wrong format. In this example, "GRS 1980" needs to be "GRS80"

%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj
from pyproj import CRS

myid = [1, 1, 1]
myorder = [1, 2, 3]
lat = [5174925.07851924, 5174890.26832387, 5174855.45812849]
long = [1521631.6994673, 1521667.11033893, 1521702.52121056]
myepsg = 2193

df = pd.DataFrame(list(zip(myid, myorder, lat, long)), columns =['myid', 'myorder', 'lat', 'long']) 
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
gdf_pt = gdf_pt.set_crs(epsg=myepsg)
myellipsoid = CRS.from_user_input(myepsg).ellipsoid


ax = gdf_pt.plot();
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);

g = pyproj.Geod(ellps=myellipsoid) 
for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = g.inv(gdf_pt.long[i], gdf_pt.lat[i], gdf_pt.long[i-1], gdf_pt.lat[i-1])
    gdf_pt.loc[i, 'az_fwd'] = myinfo[0]
    gdf_pt.loc[i, 'az_back'] = myinfo[1]
    gdf_pt.loc[i, 'dist'] = myinfo[2]
    gdf_pt.loc[i, 'bearing'] = max(myinfo[1], myinfo[0])


enter image description here

Using: Windows 10; conda 4.8.2; Python 3.8.3; shapely 1.7.0 py38hbf43935_3 conda-forge; pyproj 2.6.1.post1 py38h1dd9442_0 conda-forge


1 Answer 1


You can use CRS.get_geod:

Here is an example:

from pyproj import CRS

crs = CRS.from_epsg(epsg_code)
geod = crs.get_geod()
  • print(geod) yields Geod('+a=6378137 +f=0.0033528106811822724'), which throws an error in g = pyproj.Geod(ellps=geod): TypeError: unhashable type: 'Geod'
    – a11
    Jul 2, 2020 at 20:37
  • 1
    geod is the pyproj.Geod object. No need to pass it in.
    – snowman2
    Jul 3, 2020 at 1:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.