1

I am using pyproj inverse transform to add azimuth and distance "info" to an ordered geodataframe.

I am trying to iterate over all the rows in the geodataframe to add the "info" as a tuple to a new column in the geodataframe, but there is something amiss with my lambda function.

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

Example geodataframe:

%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

myid = [1, 1, 1, 1, 1]
myorder = [1, 2, 3, 4, 5]
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']) 
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
gdf_pt = gdf_pt.set_crs(epsg=4326)
print(gdf_pt.crs)
display(gdf_pt)
ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);

enter image description here

Expected output: I am iterating over i row number, starting at the second row, and then using that row & the previous to calc the azimuths & distances. The results are as expected.

g = pyproj.Geod(ellps='WGS84') 
for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = g.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]

display(gdf_pt)

enter image description here

Wrong output, need help here!
I am trying to iterate over i row number, starting at the second row, and then using that row & the previous to calc the azimuths & distances. However, as seen in the table display, it is not passing through i as I expected, when compared to the expected results above. It is just reading the last effort and/or over-writing

g = pyproj.Geod(ellps='WGS84') 
for i, r in gdf_pt.iloc[1:].iterrows():
    gdf_pt['mytuple'] = gdf_pt.apply(lambda x: 
                                     g.inv(gdf_pt.long[i], gdf_pt.lat[i], 
                                           gdf_pt.long[i-1], gdf_pt.lat[i-1]), axis=1)
display(gdf_pt)

enter image description here

I also tried the following, which threw the error shown:

g = pyproj.Geod(ellps='WGS84') 
for i, r in gdf_pt.iloc[1:].iterrows():
    gdf_pt.loc[i, 'mytuple'] = gdf_pt.apply(lambda x: 
                                     g.inv(gdf_pt.long[i], gdf_pt.lat[i], 
                                           gdf_pt.long[i-1], gdf_pt.lat[i-1]), axis=1)
display(gdf_pt)

ValueError: Incompatible indexer with Series
  • you over write the column with each iteration. use gdf.assign – Paul H Jul 1 at 0:37
2

As a general rule, don't iterate over the rows of a dataframe. As a more specific rule, you don't want want to try to assign values in a loop. Since each of your iterations is creating three values, I think your best bet is to make two dataframes and join them together. Note that you really don't need geopandas for this, but since you're using it, I took advantage of it:

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

myid = [1, 1, 1, 1, 1]
myorder = [1, 2, 3, 4, 5]
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')
      .join(pd.DataFrame([
          g.inv(p1.x, p1.y, p2.x, p2.y)
          if p1 and p2 else (None, None, None)
          for p1, p2 in zip(gdf.geometry, gdf.geometry.shift(1))
          ], columns=['ax_fwd', 'az_back', 'dist']
      ))
)

And that gives me:

   myid  myorder                     geometry     ax_fwd     az_back          dist
0     1        1  POINT (-118.11000 36.42000)        NaN         NaN           NaN
1     1        2  POINT (-118.12000 36.40000)  22.003207 -157.990857   2393.731315
2     1        3  POINT (-118.07000 36.32000) -26.802568  153.167789   9947.072638
3     1        4  POINT (-117.95000 36.28000) -67.582667  112.346292  11656.645203
4     1        5  POINT (-117.95000 36.08000)   0.000000  180.000000  22192.468253
| improve this answer | |
  • There is a follow-up question here on adding a GroupBy constraint: gis.stackexchange.com/questions/367098/… – AlexS1 Jul 7 at 19:00
  • This 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

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.