Update#1 I used the code provided by @GeospatialPython.com in the answer:

import os
import shapefile
from shapely.geometry import *


def cut(line_file, point_file, cut_file):
"""Cuts a line shapefile using a point shapefile."""
# Line shapefile we are going to cut
line_sf = shapefile.Reader(line_file)

# Point shapefile containing the cut points
point_sf = shapefile.Reader(point_file)

# Prefix for output shapefile file name
cut_fn = cut_file

# Create the shapefile writer for the output line shapefile
cut_sf = shapefile.Writer(shapefile.POLYLINE)

cut_sf.fields = list(line_sf.fields)

for sr in line_sf.shapeRecords():
    line = LineString(sr.shape.points)
    envelope = box(*line.bounds)
    pt = None
    for shape in point_sf.shapes():
        pt2 = Point(*shape.points[0])
        if envelope.contains(pt2):
            if not pt:
                pt = pt2
                continue
            if line.distance(pt2) < line.distance(pt):
                pt = pt2
    cut_line = None
    if not pt:
        cut_line = [line]
    else:
        coords = list(line.coords)
        distance = line.project(pt)
        if distance <= 0.0 or distance >= line.length:
            return [LineString(line)]
        for i, p in enumerate(coords):
            pd = line.project(Point(p))
            if pd == distance:
                return [LineString(coords[:i+1]), LineString(coords[i:])]
            if pd > distance:
                cp = line.interpolate(distance)
                return [
                    LineString(coords[:i] + [(cp.x, cp.y)]),
                    LineString([(cp.x, cp.y)] + coords[i:])]
    for part in cut_line:
        coords = list(part.coords)
        cut_sf.line([coords])
        cut_sf.record(*sr.record)

cut_sf.save("{}.shp".format(cut_fn))

# If the line shapefile has a cry file, copy it for the new file.
shp_name = os.path.splitext(line_sf.shp.name)[0]
line_prj = "{}.prj".format(shp_name)
if os.path.exists(line_prj):
    with open(line_prj, "r") as line_crs:
        crs = line_crs.read()
        with open("{}.prj".format(cut_fn), "w") as cut_crs:
            cut_crs.write(crs)

# Call the function with some shapefile names        
# cut("river", "points", "cut")
# cut('tc_500.shp', 'tc_stations.shp', 'cut.shp')


if __name__ == "__main__":


    cut('tc_line.shp', 'tc_stations.shp', 'tc_cut')

But returned nothing, is there some misses in this code?


I want to write a python function to split line shapefile with point shapefile (point may be near the line or falls exactly on the line), and the following is the schematic function I expect:

cut line = cut(line shapefile, point shapefile)

the function cut input line based on input point, and if the point is not one the line, it cuts the line at the nearest point (projection point) to the input point.

I now know that pyshp, shapely, etc., could handle shapefiles, and I used and modified code below, which was originally posted in the shapely manual to test the ability of pyshp and shapely.

import shapefile
from shapely.geometry import *
from pprint import pprint

def cut(line, distance):
    # Cuts a line in two at a distance from its starting point
    if distance <= 0.0 or distance >= line.length:
        return [LineString(line)]
    coords = list(line.coords)
    for i, p in enumerate(coords):
        pd = line.project(Point(p))
        if pd == distance:
            return [
                LineString(coords[:i+1]),
                LineString(coords[i:])]
        if pd > distance:
            cp = line.interpolate(distance)
            return [
                LineString(coords[:i] + [(cp.x, cp.y)]),
                LineString([(cp.x, cp.y)] + coords[i:])]



if __name__ == "__main__":

    sf = shapefile.Reader('tc_2000.shp')
    pprint([list(x.coords) for x in cut(sf, 0.5)])

and an error occurred:

AttributeError: Reader instance has no attribute 'length'

I have a polyline shapefile, which is the input shapefile 'tc_2000.shp' in the test code above, and another point shapefile loaded in QGIS, as the figure below, and obviously the points are not on the line. (The second figure is part of figure 1)

enter image description here enter image description here

I want to write a python function which inputs are point and polyline shapefile, and the output is a new polyline shapefile which is the original input polyline shapefile being split at position nearest to the point.

