# Using pyproj inverse transform with GeoPandas and GroupBy + SortBy constraints

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
gdf_pt.loc[i, 'az_back'] = myinfo
gdf_pt.loc[i, 'dist_meters'] = myinfo
gdf_pt.loc[i, 'bearing_degrees'] = max(myinfo, myinfo)

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
tempgdf.loc[i, 'az_back'] = myinfo
tempgdf.loc[i, 'dist_meters'] = myinfo
tempgdf.loc[i, 'bearing_degrees'] = max(myinfo, myinfo)
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

• Remember my previous answer that didn't use any loops? Start with that, then do `df.groupby('myid').apply(previous_answer)` Jul 7 '20 at 19:16
• IOW, don't loop through dataframes. that road only leads to hassle Jul 7 '20 at 19:17
• @PaulH Your previous answer is hard to parse for a novice; zips and lists and pipes and joins and ifs and fors, all in one line of code without any explanations.
– a11
Jul 7 '20 at 19:29
• I think you'll find that feeling to be true for any code that wasn't written by the reader. I find your code hard to parse. What's important is that you find a good interactive environment (like Jupyter Notebooks) to run and experiment with the code. Jul 7 '20 at 19:54
• @PaulH Thank you for the comments, answer, and explanations. I'll plug it into my notebook this afternoon and give it a try
– a11
Jul 7 '20 at 20:05

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
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
``````
• Thank you, this works perfectly, just a minor change for my applications: `gdf = (gdf_pt.groupby('myid').apply(shift_points).pipe(lambda df: df.join(df.apply(do_inv, g=geod, axis=1))).drop(columns=['p1', 'p2']).assign(bearing_deg=lambda df: df.loc[:, ['az_fwd', 'az_back']].max(axis=1)))`. This is because in practice, I have a geodataframe as my starting point. I only include the building of one in the example so that the questions is minimal/complete/stand-alone example per the guidelines (stackoverflow.com/help/how-to-ask). Thanks
– a11
Jul 8 '20 at 0:03