1

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? enter image description here

%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())

enter image description here

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())

enter image description here

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) – Paul H Jul 7 at 19:16
  • IOW, don't loop through dataframes. that road only leads to hassle – Paul H Jul 7 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. – AlexS1 Jul 7 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. – Paul H Jul 7 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 – AlexS1 Jul 7 at 20:05
2

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
| improve this answer | |
  • 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 – AlexS1 Jul 8 at 0:03

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.