¿Herramienta automatizada/geoprocesada para cortar, recortar o cortar polígonos utilizando polilíneas en ArcGIS Desktop?

Question

Buscando un método más simple aquí.

Estoy tratando de dividir / cortar / rebanar polígonos existentes mediante el uso de polilíneas existentes. Un ejemplo sería dividir una masa de agua o una parcela de tierra en el punto donde un puente/carretera la cruza. Pero la red de carreteras no tiene por qué dar lugar a un bucle cerrado.

Como las polilíneas no están necesariamente interconectadas o son continuas, crear un polígono a partir de ellas no es una opción (lo que elimina el uso de la herramienta de división). Además, he intentado utilizar una topología con la geometría, pero sigue fallando, probablemente debido a la geometría grande/compleja.

Flujo de trabajo actual: He logrado esto con la herramienta Feature to Polygon, combinando las líneas y los polígonos juntos, pero crea polígonos adicionales no deseados en cualquier lugar donde se cierra un bucle. He creado una máscara de los polígonos originales, y la he utilizado con la herramienta superponer->borrar para eliminar los polígonos no deseados. Esto todavía deja algo de geometría no deseada (sobre todo astillas), pero es algo factible.

Esto parece una forma extremadamente enrevesada y redonda de llevar a cabo (lo que parece que debería ser) una tarea muy sencilla.

Aparte de las ediciones manuales, o de la utilización de una topología, ¿hay alguna herramienta que pueda lograr esto en un solo paso?

Utilizando: ArcMap \ ArcInfo Desktop 10 SP5

Editar 1: En mi caso, no es realmente carreteras como se mencionó anteriormente. Tengo polígonos de agua para una zona costera, y necesito dividir los polígonos donde las presas de roca o diques se han puesto a través de los cursos de agua. Los cuales normalmente no están interconectados.

Los polígonos del agua se han simplificado y reparado, hasta el punto de que ya no llamaría a los datos "sucios", sólo complejos y grandes. Pero he conseguido que la solución mencionada anteriormente funcione para este caso.

Estoy buscando más "en general" una herramienta que pueda simplemente dividir polígonos usando polilíneas.

Editar 2:

Mapperz: Gracias por la sugerencia del Model Builder. Voy a utilizarlo como solución provisional, por ahora.

Jakub: Gracias por la sugerencia. No me opongo a una solución programática o a desarrollar una herramienta personalizada, aunque nunca he hecho una antes. Tengo experiencia en programación, pero no en conjunto con Arc. Sin embargo, preferiría algo que corte directamente la geometría, en lugar de seguir la logística indirecta anterior. En teoría, eso debería reducir las brechas resultantes, ya que no sería objeto de múltiples iteraciones de craqueo / agrupación. Aunque no estoy seguro de que eso sea tan fácil o incluso posible.

Edición 3: Estoy buscando algo que funcione como en la izquierda de la imagen.

Answer 1

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

La herramienta se llama Dividir polígonos con líneas y requiere una licencia de ArcInfo debido a algunas de las herramientas utilizadas dentro del modelo. Básicamente, lo que hice fue crear el cuadro delimitador mínimo para los polígonos y extender las líneas hasta ellos. Así que usando un poco de vudú de ModelBuilder, pude convertir las líneas en polígonos, que luego usé Identity para dividir los polígonos originales.

Por favor, pruébalo y comprueba si te funciona. En mis pruebas (limitadas) conservó los atributos de los polígonos originales, y dividió sólo los polígonos existentes.

Answer 2

Creo que "Feature to Polygon" hace exactamente lo que necesitas. Puede introducir una combinación de clases de características de polígono y polilínea. Los polígonos de salida se dividen en cada polilínea. Se requiere una licencia de ArcInfo, que usted tiene. Probado en 10.0.

Asegúrese de tener un campo poblado, para todas las características, antes de ejecutar "Feature to Polygon". Los nuevos polígonos rellenados tendrán todos los campos en blanco. Otro método es seleccionar espacialmente los polígonos "cortados" con los polígonos originales. Los polígonos rellenos "cortados" no "tendrán su centro" en los polígonos originales.

Answer 3

Es una posibilidad remota, pero ¿tienen una extensión de "cartografía de producción", "cartografía de defensa", "aviación" o "cartografía marítima" para ArcGIS? Parece que la extensión " Dados Polígonos " logrará lo que necesitas. Nunca he utilizado la herramienta personalmente, pero por la descripción, parece un ganador.

Enlace de búsqueda

---

EDITAR:

Alternativamente, he escrito un algoritmo de dados para usted basado en la discusión en " ¿Optimización del código ArcPy para cortar polígonos? " para cortar/dividir/dividir una clase de característica de polígono utilizando una clase de característica de polilínea. Créditos a los del hilo.

Sólo tienes que hacer una herramienta esqueleto que acepte como parámetros:

1. Clase de rasgo poligonal de entrada (tipo: clase de rasgo, filtro: polígono)
2. Clase de rasgo de polilínea de corte (tipo: clase de rasgo, filtro: polilínea)
3. Clase de característica de polígono de salida (Tipo: Clase de rasgo, Dirección: Salida

Hay dos supuestos que he observado al ejecutar la herramienta:

• Su polilínea debe cruzar el polígono por completo, sin detenerse a mitad de camino ni pasar a una segunda línea dentro del polígono para terminar la travesía.
• Si tiene polígonos de varias partes, los cortes resultantes también serán a veces de varias partes, ya que heredan la partición de sus padres. Los polígonos que se auto-interceptan están bien y no tienen ningún problema.

El código siguiente para guardar como archivo .py:

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

Gracias.

Answer 4

Otro ejemplo de código. Este guarda los atributos de la capa de polígonos de entrada.

La idea principal sigue siendo la misma - convertir un solo polígono en líneas, crear una capa temporal de los límites de este polígono y líneas de intersección, luego convertirlo de nuevo en polígonos, recortar por el polígono inicial y escribirlo de nuevo.

También hay una opción para incluir la consulta SQL para la capa de líneas.

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