1

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)

enter image description here

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)

enter image description here

  • You appear to be passing your projected coordinates directly as latitudes and longitudes. Values such as 5174925.07851924, 1521631.6994673are 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 at 18:16
0

I think the issue is that you are using projected coordinates, not geographic coordinates (even though you call them lat/lon in one example). The pyproj.Geod object needs geographic coordinates. So, you need the geographic equivalent for the points:

gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']), crs=myepsg)
gdf_pt_geog = gdf_pt.to_crs(gdf_pt.crs.geodetic_crs)
geod = gdf_pt_geog.crs.get_geod()

Then, use the geographic data in the rest of the process (Note: Using the properties x and y of the geometry as that contains the geographic coordinates).

for i, r in gdf_pt_geog.iloc[1:].iterrows():
    myinfo = geod.inv(gdf_pt_geog.geometry.x[i], gdf_pt_geog.geometry.y[i], gdf_pt_geog.geometry.x[i-1], gdf_pt_geog.geometry.y[i-1])
    gdf_pt_geog.loc[i, 'az_fwd'] = myinfo[0]
    gdf_pt_geog.loc[i, 'az_back'] = myinfo[1]
    gdf_pt_geog.loc[i, 'dist'] = myinfo[2]
    gdf_pt_geog.loc[i, 'bearing'] = max(myinfo[1], myinfo[0])
| improve this answer | |
  • so if I use lat/long in the geodataframe and gdf_pt = gdf_pt.set_crs(epsg=4326), then will geod = CRS.from_epsg(myUTMepsg).get_geod() cause the geod.in() function to compute the azimuth & distance metrics in UTM? In other words, there is no issue with bad distances/bearings that you would get from using the geographic coordinates? – AlexS1 Jul 5 at 1:18
  • 1
    The Geod doesn't work with projected coordinate systems like UTM as it is based on the ellipsoid alone. It uses geodesic distances: en.m.wikipedia.org/wiki/Geodesic. So, epsg:4326 would work. – snowman2 Jul 5 at 2:03
  • Your example is not recommended as the ellipsoid might be different with 4326 and the UTM projection. The example I provided above would be the safest method. – snowman2 Jul 5 at 2:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.