[GIS] Write a python function to split polyline shapefile with point shapefile

linepointpyqgispyshpshapefile

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)

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:

How to make the input shapefile directly readable by python function cut(line, point)?
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.

Best Answer

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

Related Solutions

[GIS] Checking for number of shapefiles in directory using Python

Check the glob library:

import os
import glob
path = "c:\\my\\dir"

flist = glob.glob(os.path.join(path,"shapefile*"))  <- Notice the *

>>> flist 
["c:\\my\\dir\shapefile.prj","c:\\my\\dir\shapefile.shp",...]

If you want just the fname without the path attached:

 os.path.basename(flist[0]) 
 >>> shapefile.prj

if you want just the file name without the extension:

os.path.basename(flist[0]).split(".")[0]
>>> shapefile

Also if the list is empty, then the file doesn't exist:

if flist not:
    print "File not found"
    break #or os.exit(1)

Edit: Also to check for the only filename in path :

import os
import glob
path = "c:\\my\\dir"

flist = glob.glob(os.path.join(path,"shapefile*"))  <- Notice the *
mlist = []
for file in flist:
    if os.path.basename(file).split(".")[0] not in mlist:
        mlist.append(os.path.basename(file).split(".")[0])
if len(mlist)>1: #(More than filename)
    break

To check if only one shapefile in the dir:

 def file_check(path):
      """Returns true if in path there is only one basename"""
      import os
      mlist = [] 
      for file in os.listdir(path):
          if os.path.basename(file).split(".")[0] not in mlist:
                mlist.append(os.path.basename(file).split(".")[0]
   if len(mlist)>1:
            return True
    else:
            return False

[GIS] Standalone Python script to split a polyline with a point layer

Here's a script that might do it for you. The native arcpy.Geometry class has a method called "cut" that will cut any feature using another polyline. Unfortunately, since you're using "points", we have to make "fake" lines out of these points. I essentially made a scratch polyline with the points [(Point.X+10, Point.Y+10), (Point.X-10, Point.Y-10)] - e.g. a diagonal line for every point to perform each cut. When you have more than one point intersecting a line, you have to re-cut the lines recursively smaller-and-smaller until they can't be cut up any further (without producing zero-length geometries).

A few caveats I should mention:

The script assumes there are no multi-part geometries in both feature classes. I always used Point.firstPoint for each point, so any additional multi-points thereafter will be left out. I haven't tested the effect of cutting a multi-part line, but I assume the descendants will also stay multi-part.
When you have a perfectly diagonal segment of a polyline, it won't get cut because our "fake cut line" is also diagonal. You may have to split these manually, if they exist.
The output won't have the attributes from the original polylines, but a Spatial Join can take care of that.
Lastly, of course, make sure both feature classes are in the same coordinate system.

You just have to make the skeleton tool in the toolbox that accepts 3 parameters:

Polyline Feature Class (Direction = Input, Type = Feature Class, Filter = Line)
Point Feature Class (Direction = Input, Type = Feature Class, Filter = Point)
Output Feature Class (Direction = Output, Type = Feature Class)

Here are some examples: enter image description here

And the code:

__author__ = "John K. Tran"
__contact__ = "jtran20@masonlive.gmu.edu"
__version__ = "4.0"
__created__ = "7/16/15"
__credits__ = "http://gis.stackexchange.com/questions/154708/standalone-python-script-to-split-a-polyline-with-a-point-layer"

"""
'Cuts', 'splits' or 'dices' a polyline feature class using a point feature class.
The output schema (i.e. list of fields) will be blank, but you can use a spatial join
to repopulate the attributes to the output feature class.
"""

import arcpy
import os
import sys

arcpy.env.overwriteOutput = True
arcpy.SetProgressor("default", "Firing up script...")

linefc = arcpy.GetParameterAsText(0) # Your polyline feature class
pointfc = arcpy.GetParameterAsText(1) # Your point feature class
output = arcpy.GetParameterAsText(2) # New output feature class.

spatialref1 = arcpy.Describe(pointfc).spatialReference
spatialref2 = arcpy.Describe(linefc).spatialReference
assert spatialref1.name == spatialref2.name, "Ensure both feature classes have the same projected coordinate system"

# If the Advanced license is available, just use the normal geoprocessing tools.
if arcpy.ProductInfo() in [u'ArcInfo', u'ArcServer']:
    arcpy.SetProgressorLabel("Splitting lines at points")
    arcpy.SplitLineAtPoint_management(linefc, pointfc, output, 1.0)
    arcpy.SetProgressorLabel("Deleting duplicate slices")
    outshapefieldname = arcpy.Describe(output).shapeFieldName
    arcpy.DeleteIdentical_management(output, [outshapefieldname])
    sys.exit(0)

# If not, proceed with script.
# Set up some preliminary variables to start.
arcpy.SetProgressorLabel("Gathering geometries")
points = [row[0] for row in arcpy.da.SearchCursor(pointfc, "SHAPE@")]
lines = [row[0] for row in arcpy.da.SearchCursor(linefc, "SHAPE@")]
cutlines = list()
iterations = 0

# Defining a function for the cut will let us reuse it again for subsequent cuts.
def CutLines():
    arcpy.SetProgressorLabel("Running cut: Iteration {0}".format(str(iterations)))
    for line in lines[:]:
        iscut = "Invalid"
        if line.length > 0.0: # Make sure it's not an empty geometry.
            iscut = "Not Cut"
            for point in points:
                if line.distanceTo(point) < 1.0: # Even "coincident" points can show up as spatially non-coincident in their floating-point XY values, so we set up a tolerance.
                    snappoint = line.snapToLine(point).firstPoint # To ensure coincidence, snap the point to the line before proceeding.
                    if not (snappoint.equals(line.lastPoint) and snappoint.equals(line.firstPoint)): # Make the sure the point isn't on a line endpoint, otherwise cutting will produce an empty geometry.
                        cutline1, cutline2 = line.cut(arcpy.Polyline(arcpy.Array([arcpy.Point(snappoint.X+10.0, snappoint.Y+10.0), arcpy.Point(snappoint.X-10.0, snappoint.Y-10.0)]), spatialref1)) # Cut the line.
                        if cutline1.length > 0.0 and cutline2.length > 0.0: # Make sure both descendents have non-zero geometry.
                            lines.insert(0, cutline1) # Feed the cut lines back into the "line" list to be recut.
                            lines.insert(0, cutline2) # The cut loop will only exit when all lines cannot be recut smaller and smaller (without producing zero-length geometries).
                            iscut = "Cut"
        if iscut == "Not Cut":
            cutlines.insert(0, line)
        lines.remove(line)

# Perform the pseudo-recurive Cut loop until the polyline FC can't be cut further.
while lines:
    CutLines()
    iterations += 1
    if iterations > 500: # Fail-safe to stop an infinite loop if something goes wrong (or if more than 500 points intersect a single line).
        break

# Create the output feature class.
arcpy.SetProgressorLabel("Creating output feature class")
arcpy.CreateFeatureclass_management(os.path.dirname(output), os.path.basename(output), "POLYLINE", None, None, None, spatialref2)

# Insert each cut line into the new feature class.
arcpy.SetProgressorLabel("Inserting cut lines")
with arcpy.da.InsertCursor(output, "SHAPE@") as insertcursor:
    for cutline in cutlines:
        insertcursor.insertRow((cutline,))

# Delete duplicates by comparing their geometry objects.
arcpy.SetProgressorLabel("Deleting duplicates")
crows = [row for row in arcpy.da.SearchCursor(output, ["OID@", "SHAPE@"])]
deleteOIDs = set() # Retain the OIDs for deletion, performed in the next cursor.
for drow in crows: # Compare each geometry to every other geometry in the output FC.
    dOID = drow[0]
    dgeom = drow[1]
    if dOID not in deleteOIDs:
        for srow in crows:
            sOID = srow[0]
            sgeom = srow[1]
            if sOID != dOID and sOID not in deleteOIDs: # Prevent redundant comparisons, only compare new combinations.
                dextent = dgeom.extent
                sextent = sgeom.extent
                if dextent.XMin == sextent.XMin and dextent.XMax == sextent.XMax and dextent.YMin == sextent.YMin and dextent.YMax == sextent.YMax: # See if their extents match first.
                    if dgeom.equals(sgeom): # Then see if they equal each other.
                        deleteOIDs.add(sOID) # If so, add the OID to the list of OIDs to delete.
with arcpy.da.UpdateCursor(output, "OID@") as deletecursor: # Delete the records here using your list of "DeleteOIDs".
    for delrow in deletecursor:
        delOID = delrow[0]
        if delOID in deleteOIDs:
            deletecursor.deleteRow()

# Clean up leftover variables.
namespace = dir()
for var in ["linecursor", "pointcursor", "insertcursor", "deletecursor"]:
    if var in namespace:
        exec "del " + var

arcpy.ResetProgressor()

If you get any errors, let me know.

Best Answer

Related Solutions

[GIS] Checking for number of shapefiles in directory using Python

[GIS] Standalone Python script to split a polyline with a point layer

Related Question