4

I'm trying to split a shapefile of geometry type 'line' at a shapefile of geometry type 'point. These points are snapped to the lines and are the start points/end points of other line shapefiles.

Something like this : The points are snapped to the black lines.

splitAtPoints

As of now I'm working with Geopandas/Shapely/Fiona. I've searched for a native function in Shapely that performs the split; however, I've had no luck in finding a solution. This is also true when searching through Shapely's official docs. This is odd since the split is among the basic geoprocessing tools.

My goal is to give a function of the two shapefiles (lines & points) and split the lines at points (like ArcGIS' tool). The split lines should then be stored in a different shapefile. This seems easy to do but the lack of documentation is a bump in the road.

Any ideas on how I can do this in python using Geopandas/Shapely?

6

You have many solutions and I use here a simple example

1) the easiest way

from shapely.geometry import Point, LineString
line = LineString([(1,2),(2,4),(4,5)])
point = Point(2,4)

First, you must determine if the point is within the line (Determine if shapely point is within a linestring/multilinestring)

line.distance(point) < 1e-8
True
print LineString([line.coords[0],point.coords[:][0]])
LINESTRING (1 2, 2 4)
print LineString([point.coords[:][0], line.coords[-1]])
LINESTRING (2 4, 4 5)

2) with the shapely function split

from shapely.ops import split
result = split(line, point)
result.wkt
'GEOMETRYCOLLECTION (LINESTRING (1 2, 2 4), LINESTRING (2 4, 4 5))'

3) from Get the vertices on a LineString either side of a Point

from shapely.geometry import Point,LineString    
def split(line_string, point):
    coords = line_string.coords
    j = None    
    for i in range(len(coords) - 1):
        if LineString(coords[i:i + 2]).intersects(point):
           j = i
           break    
    assert j is not None    
    # Make sure to always include the point in the first group
    if Point(coords[j + 1:j + 2]).equals(point):
        return coords[:j + 2], coords[j + 1:]
    else:
        return coords[:j + 1], coords[j:]

line1,line2 = split(line,point)
line1 = LineString(line1)
line2 = LineString(line2)
print line1, line2
LINESTRING (1 2, 2 4) LINESTRING (2 4, 4 5)

4) from Shapely Split LineStrings at Intersections with other LineStrings

# First coords of line (start + end)
coords = [line.coords[0], line.coords[-1]] 
# Add the coords from the points
coords += point.coords
# Calculate the distance along the line for each point
dists = [line.project(Point(p)) for p in coords]
# sort the coordinates
coords = [p for (d, p) in sorted(zip(dists, coords))]
lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
for lin in lines:
   print lin
LINESTRING (1 2, 2 4)
LINESTRING (2 4, 4 5)

5) you can also examine and adapt

  • Thank you @gene for your detailed answer! At first I was looking for a native function in Shapely that does that. – GeoSal Jul 22 '16 at 7:53
  • The function split – gene Jul 22 '16 at 10:02
  • The split function in Shapely is not recognized apparently. I tried importing : from shapely.ops import split But I get an unresolved error. I install the split package and nothing happens still. Very odd. I already have the points and lines in shapefiles so all I need to do is the split. The point are already projected on the lines. – GeoSal Jul 22 '16 at 20:20
  • look at the github version of shapely ops.py line 329 – gene Jul 22 '16 at 20:58
  • Yes I have seen it. That's why I don't understand why it's not recognized once imported. – GeoSal Jul 22 '16 at 22:45

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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