I am using pyproj inverse transform to add azimuth and distance "info" to an ordered geodataframe and keep getting unexpected NaN
results when I use a local UTM-based EPSG.
It works fine in WGS84 (my earlier question here), but projected systems are better than geographic for getting reliable azimuth and distance metrics, which is why I want it to work in other EPSGs.
Below are two examples for two different EPSG codes, and both return NaN
for the azimuth and distance fields. Both examples use the same syntax/approach, the point is that the issue is repeatable in two different EPSG codes.
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
from matplotlib.ticker import ScalarFormatter
import contextily as ctx
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj
from pyproj import CRS
Example 1
myid = [1, 1, 1]
myorder = [1, 2, 3]
x = [550338.0319, 550428.0048, 550523.9951, 550589.9544]
y = [3795929.972, 3795798.055, 3795659.962, 3795528.029]
myepsg = 32611
df = pd.DataFrame(list(zip(myid, myorder, y, x)), columns =['myid', 'myorder', 'y', 'x'])
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']))
gdf_pt = gdf_pt.set_crs(epsg=myepsg)
ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.ticklabel_format(style='plain', axis='both', useOffset=False)
geod = CRS.from_epsg(myepsg).get_geod()
for i, r in gdf_pt.iloc[1:].iterrows():
myinfo = geod.inv(gdf_pt.x[i], gdf_pt.y[i], gdf_pt.x[i-1], gdf_pt.y[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)
Example 2
myid = [1, 1, 1]
myorder = [1, 2, 3]
lat = [5174925.07851924, 5174890.26832387, 5174855.45812849]
long = [1521631.6994673, 1521667.11033893, 1521672.52121056]
# typo above, it says lat/long but it really is UTM-y & UTM-x
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)
geod = CRS.from_epsg(myepsg).get_geod()
ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.ticklabel_format(style='plain', axis='both', useOffset=False)
for i, r in gdf_pt.iloc[1:].iterrows():
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'] = myinfo[2]
gdf_pt.loc[i, 'bearing'] = max(myinfo[1], myinfo[0])
display(gdf_pt)
5174925.07851924, 1521631.6994673
are invalid geographic coordinates, therefore the algorithm fails to calculate geodesics. Make sure to use the transform function to convert from your projected crs to a geographic crs like WGS84 (epsg 4326) before calculating geodesics. – FSimardGIS Jul 4 '20 at 18:16