Splitting line GeoDataFrame basing on given lengths

Question

I am trying to split line GeoDataFrame into sections. I have lengths on which splits should be made stored in column. My idea was to generate points along the line using shapely.ops.interpolate and then split lines with shapely.ops.split and collect all the results in new GeoDataFrame. Unfortunatelly, shapely.ops.split returns given line as if generated points would not intersect with it. Here is my code:

lines = []
for split_distances, line in zip(suln1234['split_distances'], suln1234['geometry']):
    points = []
    if split_distances:
        for split_distance in split_distances:
            points.append(line.interpolate(split_distance))
        splitted = shapely.ops.split(line, shapely.ops.snap(shapely.geometry.MultiPoint(points), line, 1))
        for i in splitted:
            lines.append(i)
    else:
        lines.append(line)
l_split = gpd.GeoDataFrame({'geometry': lines}, crs=suln1234.crs)

As you can see, I also tried to use shapely.ops.snap, but even with crazy big tolerance, it only works for part of lines. What is the reason for such behaviour? I would imagine that if points are derived directly from line, they should allign nicely with its geometry. How to fix it?

Answer 1

Unfortunately, due to floating point issues, a point being on the linestring (due to nearest_points or interpolate) does not always result in the point intersecting with the linestring. There are several issues in the shapely repository discussing this exact pain point. The answer is just better documentation.

In order to guarantee that a point will split a line, you need to snap the Linestring to the Point (or MultiPoint), not the other way around.

from shapely.ops import split, snap
def split_line_by_point(line, point, tolerance: float=1.0e-12):
    return split(snap(line, point, tolerance), point)

The reason it works this way and not the other is that snapping a LineString to a Point inserts the point as a vertex of the linestring itself and the two geometries then reliably intersect. Snapping a Point to a LineString only places the point somewhere on the line, which is still subject to precision issues when interpolating between two vertices to test if the point is on the line. Having second thoughts. Snapping a Point to a LineString should snap the point to a vertex of the line, not just anywhere as I was thinking. But if a vertex is not nearby (in your case, within 1) even if the point is on the line it will not be placed on a vertex. Maybe that is what is going on in your case? Could the vertices of your line be more than 1 apart? Either way it seems preferable to snap the LineString to the Point not the Point to the LineString.

Answer 2

This is not the common way to iterate over rows in a (Geo)DataFrame, see How to iterate over rows in a DataFrame in Pandas, Different ways to iterate over rows in Pandas Dataframe and many more.

As you and Jeremiah England I use the snap option to split the lines (Shapely issue: Floating point issue with splitter points )

import geopandas as gpd
lines = gpd.read_file("lines.shp")
print(lines)
    dist           geometry
0    10  LINESTRING (1.00000 2.00000, 8.00000 7.00000, ...
1    15  LINESTRING (2.09821 15.04000, 0.12726 10.25341...

With Pandas itertuples() for example

from shapely.ops import snap, split
for i, row in enumerate(lines.itertuples()):
    splited = split(snap(row.geometry,row.geometry.interpolate(row.dist),0.01),row.geometry.interpolate(row.dist))
    for geom in splited:
        print(i, geom.wkt)
0 LINESTRING (1 2, 8 7, 6.749881714669379 6.374940857334689)
0 LINESTRING (6.749881714669379 6.374940857334689, 4 5, 2 4, 4 7, 8 5, 9 18, 1 2, 12 7, 4 5, 6 5, 4 9)
1 LINESTRING (2.098212290502795 15.04, 0.1272625698324044 10.25340782122905, 9.682839377571781 7.974770274768124)
1 LINESTRING (9.682839377571781 7.974770274768124, 9.888156424581005 7.925810055865924, 9.625363128491619 3.420782122905029, 13.43586592178771 3.664804469273745)
segs = []
distances = []
for row in lines.itertuples():
    splited = split(snap(row.geometry,row.geometry.interpolate(row.dist),0.01),row.geometry.interpolate(row.dist))
    for geom in splited:
      segs.append(geom)
      distance.append(row.dist)
l_split = gpd.GeoDataFrame({'dist':distance,'geometry': segs})
print(l_split)

    dist                  geometry
0    10  LINESTRING (1.00000 2.00000, 8.00000 7.00000, ...
1    10  LINESTRING (6.74988 6.37494, 4.00000 5.00000, ...
2    15  LINESTRING (2.09821 15.04000, 0.12726 10.25341...
3    15  LINESTRING (9.68284 7.97477, 9.88816 7.92581, ...

or directly

cols = ['dist','geometry']
l_split = gpd.GeoDataFrame(columns=cols)
for row in tt.itertuples():
   splited = split(snap(row.geometry,row.geometry.interpolate(row.dist),0.01),row.geometry.interpolate(row.dist))
for geom in splited:
    l_split = l_split.append({'dist': row.dist,'geometry':j },ignore_index=True)
print(l_split)     
    dist                  geometry
0    10  LINESTRING (1.00000 2.00000, 8.00000 7.00000, ...
1    10  LINESTRING (6.74988 6.37494, 4.00000 5.00000, ...
2    15  LINESTRING (2.09821 15.04000, 0.12726 10.25341...
3    15  LINESTRING (9.68284 7.97477, 9.88816 7.92581, ...

New

for dist, geom in zip(tt.dist, tt.geometry):
    splited = split(snap(geom,geom.interpolate(dist),0.01),geom.interpolate(dist))
    for i in splited:
        print(i.wkt)
LINESTRING (1 2, 8 7, 6.749881714669379 6.374940857334689)
LINESTRING (6.749881714669379 6.374940857334689, 4 5, 2 4, 4 7, 8 5, 9 18, 1 2, 12 7, 4 5, 6 5, 4 9)
LINESTRING (2.098212290502795 15.04, 0.1272625698324044 10.25340782122905, 9.682839377571781 7.974770274768124)
LINESTRING (9.682839377571781 7.974770274768124, 9.888156424581005 7.925810055865924, 9.625363128491619 3.420782122905029, 13.43586592178771 3.664804469273745)

And with list comprehension;

segs = [split(snap(geom,geom.interpolate(dist),0.01),geom.interpolate(dist)) for dist, geom in zip(tt.dist, tt.geometry)]