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

  • Maybe this question gis.stackexchange.com/questions/61474/… helps you. – jgrocha Jul 16 '15 at 19:30
  • That uses Fiona/Shapely, which are two things I don't have downloaded. I was wondering if anyone could help me out with a Python standalone script. – Tiffany Jul 16 '15 at 20:48
  • You're going to want to use some sort of API/Library. The format of shapefiles is published but I, wouldn't like to tackle binary I/O of a shapefile in python. I'd give it a go in C# or C++ but why bark if you've got a dog there are a few good APIs' like OGR and Shapely that will do all the hard work, you don't need to stick to arcpy just because you've got ArcGis! It's only one (very good) tool in a complete GIS toolkit. If you're sure you want to do this in arcpy please post what you've got so far so we can see where you're stuck. – Michael Stimson Jul 16 '15 at 23:39
  • There is a thread on Cutting Polygons by Polylines that may be helpful, at least structure-wise, to writing a script the Cuts Polylines by Points. – John Jul 17 '15 at 4:16
up vote 6 down vote accepted

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: enter image description here 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.

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.