11

In looking at Buffering with physical barrier using ArcGIS for Desktop?, it occurred to me that I am not sure how one would go about using geoprocessing tools in ArcGIS to split a polygon with a line programatically.

Manually, you would use the Cut Polygons tool or the Split Polygons tool on the Topology toolbar, but how would you accomplish the same task using modelbuilder or python groprocessing scripting tools?

Right off the bat I think of all of the tools in the Analysis toobox like Union, Identity, etc, but those are all Polygon-Polygon tools, NOT Polygon-Line tools. Even the Split tool is Polygon-Polygon.

Any ideas?

6

Using ET Geowizard you can access the code for the Split Polygons with Polylines tool:

enter image description here

Here is the link to the script.

Alternatively, you can use ArcObjects to do this:

Cut Polygon Snippet

You may also use the one side buffer method described here.

4

After the fact, I ended up creating my own ModelBuilder tool. I had forgot about this question and posted my solution to another similar question. For completeness, this is a repost of the answer:

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

4

If you want to go outside of ArcGIS then use geom.splitpolysbylines.

Personally i have never use it in a programme but i think you can access this commondline with python, please see the help for more details.

3

if you haven't some high accuracy issues , you would buffer the line with the minimum distance for eg (0.002 i think that should be superior than the accuracy of you feature class) , then apply a erase tool to polygon by the buffered line.

  • Good idea, which I originally considered too, but this would actually result in either 3 polygons from the original (left side of buffered line, sliver polygons in the middle, and right side of original polygon). I guess you qualified that with "high accuracy issues", but I don't think this workaround qualify's as a true "solution" for this question. – RyanDalton Jul 30 '12 at 23:14
  • sorry i've made a mistake, i meant the erase tool , i've updated to the answer – geogeek Jul 30 '12 at 23:18
  • yes i see, that's weird that we lack this tool from the toolbox, i wish that we could find a better workaround – geogeek Jul 30 '12 at 23:20
0

Updated arcpy code to split polygons in horizontal or vertical direction using percent values

'''-------------------------------------------------------------------------------------------
 Tool Name:   Polygon Bisector in percent parts
 Source Name: polygonbisector.py and splitPolygonsWithLines
 Version:     ArcGIS 10.0
 Author:      ESRI, Inc.
 Required Arguments:
              Input Features (Feature Layer)
              Output Feature Class (Feature Class)
 Optional Arguments:
              Axis (X|Y)
              Group Field(s) (Field)
              Acceptable Error Percent (Double)

 Description: Computes a line that bisects, or divides in Commulative Percent value, a polygon area along a line
              of constant latitude or longitude. each successive input polygon's area will be
              on either side of the bisecting line.
https://www.arcgis.com/home/item.html?id=9aadb577ccb74f0e88b13a0e3643ca4d
Credits (Attribution)
Esri, Inc. dflater@esri.com

http://www.arcgis.com/home/item.html?id=cd6b2d45df654245b7806a896670a431
Split Polygons using Line features
Shapefile by daltongis
Credits (Attribution)
GIS.StackExchange.com

updated by : Sham Davande, GIS Analyst - sham.davande@gmail.com
-------------------------------------------------------------------------------------------'''

# Import system modules
import arcpy
from arcpy import env
import os
import sys
import math
import time
arcpy.env.overwriteOutput = True


# Change input polygon feature class name and path here
split_poly = "D:\\temp\\polygon2.shp"

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
#   Split percent   Commulative Percent value for python code
#   3.7         3.7
#   25.5            29.2
#   20.8            50.0
#   10.5            60.5
#   31.5            92.0
#   8.0 

# its cumulative percent values to split polygon with N number percent parts
list = ["3.7","29.2","50.0","60.5","92.0"]

# Set local variables
out_folder_path = "C:\\" 
out_name = "temp71"
# Set workspace
env.workspace = "C:\\temp71"

# Execute CreateFolder
arcpy.CreateFolder_management(out_folder_path, out_name)

