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

Question

I have a polyline layer and a layer of points that are on top of the polylines. I do not have the ArcGis advanced license so I can't use the Split Line at Point tool. Can someone help me with a standalone Python script that will split my polyline by my point layer? (By the way, I know that there is a tool like this in ET Geowizards, but I don't have the license for that so I can't make that tool into a standalone script).

Answer 1

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:

1. Polyline Feature Class (Direction = Input, Type = Feature Class, Filter = Line)
2. Point Feature Class (Direction = Input, Type = Feature Class, Filter = Point)
3. Output Feature Class (Direction = Output, Type = Feature Class)

Here are some examples:

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.