pyproj unexpectedly returns nan for azimuth (bearing)
I am trying to use pyproj inverse to get the forward bearing, backward bearing, and distance between two points, as discussed in the first answer here. However, I get 'nan' for any tries. As far as I know, I'm using the right ellipsoid. A related curiosity is, how can I extract ellipsoid info from a geodataframe, in the format that is needed for the inv input, using pyproj CRS?
Thank you for any suggestions
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
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj
myid = [1, 1, 1, 1, 1]
myorder = [1, 2, 3, 4, 5]
lat = [36.42, 36.4, 36.32, 36.28,36.08]
long = [-118.11, -118.12, -118.07, -117.95, -117.95]
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=4326)
print(gdf_pt.crs)
display(gdf_pt)
ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
######
geodesic = pyproj.Geod(ellps='WGS84')
# is there a way to get the ellipsoid info from "gdf_pt.crs"
fwd_azimuth,back_azimuth,distance = geodesic.inv(gdf_pt.lat[0], gdf_pt.long[0], gdf_pt.lat[1], gdf_pt.long[1])
print(fwd_azimuth,back_azimuth,distance)
## it returns nan?
This question was answered on GIS Stack Exchange here and here. In case someone comes across this post, I'll summarize the answer:
The NaN error returned in the example in the question was due to mixing the input order: it should be fwd_azimuth,back_azimuth,distance = geodesic.inv(x1, y1, x2, y2) but was mistakenly entered as y1, x1, y2, x2.
The syntax to extract the geoid input is geod = CRS.from_epsg(myepsg).get_geod() where myepsg is an integer for the EPSG code (credit to the answer here).
Here is the working code and output:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj
from pyproj import CRS
### Build example Point GeoDataFrame ###
myid = [1, 1, 1, 1, 1]
myorder = [1, 2, 3, 4, 5]
lat = [36.42, 36.4, 36.32, 36.28,36.08]
long = [-118.11, -118.12, -118.07, -117.95, -117.95]
df = pd.DataFrame(list(zip(myid, myorder, lat, long)), columns =['myid', 'myorder', 'lat', 'long'])
df.sort_values(by=['myid', 'myorder'], inplace=True)
df.reset_index(drop=True, inplace=True)
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
gdf_pt = gdf_pt.set_crs(epsg=4326)
display(gdf_pt.style.hide_index())
ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
ctx.add_basemap(ax, crs='epsg:4326', source=ctx.providers.Esri.WorldShadedRelief)
# label the points just to double check
# zip joins x and y coordinates in pairs
for x,y,z in zip(gdf_pt.long, gdf_pt.lat, gdf_pt.myorder):
label = int(z)
# this method is called for each point
plt.annotate(label, # this is the text
(x,y), # this is the point to label
textcoords="offset points", # how to position the text
xytext=(0,10), # distance from text to points (x,y)
ha='center') # horizontal alignment can be left, right or center
### do analysis
myUTMepsg = 32611
geod = CRS.from_epsg(myUTMepsg).get_geod()
for i, r in gdf_pt.iloc[1:].iterrows():
#print(i, gdf_pt.long[i], gdf_pt.long[i-1])
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_meters'] = myinfo[2]
gdf_pt.loc[i, 'bearing_degrees'] = max(myinfo[1], myinfo[0])
display(gdf_pt.style.hide_index())
