4

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

enter image description here

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

1 Answer 1

5
+50

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))
3
  • Thank you, I was having trouble updating the geometry using the function from linked post. thanks for your help
    – a11
    Jul 17, 2020 at 0:39
  • What is the name of the kind of pythonic logic being used here where you refer to the function inside of itself? Can ask another question on S/O if it's more appropriate/deserving of a full question.
    – Heymans
    Jun 24 at 0:21
  • 1
    The term is recursive function, but I wouldn't say there is anything specifically pythonic about it, more of a general CS concept.
    – Nick
    Jun 24 at 1:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.