The schematic function may look like:

splited polyline = cut(line, point)

and here's the code I have written:

import shapefile
from shapely.geometry import *
from pprint import pprint

def cut(line, point):
    coords = list(line.coords)
    distance = line.project(point)
    if distance <= 0.0 or distance >= line.length:
        return [LineString(line)]
    for i, p in enumerate(coords):
        pd = line.project(Point(p))
        print pd
        if pd == distance:
            print 'pd=distance'
            return [LineString(coords[:i+1]), LineString(coords[i:])]
        if pd > distance:
            print 'pd>distance'
            cp = line.interpolate(distance)
            return [
                LineString(coords[:i] + [(cp.x, cp.y)]),
                LineString([(cp.x, cp.y)] + coords[i:])]




if __name__ == "__main__":

    line = LineString([(0, 0), (0, 1), (1, 1)])
    p = Point([(0.5, 0.5)])
    pprint([list(x.coords) for x in cut(line, p)])

and its result is correct:

0.0
1.0
pd>distance
[[(0.0, 0.0), (0.0, 0.5)], [(0.0, 0.5), (0.0, 1.0), (1.0, 1.0)]]

So now I have 2 major questions:

  1. How to make the input shapefile directly readable by python function cut(line, point)?
  2. Could this kind of function cut(line, point) I wrote handle point shapefile with many points and complicated line shapefile with tributaries (as the figure above)?

How to achieve the ultimate goal of getting a splited line?

I am a newbie trying to use python to handle shapefiles.

In your first example, the problem is the cut function expects a shapely LineString object and you passed it a shapefile object.

For the second part of your question, the following example takes a point shapefile with one or more points, and cuts a line shapefile using those points. The line shapefile can be as complex as you want:

import os
import shapefile
from shapely.geometry import *

def cut(line_file, point_file, cut_file):
    """Cuts a line shapefile using a point shapefile."""
    # Line shapefile we are going to cut
    line_sf = shapefile.Reader(line_file)

    # Point shapefile containing the cut points
    point_sf = shapefile.Reader(point_file)

    # Prefix for output shapefile file name
    cut_fn = cut_file

    # Create the shapefile writer for the output line shapefile
    cut_sf = shapefile.Writer(shapefile.POLYLINE)

    cut_sf.fields = list(line_sf.fields)

    for sr in line_sf.shapeRecords():
        line = LineString(sr.shape.points)
        envelope = box(*line.bounds)
        pt = None
        for shape in point_sf.shapes():
            pt2 = Point(*shape.points[0])
            if envelope.contains(pt2):
                if not pt:
                    pt = pt2
                    continue
                if line.distance(pt2) < line.distance(pt):
                    pt = pt2
        cut_line = None
        if not pt:
            cut_line = [line]
        else:
            coords = list(line.coords)
            distance = line.project(pt)
            if distance <= 0.0 or distance >= line.length:
                return [LineString(line)]
            for i, p in enumerate(coords):
                pd = line.project(Point(p))
                if pd == distance:
                    return [LineString(coords[:i+1]), LineString(coords[i:])]
                if pd > distance:
                    cp = line.interpolate(distance)
                    return [
                        LineString(coords[:i] + [(cp.x, cp.y)]),
                        LineString([(cp.x, cp.y)] + coords[i:])]
        for part in cut_line:
            coords = list(part.coords)
            cut_sf.line([coords])
            cut_sf.record(*sr.record)

    cut_sf.save("{}.shp".format(cut_fn))

    # If the line shapefile has a cry file, copy it for the new file.
    shp_name = os.path.splitext(line_sf.shp.name)[0]
    line_prj = "{}.prj".format(shp_name)
    if os.path.exists(line_prj):
        with open(line_prj, "r") as line_crs:
            crs = line_crs.read()
            with open("{}.prj".format(cut_fn), "w") as cut_crs:
                cut_crs.write(crs)

# Call the function with some shapefile names        
cut("river", "points", "cut")
  • I used your script and try this: cut('line.shp', 'point.shp', 'cut.shp'), but outputed nothing, is that I missunderstanding the third arg in your function? – Heinz Oct 17 '16 at 5:15

Your Answer

 
discard

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.