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
print(myellipsoid)
print(gdf_pt.crs)
display(gdf_pt)
ax = gdf_pt.plot();
ax.set_aspect('equal')
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])
display(gdf_pt)
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
