pyproj неожиданно возвращает nan для азимута (азимута)

Я пытаюсь использовать pyproj inverse, чтобы получить направление вперед, назад направление и расстояние между двумя точками, как обсуждалось в первом ответе здесь. Тем не менее, я получаю «нан» за любые попытки. Насколько я знаю, я использую правильный эллипсоид. Связанное с этим любопытство заключается в том, как я могу извлечь информацию об эллипсоиде из кадра геоданных в формате, необходимом для ввода inv, с помощью pyproj CRS?

Спасибо за любые предложения

Выполняется следующее:
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 -кузница

%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?

введите здесь описание изображения


pyproj azimuth bearing
person a11    schedule 30.06.2020    source источник
comment
NaN был глупой ошибкой в ​​упорядочении широта/долгота и исправлен с помощью этого: fwd_azimuth,back_azimuth,distance = geodesic.inv(gdf_pt.long[0], gdf_pt.lat[0], gdf_pt.long[1], gdf_pt.lat[1])   -  person a11    schedule 30.06.2020
comment
Однако я до сих пор не знаю, как извлечь ввод для pyproj.Geod(ellps=?) из информации .crs; жду советов сообщества по этому поводу   -  person a11    schedule 30.06.2020

Ответы (1)


arrow_upward
0
arrow_downward

Ответ на этот вопрос был дан на GIS Stack Exchange здесь и здесь. Если кто-то наткнется на этот пост, я резюмирую ответ:

Ошибка 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 (кредит на ответ здесь).

Вот рабочий код и вывод:

%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()) введите здесь описание изображения

person a11    schedule 07.07.2020