I have a GeoPandas LineString GeoDataFrame. I am trying to evenly interpolate the spacing of the nodes that make up each line, with a couple constraints: keep endpoints unchanged, and specify maximum spacing interval. The original nodes should be removed if they are not needed (i.e., a node on a straight line versus a node at an angle on the line).
Desired output format is GeoPandas LineString GeoDataFrame with the following attributes: lineid, geometry.
I am using GeoPandas and Shapely libraries and would prefer to find a solution using Shapely. There is an arcpy tool that can do something similar (GeneratePointsAlongLines()
), but I am looking for an open source solution.
I am trying to implement the solution here because it only requires Shapely (as opposed to solution here that requires OGR), but I am stuck on a couple parts:
- my test script adds a column for every linestring. In the screenshot below, the first line should have nverts =22, and the second line should be 16 (in other words, the 0-column should read [22,16] and there should not be a 1-column)
- when the
return LineStr...
part of my function is not commented out, I get an error:TypeError: 'Series' object cannot be interpreted as an integer
I think this is maybe because it is not handling multiple lineids correctly.
Where am I going wrong in my script, or is there a better open-source, out-of-the-box solution for my line interpolation?
Using Windows 10, conda 4.8.2, Python 3.8.3
Libraries
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
from shapely import wkt
import pyproj
from pyproj import CRS
Setup repeatable example
##### BUILD EXAMPLE LINESTRING GEODATAFRAME #####
myid = [1, 1, 1, 2, 2, 2]
myorder = [1, 2, 3, 1, 2, 3]
x = [550338, 550428, 550523, 550529, 550689, 550703]
y = [3795929, 3795798, 3795659, 3795419, 3795290, 3795257]
myepsg = 32611
df = pd.DataFrame(
list(zip(myid, myorder, x, y)),
columns =['myid', 'myorder', 'x', 'y']
)
gdf_pt = gpd.GeoDataFrame(
df, geometry=gpd.points_from_xy(df['x'], df['y'])
)
gdf_pt = gdf_pt.set_crs(epsg=myepsg)
gdf_line = (
gdf_pt.sort_values(by='myorder')
.groupby('myid')['geometry']
.apply(lambda x: LineString(x.tolist()))
)
gdf_line = gpd.GeoDataFrame(gdf_line, geometry='geometry')
gdf_line.reset_index(drop=False, inplace=True)
## delete stuff used to build because I don't
## have access to it in a real application
del myid, myorder, x, y, myepsg
Plotting
display(gdf_pt.style.hide_index())
display(gdf_line.style.hide_index())
ax = gdf_line.plot();
gdf_pt.plot(ax=ax)
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)
## this is not necessary, but I like to add labels
for x,y,z1,z2 in zip(gdf_pt.x, gdf_pt.y, 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')
Analysis & Problem Area
def mytest(row, gdf, maxlen):
nverts = round(gdf.length / maxlen).astype(int)
return nverts
#return LineString(
#[gdf.interpolate(float(n) / nverts, normalized=True)
#for n in range(nverts + 1)])
lineid = gdf_line.myid
geometry = gdf_line.geometry
testout = (
pd.DataFrame(list(zip(lineid, geometry)), columns =['lineid', 'geometry'])
.pipe(gpd.GeoDataFrame, geometry=geometry)
.pipe(lambda df: df.join(df.apply(mytest, gdf = gdf_line, maxlen = 15, axis=1)))
)
display(testout)
print(type(testout))