¿Dividir polígono con línea usando ArcPy?

¿Al mirar Buffering con barrera física usando ArcGIS for Desktop? , se me ocurrió que no estoy seguro de cómo usar las herramientas de geoprocesamiento en ArcGIS para dividir un polígono con una línea mediante programación.

Manualmente, usaría la herramienta Cortar polígonos o la herramienta Dividir polígonos en la barra de herramientas de Topología , pero ¿cómo realizaría la misma tarea utilizando las herramientas de creación de scripts de modelador de Python o de modelbuilder?

De buenas a primeras pienso en todas las herramientas en el cuadro de análisis también como Unión, Identidad, etc., pero esas son todas las herramientas Polígono-Polígono, NO herramientas de Línea Poligonal. Incluso la herramienta Dividir es Polígono-Polígono.

¿Algunas ideas?

Usando ET Geowizard puede acceder al código para la herramienta Dividir polígonos con polilíneas:

Aquí está el enlace al guión .

Alternativamente, puede usar ArcObjects para hacer esto:

Cortar fragmento de polígono

También puede usar el método de un lado del búfer descrito aquí .

Después del hecho, terminé creando mi propia herramienta ModelBuilder. Me había olvidado de esta pregunta y publiqué mi solución a otra pregunta similar . Para completar, este es un reenvío de la respuesta:

Pensé que debía haber una manera de hacer esto, así que creé lo que creo que es una solución bastante buena. Lo publiqué en el sitio de Recursos de ArcGIS en la Comunidad-> Técnico-> Análisis y Geoprocesamiento-> Análisis-> Galería.

La herramienta se llama Split Polygons With Lines y requiere una licencia ArcInfo debido a algunas de las herramientas utilizadas en el modelo. Esencialmente, lo que hice fue crear el cuadro de límite mínimo para los polígonos y extenderles las líneas. Entonces, usando un poco de vudú ModelBuilder, pude convertir la línea en polígonos, que luego usé Identity para dividir los polys originales.

Pruébelo y vea si funciona para usted. En mis pruebas (limitadas) conservó los atributos de los polígonos originales y dividió solo los polígonos existentes.

Si desea salir de ArcGIS, use geom.splitpolysbylines .

Personalmente, nunca lo he usado en un programa, pero creo que puede acceder a este commondline con python, consulte la ayuda para obtener más detalles.

si no tiene problemas de alta precisión, debería tamponar la línea con la distancia mínima, por ejemplo (0.002, creo que debería ser superior a la precisión de su clase de entidad), luego aplicar una herramienta de borrado al polígono por la línea tamponada.

Código de arcpy actualizado para dividir polígonos en dirección horizontal o vertical utilizando valores de porcentaje

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