¿Herramienta automatizada / de geoprocesamiento para cortar, recortar o cortar polígonos usando polilíneas usando ArcGIS Desktop? [duplicar]

16

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

Estoy tratando de dividir / cortar / cortar polígonos existentes mediante el uso de polilíneas existentes. Un ejemplo sería dividir un cuerpo de agua o una parcela de tierra en el punto donde un puente / carretera lo cruza. Pero la red de carreteras no necesariamente da como resultado un circuito cerrado.

Dado que las líneas de polietileno 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 usar una topología con la geometría, pero sigue fallando, probablemente debido a la geometría grande / compleja.

Flujo de trabajo actual: Lo he logrado con la herramienta Característica a polígono, combinando las líneas y los polígonos, pero crea polígonos no deseados adicionales en cualquier lugar donde se cierre un bucle. Creé una máscara de los polígonos originales y la utilicé con la herramienta de superposición y borrar para eliminar los polígonos no deseados. Esto todavía deja algo de geometría no deseada (en su mayoría, astillas), pero es algo viable.

Esto parece ser una forma extremadamente complicada y redondeada de lograr (lo que parece que debería ser) una tarea muy simple.

Aparte de las ediciones manuales o el uso de una topología, ¿existe una herramienta que pueda lograr esto en un solo paso?

Utilizando: ArcMap \ ArcInfo Desktop 10 SP5

Edit 1: En mi caso, no es en realidad carreteras como se mencionó anteriormente. Tengo polígonos de agua para un área costera, y necesito dividir los polígonos donde se han colocado diques o diques rocosos a través de las vías fluviales. Los cuales típicamente no están interconectados.

Los polígonos de agua se han simplificado y reparado, hasta el punto en que ya no llamaría a los datos "sucios", simplemente complejos y grandes. Pero tengo la solución alternativa mencionada anteriormente para este caso.

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

Edit 2:

Mapperz: Gracias por la sugerencia de Model Builder. Voy a usar eso como una solución de interrupción, por ahora.

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

Edición 3: estoy buscando algo para funcionar como en la imagen de la izquierda.

    
pregunta brnt 22.08.2012 - 16:44

4 respuestas

11

Pensé que debía haber una manera de hacer esto, así que creé lo que creo que es una buena solución. Lo publiqué en el sitio ArcGIS Resources y en el sitio Community- > Technical- > Analysis & Geoprocesamiento- > Análisis- > Galería.

La herramienta se llama Dividir polígonos con líneas y requiere una licencia de ArcInfo debido a algunos de Las herramientas utilizadas dentro del modelo. Esencialmente, lo que hice fue crear el cuadro de límite mínimo para los polígonos y extenderles las líneas. Así que usando un voodoo de ModelBuilder, pude convertir las líneas en polígonos, que luego usé Identity para dividir las polis originales.

Por favor, pruébalo y ve si funciona para ti. En mis (limitadas) pruebas, conservó los atributos de los polígonos originales y dividió solo los polígonos existentes.

    
respondido por el RyanDalton 24.08.2012 - 00:51
5

Creo que "Feature to Polygon" hace exactamente lo que necesitas. Puede ingresar una combinación de clases de entidad 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 que se complete, para todas las funciones, antes de ejecutar "Característica en Polígono. Los nuevos Rellenos de Polígonos tendrán todos los campos en blanco. Otro método es seleccionar espacialmente los polígonos" cortados "con los Polígonos originales. Los polis rellenos "cortados" no tendrán "su centro" en los polígonos originales.

    
respondido por el klewis 23.08.2012 - 20:36
2

Esta es una posibilidad remota, pero ¿tiene una extensión de 'Mapeo de producción', 'Mapeo de defensa', 'Aviación' o 'Trazado marítimo' para ArcGIS? Parece que la herramienta " Dice Polygons " logrará lo que necesites. Nunca he usado la herramienta personalmente, pero según la descripción, parece un ganador.

Enlace de búsqueda

EDITAR:

Alternativamente, he escrito un algoritmo de dados basado en la discusión en " ¿Optimizando el código de ArcPy para cortar un polígono? " para cortar / dividir / cortar una clase de entidad poligonal usando una clase de entidad de polilínea. Créditos a los que están en el hilo.

Solo necesitas hacer una herramienta de esqueleto que acepte como parámetros:

  1. Clase de entidad poligonal de entrada (Tipo: Clase de entidad, Filtro: Polígono)
  2.   
  3. clase de entidad de corte de polilínea (tipo: clase de entidad, filtro: polilínea)
  4.   
  5. Clase de entidad poligonal de salida (Tipo: Clase de entidad, Dirección: Salida

Observé dos suposiciones al ejecutar la herramienta:

  • Su polilínea debe cruzar el polígono por completo, sin detenerse a mitad de camino o pasar a una segunda línea dentro del polígono para terminar el recorrido.
  •   
  • Si tiene polígonos de múltiples partes, los segmentos resultantes también serán a veces múltiples, ya que heredan la partición de sus padres. Los polígonos que se intersecan están bien y no tienen problemas.

Código abajo para guardar como archivo .py:

__author__ = "John K. Tran, Tristan Forward"
__contact__ = "[email protected], 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, ["[email protected]", "[email protected]"]) 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, ["[email protected]", "[email protected]"]) 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, ["[email protected]", "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!

    
respondido por el John 03.07.2015 - 19:32
0

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

La idea principal sigue siendo la misma: convertir un solo polígono en líneas, crear una capa temporal del límite de este polígono e intersectar sus líneas, luego convertirlo de nuevo en polígonos, cortar por polígono inicial y escribirlo de nuevo.

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

# 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("[email protected]")
    totalFeatures=int(arcpy.GetCount_management(polySelect).getOutput(0))

    count=0
    timePrev=time.time()
    with arcpy.da.SearchCursor(polySelect,["[email protected]"]+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","[email protected]") 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()
    
respondido por el Serge Norin 24.08.2016 - 14:18

Lea otras preguntas en las etiquetas