Using pyproj inverse transform with GeoPandas and GroupBy + SortBy constraints

Question

I am using pyproj inverse transform to add azimuth and distance "info" to a grouped and ordered geodataframe. I asked a similar question earlier here, but I didn't realize grouping was a constraint, and I cannot figure out how to incorporate GroupBy into that answer. (The code in that answer is very compact and a bit hard to understand for a novice).

In the example below, I want to perform pyproj inverse transform in 'myid' groupings. In other words, I don't care about the distance and azimuth between point "1-3" and point "2-1"; that is to say, the distances and azimuths for every "X-1" point should be nan. How can I achieve this?

%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
from pyproj import CRS

#### Build example GeoDataFrame for easily reproducible example ####
myUTMepsg = 32611
myid = [1, 1, 1, 2, 2]
myorder = [3, 2, 1, 2, 1]
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);
#### End

#### this is not necessary, but I like to add labels to the
#### points on the plot so that I can have a visual check
for x,y,z1,z2 in zip(gdf_pt.long, gdf_pt.lat, gdf_pt.myid, gdf_pt.myorder):
    label = str(int(z1)) + '-' + str(int(z2))
    plt.annotate(label, (x,y), textcoords = 'offset points', 
                 xytext = (12,-5), ha = 'center')
#### End    

#### Analysis & Problem Area
geod = CRS.from_epsg(myUTMepsg).get_geod()
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_meters'] = myinfo[2]
    gdf_pt.loc[i, 'bearing_degrees'] = max(myinfo[1], myinfo[0])
    
display(gdf_pt.style.hide_index())

I tried breaking out a working geodataframe inside the loop, but I ran into trouble trying to combine the working and final geodataframe; I kept getting duplicate entries (yellow below; red arrow is the one I want to keep), even when specifying keys, using Pandas merge, concatenate, and join.

#### Analysis & Problem Area
geod = CRS.from_epsg(myUTMepsg).get_geod()
for idval in gdf_pt.myid.unique():
    tempgdf = gdf_pt[gdf_pt.myid == idval].copy()
    tempgdf.sort_values(by=['myid', 'myorder'], inplace=True)
    tempgdf.reset_index(drop=True, inplace=True)
    for i, r in tempgdf.iloc[1:].iterrows():
        myinfo = geod.inv(tempgdf.long[i], tempgdf.lat[i], tempgdf.long[i-1], tempgdf.lat[i-1])
        tempgdf.loc[i, 'az_fwd'] = myinfo[0]
        tempgdf.loc[i, 'az_back'] = myinfo[1]
        tempgdf.loc[i, 'dist_meters'] = myinfo[2]
        tempgdf.loc[i, 'bearing_degrees'] = max(myinfo[1], myinfo[0])
    gdf_pt = pd.merge(gdf_pt, tempgdf, how='outer') ## I tried various key controls here without success
    del tempgdf

display(gdf_pt.style.hide_index())

Using Windows 10, conda 4.8.2, Python 3.8.3

Answer 1

I think a lot of the hassle comes from trying to loop through your dataframes. It's a difficult mental paradigm shift to make, but it's worth working through.

So that steps you need to solve your problem:

1. write a function to precompute your shifted geometries for each group
2. write a function to compute your metrics given a row with two points
3. pre-sort your dataframe
4. group by your myid columns + compute your shifted points
5. apply the function to compute your metrics

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj


def shift_points(gdf):
    g = pyproj.Geod(ellps='WGS84')
    
    return (
        gdf.assign(p1=gdf['geometry'], p2=gdf['geometry'].shift(1)) 
    )

def do_inv(row, g):
    if row['p1'] and row['p2']:
        vals = g.inv(row['p1'].x, row['p1'].y, row['p2'].x, row['p2'].y)
    else:
        vals = [None, None, None]
    return pd.Series(vals, index=['az_fwd', 'az_back', 'dist'])


myid = [1, 1, 1, 2, 2]
myorder = [3, 2, 1, 2, 1]
lat = [36.42, 36.4, 36.32, 36.28,36.08]
long = [-118.11, -118.12, -118.07, -117.95, -117.95]
g = pyproj.Geod(ellps='WGS84') 

gdf = (
    pd.DataFrame(list(zip(myid, myorder)), columns =['myid', 'myorder']) 
      .pipe(gpd.GeoDataFrame, geometry=gpd.points_from_xy(long, lat), crs='epsg:4326')
      .sort_values(by=['myid', 'myorder'])
      .reset_index(drop=True)
      .groupby('myid')
      .apply(shift_points)
      .pipe(lambda df: df.join(df.apply(do_inv, g=g, axis=1)))
      .drop(columns=['p1', 'p2'])
      .assign(bearing_deg=lambda df: df.loc[:, ['az_fwd', 'az_back']].max(axis=1)) 
)

And that gives me:

   myid  myorder                     geometry  az_fwd  az_back      dist  bearing_deg
0     1        1  POINT (-118.07000 36.32000)     NaN      NaN       NaN          NaN
1     1        2  POINT (-118.12000 36.40000)  153.17    -26.8   9947.07       153.17
2     1        3  POINT (-118.11000 36.42000) -157.99     22.0   2393.73        22.00
3     2        1  POINT (-117.95000 36.08000)     NaN      NaN       NaN          NaN
4     2        2  POINT (-117.95000 36.28000)  180.00      0.0  22192.47       180.00