This is my script:

def cut_line_at_points(line, points):
    line is gdf format
    points are series of points in gdf format
    # union all geometries - convert the gdf to shpely (multi)line and (multi)point
    line = line.geometry.unary_union
    points = points.geometry.unary_union

    # check if line in multi-part and merge all parts to one linestring
    if line.geom_type != 'LineString':
        # Put the sub-line coordinates into a list of sublists
        outcoords = [list(i.coords) for i in line]
        # Flatten the list of sublists and use it to make a new line
        line = LineString([i for sublist in outcoords for i in sublist])

    # First coords of line
    coords = list(line.coords)

    # Keep list coords where to cut (cut = 1)
    cuts = [0] * len(coords)
    cuts[0] = 1
    cuts[-1] = 1

    # Add the coords from the points
    coords += [list(p.coords)[0] for p in points]
    cuts += [1] * len(points)

    # Calculate the distance along the line for each point
    dists = [line.project(Point(p)) for p in coords]

    # create new df and sort again if geom_type was LineString

    # sort the coords/cuts based on the distances
    coords = [p for (d, p) in sorted(zip(dists, coords))]
    cuts = [p for (d, p) in sorted(zip(dists, cuts))]

    # generate the Lines
    lines = []
    for i in range(len(coords)-1):
        if cuts[i] == 1:
            # find next element in cuts == 1 starting from index i + 1
                j = cuts.index(1, i + 1)
                cuts[-1] = 1
                j = cuts.index(1, i + 1)

    # transform Geometry Collection to GeoDataFrame
    segments = [feature for feature in lines]
    # set the geometry and other columns
    gdf_segments = gpd.GeoDataFrame(list(range(len(segments))), geometry=segments)
    gdf_segments.columns = ['index', 'geometry']

    return gdf_segments

everything is working really well when it's a shapely.LineString

enter image description here

but when I'm trying to cut the Multi-linestring after the flattening it all goes south:

enter image description here

I found out that when Calculating the distance along the line it scrambles the points distances, in this line:

dists = [line.project(Point(p)) for p in coords]

but I don't have any idea how to fix it. I use Python 3.7 and Shapely.


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.

Browse other questions tagged or ask your own question.