pyproj 意外 returns 方位角（方位角）的 nan

Question

我正在尝试使用 pyproj inverse to get the forward bearing, backward bearing, and distance between two points, as discussed in the first answer 。但是，我得到 'nan' 任何尝试。据我所知，我使用的是正确的椭圆体。一个相关的好奇心是，如何使用 pyproj CRS 以 inv 输入所需的格式从地理数据框中提取椭球体信息？

感谢您的任何建议

运行以下内容：
Windows 10
康达 4.8.2
Python 3.8.3
匀称 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?

Answer 1

此问题已在 GIS Stack Exchange here and here 上得到解答。如果有人遇到这个 post，我会总结一下答案：

题目中的例子返回NaN错误是因为输入顺序错乱：本应是fwd_azimuth,back_azimuth,distance = geodesic.inv(x1, y1, x2, y2)却被误输入为y1, x1, y2, x2.

提取大地水准面输入的语法是 geod = CRS.from_epsg(myepsg).get_geod()，其中 myepsg 是 EPSG 代码的整数（归功于答案 here）。

这是工作代码和输出：

%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])

显示(gdf_pt.style.hide_index())

pyproj 意外 returns 方位角（方位角）的 nan

pyproj unexpectedly returns nan for azimuth (bearing)

python

pyproj

bearing

azimuth