Pyproj inverse transform returns NaN with EPSG that is not 4326

Question

I am using pyproj inverse transform to add azimuth and distance "info" to an ordered geodataframe and keep getting unexpected NaN results when I use a local UTM-based EPSG.

It works fine in WGS84 (my earlier question here), but projected systems are better than geographic for getting reliable azimuth and distance metrics, which is why I want it to work in other EPSGs.

Below are two examples for two different EPSG codes, and both return NaN for the azimuth and distance fields. Both examples use the same syntax/approach, the point is that the issue is repeatable in two different EPSG codes.

Running the following:
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

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

Example 1

myid = [1, 1, 1]
myorder = [1, 2, 3]
x = [550338.0319, 550428.0048, 550523.9951, 550589.9544]
y = [3795929.972, 3795798.055, 3795659.962, 3795528.029]
myepsg = 32611

df = pd.DataFrame(list(zip(myid, myorder, y, x)), columns =['myid', 'myorder', 'y', 'x']) 
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']))
gdf_pt = gdf_pt.set_crs(epsg=myepsg)

ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.ticklabel_format(style='plain', axis='both', useOffset=False)

geod  = CRS.from_epsg(myepsg).get_geod()
for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = geod.inv(gdf_pt.x[i], gdf_pt.y[i], gdf_pt.x[i-1], gdf_pt.y[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])

display(gdf_pt)

Example 2

myid = [1, 1, 1]
myorder = [1, 2, 3]
lat = [5174925.07851924, 5174890.26832387, 5174855.45812849]
long = [1521631.6994673, 1521667.11033893, 1521672.52121056]
# typo above, it says lat/long but it really is UTM-y & UTM-x
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)

geod  = CRS.from_epsg(myepsg).get_geod()

ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.ticklabel_format(style='plain', axis='both', useOffset=False)

for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = geod.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])

display(gdf_pt)

Answer 1

I think the issue is that you are using projected coordinates, not geographic coordinates (even though you call them lat/lon in one example). The pyproj.Geod object needs geographic coordinates. So, you need the geographic equivalent for the points:

gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']), crs=myepsg)
gdf_pt_geog = gdf_pt.to_crs(gdf_pt.crs.geodetic_crs)
geod = gdf_pt_geog.crs.get_geod()

Then, use the geographic data in the rest of the process (Note: Using the properties x and y of the geometry as that contains the geographic coordinates).

for i, r in gdf_pt_geog.iloc[1:].iterrows():
    myinfo = geod.inv(gdf_pt_geog.geometry.x[i], gdf_pt_geog.geometry.y[i], gdf_pt_geog.geometry.x[i-1], gdf_pt_geog.geometry.y[i-1])
    gdf_pt_geog.loc[i, 'az_fwd'] = myinfo[0]
    gdf_pt_geog.loc[i, 'az_back'] = myinfo[1]
    gdf_pt_geog.loc[i, 'dist'] = myinfo[2]
    gdf_pt_geog.loc[i, 'bearing'] = max(myinfo[1], myinfo[0])