# Set local variables
out_folder = "C:\\temp71" 
out_name = "NewGdb.gdb"
out_name2 = "C:\\temp71\\NewGdb.gdb"
geometry_type = "POLYLINE"
template = ""
has_m = "DISABLED"
has_z = "DISABLED"

# Execute CreateFileGDB
arcpy.CreateFileGDB_management(out_folder, out_name)

percent_lines = "C:\\temp71\\NewGdb.gdb\\line1"
line2 = "C:\\temp71\\NewGdb.gdb\\line2"
line3 = "C:\\temp71\\NewGdb.gdb\\line3"
allLines = "C:\\temp71\\NewGdb.gdb\\all_lines"
outFeatureClass1 = "polygon"
outFeatureClass2 = "C:\\temp71\\NewGdb.gdb\\polygon"
expression = ""

#############################################
# No hard coded input files path mentioned in the code below here after
# and don't want interfere somebodies C:\Temp folder, so C:\temp71 folder automatically created by code below. 
#############################################


file_prj1 = os.path.splitext(split_poly)[0]
file_prj2 = str(file_prj1 + ".prj")
out_coordinate_system = str(file_prj2)

# Execute CreateFeatureclass
arcpy.CreateFeatureclass_management(out_name2,"all_lines", "POLYLINE", "", "DISABLED", "DISABLED", out_coordinate_system)
# Execute FeatureClassToFeatureClass
arcpy.FeatureClassToFeatureClass_conversion(split_poly, out_name2,
                                            outFeatureClass1, expression)

