0

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)

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

a11
  • 940
  • 10
  • 22
  • 1
    Side note: geopandas 7+ uses the pyproj CRS object: https://geopandas.readthedocs.io/en/latest/projections.html – snowman2 Jul 03 '20 at 16:23

1 Answers1

0

You can use CRS.get_geod:

Here is an example:

from pyproj import CRS

crs = CRS.from_epsg(epsg_code)
geod = crs.get_geod()
snowman2
  • 7,321
  • 12
  • 29
  • 54
  • 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 02 '20 at 20:37
  • 1
    geod is the pyproj.Geod object. No need to pass it in. – snowman2 Jul 03 '20 at 01:36