python - pyproj 意外返回 nan 方位角（方位角）

Question

我正在尝试使用pyproj inverse来获取前向轴承、后向轴承和两点之间的距离，如第一个答案中所述。但是，任何尝试我都会得到“nan”。据我所知，我使用的是正确的椭球体。一个相关的好奇心是，如何使用 pyproj CRS 从地理数据框中以 inv 输入所需的格式提取椭球信息？

感谢您的任何建议

运行以下命令：
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?

标签： pythonpyprojbearingazimuth

Answer 1

此处和此处在 GIS Stack Exchange 上回答了这个问题。万一有人看到这篇文章，我将总结答案：

问题示例中返回的 NaN 错误是由于混合了输入顺序：它应该是fwd_azimuth,back_azimuth,distance = geodesic.inv(x1, y1, x2, y2)但被错误地输入为y1, x1, y2, x2.

提取大地水准面输入的语法是EPSG 代码的整数geod = CRS.from_epsg(myepsg).get_geod()where （感谢此处的答案）。myepsg

这是工作代码和输出：

%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()）