for i in list:
    if i > 0:
        percent_area = float(i)

        ith_part = 100/percent_area

        print "Generating a polyline to split a polygon into two horizontal parts of", percent_area, "% and",100-percent_area, "% areas"
        print ith_part, "ith part of polygon"

        schemaType = "NO_TEST"
        fieldMappings = ""
        subtype = ""

        # Main function, all functions run in GravityModel
        def PolygonBisector(in_features, out_fc, axis="x", groupfields=[], error=0.001):
        # Error if sufficient license is not available
            if arcpy.ProductInfo().lower() not in ['arcinfo']:
                arcpy.AddError("An ArcInfo/Advanced license is required.")
                sys.exit()

            # Set geoprocessing environments
            arcpy.env.overwriteOutput = True
            arcpy.env.qualifiedFieldNames = False

            shapefield = arcpy.Describe(in_features).shapeFieldName
            rounder = GetRounder(in_features)

        # If group fields are specified, dissolve by them
            if groupfields:
                in_features = arcpy.management.Dissolve(in_features, "in_memory/grouped", groupfields)
            else:
                groupfields = [arcpy.Describe(in_features).OIDFieldName]
            fields = [shapefield] + groupfields

        # Create output feature class and set up cursor
            icur = irow = scur = None
            arcpy.management.CreateFeatureclass(os.path.dirname(out_fc), os.path.basename(out_fc), "POLYLINE", "", "", "", arcpy.Describe(in_features).spatialReference)
            arcpy.management.AddField(out_fc, "Group_", "TEXT", "", "", "", "Group: {0}".format(", ".join(groupfields)))
            icur = arcpy.InsertCursor(out_fc)
            scur = arcpy.SearchCursor(in_features, "", "", ";".join(fields))
            count = int(arcpy.management.GetCount(in_features).getOutput(0))
            arcpy.SetProgressor("step", "Processing polygons...", 0, count, 1)
            bigi = 1

            # Begin processing
            try:
                for row in scur:
                    minx = miny = float("inf")
                    maxx = maxy = float("-inf")
                    totalarea = 0
                    feat = row.getValue(shapefield)
                    totalarea = row.getValue(shapefield).area
                    group = []
                    for field in groupfields:
                        group.append(str(row.getValue(field)))
                    partnum = 0
                    # Get the min and max X and Y
                    for part in feat:
                        for point in feat.getPart(partnum):
                            if point:
                                minx = point.X if point.X < minx else minx
                                miny = point.Y if point.Y < miny else miny
                                maxx = point.X if point.X > maxx else maxx
                                maxy = point.Y if point.Y > maxy else maxy
                        partnum += 1

                    # Process the polygon
                    # Some variables
                    conditionmet = False
                    difference = 0
                    lastdifference = float("inf")
                    differences = {}
                    itys = {}
                    i = 1
                    strike = 0

                    # The starting bisector (half the distance from min to max)
                    if axis == "x":
                        ity = (miny + maxy)/2.0
                    else:
                        ity = (minx + maxx)/2.0

                    while not conditionmet:
                        # Construct a line through the middle
                        if axis == "x":
                            line = MakeBisector(minx, maxx, ity, in_features, axis)
                        else:
                            line = MakeBisector(miny, maxy, ity, in_features, axis)
                        # The FeatureToPolygon function does not except a geometry object, so make a temporary feature class
                        templine = arcpy.management.CopyFeatures(line, "in_memory/templine")
                        temppoly = arcpy.management.CopyFeatures(feat, "in_memory/temppoly")
                        # Intersect then Feature To Polygon
                        bisected = arcpy.management.FeatureToPolygon([temppoly, templine], "in_memory/bisected")
                        clip = arcpy.analysis.Clip(bisected, in_features, "in_memory/clip")

                        # Group bisected polygons according to above or below the bisector
                        arcpy.management.AddField(clip, "FLAG", "SHORT")
                        ucur = arcpy.UpdateCursor(clip, "", "")
                        flag = 0
                        try:
                            for urow in ucur:
                                ufeat = urow.getValue(arcpy.Describe(clip).shapeFieldName)
                                partnum = 0
                                for upart in ufeat:
                                    for upoint in ufeat.getPart(partnum):
                                        if upoint:
                                            if axis == "x":
                                                if round(upoint.Y, rounder) > round(ity, rounder):
                                                    flag = 1
                                                    break
                                                elif round(upoint.Y, rounder) < round(ity, rounder):
                                                    flag = -1
                                                    break
                                            else:
                                                if round(upoint.X, rounder) > round(ity, rounder):
                                                    flag = 1
                                                    break
                                                elif round(upoint.X, rounder) < round(ity, rounder):
                                                    flag = -1
                                                    break
                                    partnum += 1
                                urow.setValue("FLAG", flag)
                                ucur.updateRow(urow)
                        except:
                            raise
                        finally:
                            if ucur:
                                del ucur

                        # Check if the areas are halved
                        dissolve = arcpy.management.Dissolve(clip, "in_memory/dissolve", "FLAG")
                        scur2 = arcpy.SearchCursor(dissolve)
                        try:
                            for row2 in scur2:
                                firstarea = row2.getValue(arcpy.Describe(dissolve).shapeFieldName).area
                                firstflag = row2.getValue("FLAG")
                                break
                        except:
                            raise
                        finally:
                            if scur2:
                                del scur2

                        difference = abs(firstarea - (totalarea/ith_part))
                        differences[i] = difference
                        itys[i] = ity
                        print round(100*(difference/(totalarea/ith_part)),5)
                        #arcpy.AddWarning(round(100*(difference/(totalarea/ith_part)),5))
                        # Stop if tolerance is achieved
                        if (difference/(totalarea/ith_part))*100 <= error:
                            conditionmet = True
                            break
                        # Moving the line in the wrong direction? due to coordinate system origins or over-compensation
                        if difference > lastdifference:
                            firstflag = firstflag*-1.0
                        # If we're not improving
                        if abs(difference) > min(differences.values()):
                            strike+=1
                        # Or if the same values keep appearing
                        if differences.values().count(difference) > 3 or strike >=3:
                            arcpy.AddWarning("Tolerance could not be achieved. Output will be the closest possible.")
                            # Reconstruct the best line
                            if axis == "x":
                                line = MakeBisector(minx, maxx, itys[min(differences,key = lambda a: differences.get(a))], in_features, axis)
                            else:
                                line = MakeBisector(miny, maxy, itys[min(differences,key = lambda a: differences.get(a))], in_features, axis)
                            break
                        # Otherwise move the bisector so that the areas will be more evenly split
                        else:
                            if firstflag == 1:
                                if axis == "x":
                                    ity = ((ity-miny)/((totalarea/ith_part)/firstarea)) + miny
                                else:
                                    ity = ((ity-minx)/((totalarea/ith_part)/firstarea)) + minx
                            elif firstflag == -1:
                                if axis == "x":
                                    ity = ((ity-miny)*math.sqrt((totalarea/ith_part)/firstarea)) + miny
                                else:
                                    ity = ((ity-minx)*math.sqrt((totalarea/ith_part)/firstarea)) + minx
                            lastdifference = difference
                        i +=1
                    irow = icur.newRow()
                    irow.setValue(arcpy.Describe(out_fc).shapeFieldName, line)
                    irow.setValue("Group_", ", ".join(group))
                    icur.insertRow(irow)
                    arcpy.SetProgressorPosition()
                    arcpy.AddMessage("{0}/{1}".format(bigi, count))
                    bigi +=1

            except:
                if arcpy.Exists(out_fc):
                    arcpy.management.Delete(out_fc)
                raise
            finally:
                if scur:
                    del scur
                    if icur:
                        del icur
                if irow:
                    del irow
                for data in ["in_memory/grouped", temppoly, templine, clip, bisected, dissolve]:
                    if data:
                        try:
                            arcpy.management.Delete(data)
                        except:
                            ""

        def MakeBisector(min,max,constant, templatefc, axis):
            if axis == "x":
                array = arcpy.Array()
                array.add(arcpy.Point(min, constant))
                array.add(arcpy.Point(max, constant))
            else:
                array = arcpy.Array()
                array.add(arcpy.Point(constant, min))
                array.add(arcpy.Point(constant, max))
            line = arcpy.Polyline(array, arcpy.Describe(templatefc).spatialReference)
            return line

        def GetRounder(in_features):
            try:
                unit = arcpy.Describe(in_features).spatialReference.linearUnitName.lower()
            except:
                unit = "dd"
            if unit.find("foot") > -1:
                rounder = 1
            elif unit.find("kilo") > -1:
                rounder = 3
            elif unit.find("meter") > -1:
                rounder = 1
            elif unit.find("mile") > -1:
                rounder = 3
            elif unit.find("dd") > -1:
                rounder = 5
            else:
                rounder = 3
            return rounder

        # Run the script
        if __name__ == '__main__':
            # Get Parameters
            in_features = str(outFeatureClass2)
            out_fc = percent_lines
            axis = arcpy.GetParameterAsText(2).lower() or "x"
            groupfields = arcpy.GetParameterAsText(3).split(";") if arcpy.GetParameterAsText(3) else []
            error = float(arcpy.GetParameter(4)) if arcpy.GetParameter(4) else 0.001
            out_data = line2
            # Run the main script
            PolygonBisector(in_features, out_fc, axis, groupfields, error)
            arcpy.Copy_management(out_fc, out_data)
            # Use Append tool to move features into single dataset
            arcpy.Append_management(out_data, allLines, schemaType, fieldMappings, subtype)
            print "again generating a polyline to split same polygon"
            print ""

        else:
            print "Error"

print "finished"
print ""
print "Splitting polygon using multiple lines generated for a percent area values provided by you..."
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=outFeatureClass2
    if arcpy.Exists(inputPoly):
        arcpy.AddMessage("Input polygons: "+inputPoly)
    else:
        arcpy.AddError("Input polygons layer %s is invalid" % (inputPoly))

    if inputLines=="":
        inputLines=allLines
    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()
print ""
print "Done"

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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