Dividir polígono com linha usando ArcPy?

Ao analisar o Buffer com barreira física usando o ArcGIS for Desktop? , ocorreu-me que não tenho certeza de como alguém poderia usar as ferramentas de geoprocessamento no ArcGIS para dividir um polígono com uma linha programaticamente.

Manualmente, você usaria a ferramenta Cortar polígonos ou a ferramenta Dividir polígonos na barra de ferramentas Topologia , mas como você executaria a mesma tarefa usando as ferramentas de script modelbuilder ou groprocessamento de python?

Logo de cara , penso em todas as ferramentas da caixa de análise, como Union, Identity, etc., mas essas são todas ferramentas Polygon-Polygon, NÃO ferramentas Polygon-Line. Até a ferramenta Dividir é Polígono-Polígono.

Alguma ideia?

Usando o ET Geowizard, você pode acessar o código da ferramenta Split Polygons with Polylines:

Aqui está o link para o script .

Como alternativa, você pode usar o ArcObjects para fazer isso:

Cortar fragmento de polígono

Você também pode usar o método de buffer de um lado descrito aqui .

Depois disso, acabei criando minha própria ferramenta ModelBuilder. Esqueci essa pergunta e postei minha solução em outra pergunta semelhante . Para ser completo, este é um repost da resposta:

Eu pensei que deveria haver uma maneira de fazer isso, então criei o que acredito ser uma solução muito boa. Eu o publiquei no site ArcGIS Resources na Comunidade -> Técnico -> Análise e geoprocessamento -> Análise -> Galeria.

A ferramenta é chamada de dividir polígonos com linhas e requer uma licença do ArcInfo devido a algumas das ferramentas usadas no modelo. Essencialmente, o que eu fiz foi criar a caixa delimitadora mínima para os polígonos e estender as linhas para eles. Então, usando um pouco de vodu do ModelBuilder, consegui transformar o trabalho de linha em polígonos, que depois usei Identity para dividir os polys originais.

Teste-o e veja se funciona para você. Nos meus testes (limitados), ele preservou os atributos dos polígonos originais e dividiu apenas os polígonos existentes.

Se você quiser sair do ArcGIS, use geom.splitpolysbylines .

Pessoalmente, nunca o usei em um programa, mas acho que você pode acessar esta linha de comando com python, consulte a ajuda para obter mais detalhes.

se você não tiver problemas de alta precisão, você armazenará em buffer a linha com a distância mínima, por exemplo (0,002, acho que deve ser superior à precisão da classe de recursos) e, em seguida, aplique uma ferramenta de apagamento ao polígono pela linha em buffer.

Código arcpy atualizado para dividir polígonos na direção horizontal ou vertical usando valores percentuais

'''-------------------------------------------------------------------------------------------
 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. [email protected]

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 - [email protected]
-------------------------------------------------------------------------------------------'''

# 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"