16

This question already has an answer here:

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.

Cut Polygon vs Feature to Polygon

marked as duplicate by PolyGeo Oct 18 '16 at 6:06

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • Use can recreate your actions with the Model Builder - help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… - ArcGIS is always convoluted... – Mapperz Aug 22 '12 at 14:47
  • It sounds like you have very dirty/ugly data. It is difficult to automate processing on dirty data. It usually requires some manual cleanup before processing. For example, if your roads layer isn't continuous and has random gaps or missing vertices, it could be rough to chop things up with it. You might want to play with the Extend Line tool to clean up danglers. – Baltok Aug 22 '12 at 15:33
  • I think adding a Cut Polygons tool, that divides an input polygon feature class using an input line feature class, to ArcToolbox would make an excellent ArcGIS Idea. If you submit it be sure to place a link in your original Question to try and attract some upvotes. – PolyGeo Aug 22 '12 at 22:54
  • These are the type of situations where you could develop a custom tool with ArcObjects. I don't have a tool for your particular situation but this can easily be handled programatically. You pretty much have all the logistics figured out you just need to write it in code. I see you want an existing tool so I presume programmatic solution is not acceptable here (?) There would have to be a function that handles the gaps in your lines which would add complexity but not impossible. Sounds like a good custom add-in project. – Jakub Sisak GeoGraphics Aug 23 '12 at 13:32
  • Hi, I can’t add the new answer here. How I can do it? But I had the same task and found a way to solve it directly in ArcGIS by ordinary tools. It needs 3 steps. If it is ok, maybe you can move it as answer here? I suppose this method is very simple and useful. – natalia Apr 17 '17 at 17:00
11

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.

5

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.

  • Not quite. Feature to Polygon fills in any closed geometry with a new polygon. I only want to slice existing polygons. I am currently using this as part of my workaround, though. (See orig question) – brnt Aug 23 '12 at 18:56
  • 1
    This is simple to remedy. 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. – klewis Aug 27 '12 at 16:46
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__ = "jtran20@masonlive.gmu.edu, 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!

  • Nice work! Provided it as an answer to a Geonet post (geonet.esri.com/message/627051#comment-627051) – capie69 Aug 10 '16 at 6:46
  • @capie69 how has this worked on large (>1,000) number of features? I just wanted to sense before trying this out. It seems easy enough to test out that I might just do that instead! – jwx Jan 25 '17 at 16:48
  • @jwx haven't tested with lots of features, but should be fast enough. – capie69 Jan 26 '17 at 8:17
  • From what I recall, it works reasonably well for anything less than ~50,000 features. The progressor will keep tabs on how many features have been processed so you can get a sense of speed. – John Feb 2 '17 at 1:31
  • [There's a python toolbox here that can split polygons by polylines with regular licence][1] [1]: gis.stackexchange.com/questions/124198/… – Tristan Forward Aug 17 '18 at 22:33
0

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.