Using shapely interpolate to evenly re-sample points on a linestring geodataframe at specified interval

Question

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:

1. 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)
   
   
2. 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))

Answer 1

Why not use the function defined in Mike T's answer from the post you linked to?

def redistribute_vertices(geom, distance):
    if geom.geom_type == 'LineString':
        num_vert = int(round(geom.length / distance))
        if num_vert == 0:
            num_vert = 1
        return LineString(
            [geom.interpolate(float(n) / num_vert, normalized=True)
             for n in range(num_vert + 1)])
    elif geom.geom_type == 'MultiLineString':
        parts = [redistribute_vertices(part, distance)
                 for part in geom]
        return type(geom)([p for p in parts if not p.is_empty])
    else:
        raise ValueError('unhandled geometry %s', (geom.geom_type,))

If you want to update the geometry of each line you can do ...

gdf_line_interpolate = gdf_line.copy()
gdf_line_interpolate.geometry = gdf_line.geometry.apply(redistribute_vertices,distance=15)
gdf_line_interpolate['nverts'] = gdf_line_interpolate.geometry.apply(lambda x: len(x.coords))

or if you just wanted to add the attribute column you could do

gdf_line_interpolate = gdf_line.copy()
gdf_line_interpolate['nverts'] = gdf_line.geometry.apply(redistribute_vertices,distance=15).apply(lambda x: len(x.coords))