Automated/geoprocesing tool to slice, clip, or cut polygons using polylines using ArcGIS Desktop? [duplicate]

Question

Looking for a simpler method here.

I am trying to split up / cut / slice existing polygons by using existing poly-lines. An example would be splitting a water body or land parcel at the point where a bridge/road crosses it. But the road network may not necessarily result in a closed loop.

Since the poly-lines are not necessarily interconnected or continuous, creating a polygon from them is not an option (which eliminates using the split tool). Also, I have tried using a topology with the geometry, but it keeps failing, probably due to the large/complex geometry.

Current Workflow: I have accomplished this with the Feature to Polygon tool, combining the lines and polygons together, but it creates additional unwanted polygons anywhere a loop closes. I created a mask of the original polygons, and used this with the overlay->erase tool to remove the unwanted polygons. This still leaves some unwanted geometry (mostly slivers), but is somewhat workable.

This seems like an extremely convoluted and round-a-bout way to accomplish (what seems like it should be) a very simple task.

Other than manual edits, or using a topology, is there is tool that could accomplish this in a single step?

Using: ArcMap \ ArcInfo Desktop 10 SP5

Edit 1: In my case, it is not actually roads as mentioned above. I have water polygons for a coastal area, and need to split the polygons where rock dams or levees have been put across the waterways. Which are typically not interconnected.

The water polygons have been simplified and repaired, to the point where I would no longer call the data "dirty" just complex and large. But I've got the workaround mentioned above working for this case.

I'm looking more "in general" for a tool that can simply divide polygons using polylines.

Edit 2:

Mapperz:Thanks for the Model Builder suggestion. I am going to use that as a stop gap solution, for now.

Jakub: Thanks for the suggestion. I'm not opposed to a programmatic solution or developing a custom tool, although I've never made one before. I have programming experience, but not in conjunction with Arc. I'd prefer something that directly cuts the geometry though, instead of following the roundabout logistics above. In theory that should reduce resultant gaps, since it wouldn't be subject to multiple iterations of cracking/clustering. Although, I'm not sure that would be as easy or even possible.

Edit 3: I'm looking for something to function like on the left in the image.

Answer 1

I thought there must be a way to do this, so I created my what I believe to be a pretty good solution. I have posted it on the ArcGIS Resources site and in the Community->Technical->Analysis & Geoprocessing->Analysis->Gallery.

The tool is called Split Polygons With Lines and requires an ArcInfo license because of some of the tools used within the model. Essentially what I did was create the minimum bounding box for the polygons and extend the lines to them. So using some ModelBuilder voodoo, I was able to turn the linework into polygons, which then I used Identity to split the original polys.

Please test it out and see if it works for you. In my (limited) tests it preserved attributes of the original polygons, and split only the existing polygons.

Answer 2

I believe "Feature to Polygon" does exactly what you need. You can input a combination of Polygon and Polyline feature classes. The output Polygons are split at each Polyline. An ArcInfo license is required, which you have. Tested in 10.0.

Make sure you have a field that is populated, for all features, before running "Feature to Polygon. The new Polygon Fills will have all blank fields. Another method is to spatially select the "cut" polygons with the original Polygons. The "cut" filled polys will not "have their center" in the original polygons.

Answer 3

This is a long shot, but do you have a 'Production Mapping', 'Defense Mapping', 'Aviation', or 'Maritime Charting' extension for ArcGIS? It seems the "Dice Polygons" tool will accomplish what you need. I've never used the tool personally, but from the description, it sounds like a winner.

Search Link

---

EDIT:

Alternatively, I've written a dice algorithm for you based on the discussion in "Optimizing ArcPy code to cut polygon?" to cut/split/dice a polygon feature class using a polyline feature class. Credits to those on the thread.

You just need to make a skeleton tool that accepts as parameters:

1. Input Polygon Feature Class (Type: Feature Class, Filter: Polygon)
2. Cutting Polyline Feature Class (Type: Feature Class, Filter: Polyline)
3. Output Polygon Feature Class (Type: Feature Class, Direction: Output

There are two assumptions I noticed while running the tool:

• Your polyline must cross the polygon entirely, no stopping halfway or passing off to a second line inside the polygon to finish off the traverse.
• If you have multi-part polygons, the resulting slices will also sometimes be multi-part, since they inherit the partitioning of their parents. Self-intersecting polygons are fine and have no issue.

Code below for saving as .py file:

__author__ = "John K. Tran, Tristan Forward"
__contact__ = "[email protected], https://gis.stackexchange.com/users/6996/tristan-forward"
__version__ = "2.0"
__created__ = "7/3/15"
__credits__ = "https://gis.stackexchange.com/questions/124198/optimizing-arcpy-code-to-cut-polygon"

"""
Cut polygons by polylines, splitting each polygon into slices.
:param to_cut: The polygon to cut.
:param cutter: The polylines that will each polygon.
:param out_fc: The output with the split geometry added to it.
"""

import os
import sys
import arcpy

arcpy.SetProgressor("default", "Firing up script...")

to_cut = arcpy.GetParameterAsText(0)
cutter = arcpy.GetParameterAsText(1)
out_fc = arcpy.GetParameterAsText(2)

spatialref = arcpy.Describe(to_cut).spatialReference
polygons = []
lines = []
slices = []
gcount = 0

pcount = 0
with arcpy.da.SearchCursor(to_cut, ["SHAPE@", "OID@"]) as pcursor:
    for prow in pcursor:
        arcpy.SetProgressorLabel("Generating slices: {0} rows complete".format(str(pcount)))
        polygon = prow[0]
        polyid = prow[1]
        polygons.append((polygon, polyid))
        pcount += 1
del pcursor

lcount= 0
with arcpy.da.SearchCursor(cutter, ["SHAPE@", "OID@"]) as lcursor:
    for lrow in lcursor:
        line = lrow[0]
        lineid = lrow[1]
        lines.append((line, lineid))
        lcount += 1
del lcursor

def cut_geometry():
    global polygons
    global lines
    global slices
    global gcount
    for eachpoly, eachpolyid in polygons:
        iscut = False
        for eachline, eachlineid in lines:
            if eachline.crosses(eachpoly):
                try:
                    slice1, slice2 = eachpoly.cut(eachline)
                    polygons.insert(0, (slice1, eachpolyid))
                    polygons.insert(0, (slice2, eachpolyid))
                    iscut = True
                except RuntimeError:
                    continue
        if iscut == False:
            slices.append((eachpoly, eachpolyid))
            gcount += 1
            if gcount % 10 == 0:
                arcpy.SetProgressorLabel("Cutting polygons: {0} rows complete".format(str(gcount)))
        polygons.remove((eachpoly, eachpolyid))

while polygons:
    cut_geometry()

arcpy.SetProgressorLabel("Creating output feature class")
arcpy.CreateFeatureclass_management(os.path.dirname(out_fc), os.path.basename(out_fc), "POLYGON",
                                    spatial_reference = spatialref)
arcpy.AddField_management(out_fc, "SOURCE_OID", "LONG")
scount = 0
with arcpy.da.InsertCursor(out_fc, ["SHAPE@", "SOURCE_OID"]) as icursor:
    for eachslice in slices:
        if scount % 10 == 0:
            arcpy.SetProgressorLabel("Inserting slices: {0} rows complete".format(str(scount)))
        icursor.insertRow(eachslice)
        scount += 1
del icursor

arcpy.SetProgressorLabel("Deleting duplicate slices")
shapefieldname = arcpy.Describe(out_fc).shapeFieldName
arcpy.DeleteIdentical_management(out_fc, [shapefieldname, "SOURCE_OID"])

arcpy.ResetProgressor()

Thanks!

Answer 4

Another code sample. This one saves attributes from input polygon layer.

The main idea stills the same - convert single polygon to lines, create temporary layer of this polygon's boundary and intersecting it lines, then convert it back to polygons, clip by initial polygon and write it back.

Also there is an option to include SQL-query for the line layer.

# coding: utf-8

import arcpy
import os
import sys
import time

def splitPolygonsWithLines(Poly, Lines, LinesQuery="", outPoly=""):
    inputPoly=Poly
    inputLines=Lines
    query=LinesQuery

    inputPolyName=os.path.basename(inputPoly)
    inputLinesName=os.path.basename(inputLines)

    parDir=os.path.abspath(inputPoly+"\..")
    if outPoly=="":
        outputPolyName=os.path.splitext(inputPolyName)[0]+u"_Split"+os.path.splitext(inputPolyName)[1]
        outputPoly=os.path.join(parDir,outputPolyName)
    else:
        outputPolyName=os.path.basename(outPoly)
        outputPoly=outPoly

    sr=arcpy.Describe(inputPoly).spatialReference
    fieldNameIgnore=["SHAPE_Area", "SHAPE_Length"]
    fieldTypeIgnore=["OID", "Geometry"]

    #############################################################################################################################
    arcpy.CreateFeatureclass_management (parDir, outputPolyName, "POLYGON", "", "", "", sr)

    arcpy.AddField_management(outputPoly, "OLD_OID", "LONG")
    for field in arcpy.ListFields(inputPoly):
        if (field.type not in fieldTypeIgnore and field.name not in fieldNameIgnore):
            arcpy.AddField_management (outputPoly, field.name, field.type)


    arcpy.MakeFeatureLayer_management(inputLines,inputLinesName+"Layer",query)
    arcpy.MakeFeatureLayer_management(inputPoly,inputPolyName+"Layer")

    arcpy.SelectLayerByLocation_management(inputPolyName+"Layer","INTERSECT",inputLinesName+"Layer","","NEW_SELECTION")
    arcpy.SelectLayerByAttribute_management(inputPolyName+"Layer", "SWITCH_SELECTION")

    fieldmappings = arcpy.FieldMappings()
    for field in arcpy.ListFields(inputPoly):
        if (field.type not in fieldTypeIgnore and field.name not in fieldNameIgnore):
            fm=arcpy.FieldMap()
            fm.addInputField(outputPoly, field.name)
            fm.addInputField(inputPolyName+"Layer", field.name)

            fm_name = fm.outputField
            fm_name.name = field.name
            fm.outputField = fm_name

            fieldmappings.addFieldMap (fm)

    fm=arcpy.FieldMap()
    fm.addInputField(outputPoly, "OLD_OID")
    fm.addInputField(inputPolyName+"Layer", "OBJECTID")

    fm_name = fm.outputField
    fm_name.name = "OLD_OID"
    fm.outputField = fm_name

    fieldmappings.addFieldMap (fm)

    arcpy.Append_management(inputPolyName+"Layer", outputPoly, "NO_TEST", fieldmappings)

    polySelect=arcpy.SelectLayerByLocation_management(inputPolyName+"Layer","INTERSECT",inputLinesName+"Layer","","NEW_SELECTION")
    lineSelect=arcpy.SelectLayerByLocation_management(inputLinesName+"Layer","INTERSECT",inputPolyName+"Layer","","NEW_SELECTION")
    #############################################################################################################################

    fields=[f.name for f in arcpy.ListFields(inputPoly) if (f.type not in fieldTypeIgnore and f.name not in fieldNameIgnore)]
    fields.append("SHAPE@")
    totalFeatures=int(arcpy.GetCount_management(polySelect).getOutput(0))

    count=0
    timePrev=time.time()
    with arcpy.da.SearchCursor(polySelect,["OID@"]+fields) as curInput:
        for rowInput in curInput:
            linesTemp=arcpy.SelectLayerByLocation_management(lineSelect,"INTERSECT",rowInput[-1],"","NEW_SELECTION")
            geometry=arcpy.CopyFeatures_management(linesTemp,arcpy.Geometry())
            geometry.append(rowInput[-1].boundary())
            arcpy.FeatureToPolygon_management (geometry, "in_memory\polygons_init")
            arcpy.Clip_analysis ("in_memory\polygons_init", rowInput[-1], "in_memory\polygons_clip")
            with arcpy.da.SearchCursor("in_memory\polygons_clip","SHAPE@") as curPoly:
                newGeom=[]
                for rowP in curPoly:
                    if not rowP[0].disjoint(rowInput[-1]):
                        newGeom.append(rowP[0])
            arcpy.Delete_management("in_memory")

            with arcpy.da.InsertCursor(outputPoly, ["OLD_OID"]+fields) as insCur:
                for geom in newGeom:
                    insertFeature=[r for r in rowInput[:-1]]
                    insertFeature.append(geom)
                    insCur.insertRow(insertFeature)
            count+=1
            if int(time.time()-timePrev)%5==0 and int(time.time()-timePrev)>0:
                timePrev=time.time()
                arcpy.AddMessage("\r{0}% done, {1} features processed".format(int(count*100/totalFeatures),int(count)))

def main():
    arcpy.env.overwriteOutput = True
    arcpy.env.XYTolerance = "0.1 Meters"

    inputPoly = arcpy.GetParameterAsText(0) # required
    inputLines = arcpy.GetParameterAsText(1) # required
    linesQuery = arcpy.GetParameterAsText(2) # optional
    outPoly = arcpy.GetParameterAsText(3) # optional

    if inputPoly=="":
        inputPoly="D:/Projects/MOVE/MOVE.gdb/poly"
    if arcpy.Exists(inputPoly):
        arcpy.AddMessage("Input polygons: "+inputPoly)
    else:
        arcpy.AddError("Input polygons layer %s is invalid" % (inputPoly))

    if inputLines=="":
        inputLines="D:/Projects/MOVE/MOVE.gdb/line"
    if arcpy.Exists(inputLines):
        arcpy.AddMessage("Input lines: "+inputPoly)
    else:
        arcpy.AddError("Input lines layer %s is invalid" % (inputLines))

    if linesQuery=="":
        arcpy.AddMessage("Performing without query")

    if outPoly == "":
        arcpy.AddMessage("Output will be created at the same location as input polygons layer is.")

    splitPolygonsWithLines(inputPoly, inputLines, linesQuery, outPoly)

if __name__ == "__main__":
    main()