¿Cómo crear un polígono conectando los extremos de líneas múltiples utilizando python?

Question

Estoy tratando de averiguar cómo crear un polígono que conecta a todos los puntos extremos de un shapefile que contiene un conjunto de polyilnes con pythonscript en ArcGIS, estoy teniendo problemas para hacerlo, ya que el orden de los nodos en el polígono es importante. Quiero lograr el gris polígono, en la imagen de las líneas de color verde

Answer 1

PASOS:

Calcular las secciones de centro de los puntos:

Construido sus Euclidiana mínimo árbol de expansión, disolverlo y calcular buffer, distancia equivalente a la mitad de la menor longitud de la sección:

Crear extremo de la sección de puntos y calcular su pk (distancia a lo largo de la línea) en el límite de la memoria intermedia (polilínea cerrada de la versión de búfer):

Ordenar los puntos finales en orden ascendente mediante pk campo. Los puntos a continuación marcados por sus FID:

Crear polígono de conjunto ordenado de puntos:

Secuencia de comandos:

import arcpy, traceback, os, sys,time
from heapq import *
from math import sqrt
import itertools as itt
from collections import defaultdict

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)
    # MST by PRIM's
    def prim( nodes, edges ):
        conn = defaultdict( list )
        for n1,n2,c in edges:
            conn[ n1 ].append( (c, n1, n2) )
            conn[ n2 ].append( (c, n2, n1) )
        mst = []
        used = set( nodes[ 0 ] )
        usable_edges = conn[ nodes[0] ][:]
        heapify( usable_edges )

        while usable_edges:
            cost, n1, n2 = heappop( usable_edges )
            if n2 not in used:
                used.add( n2 )
                mst.append( ( n1, n2, cost ) )

                for e in conn[ n2 ]:
                    if e[ 2 ] not in used:
                        heappush( usable_edges, e )
        return mst        


    mxd = arcpy.mapping.MapDocument("CURRENT")
    SECTIONS=arcpy.mapping.ListLayers(mxd,"SECTION")[0]
    PGONS=arcpy.mapping.ListLayers(mxd,"RESULT")[0]
    d=arcpy.Describe(SECTIONS)
    SR=d.spatialReference

    cPoints,endPoints,lMin=[],[],1000000
    with arcpy.da.SearchCursor(SECTIONS, "Shape@") as cursor:
        # create centre and end points
        for row in cursor:
            feat=row[0]
            l=feat.length
            lMin=min(lMin,feat.length)
            theP=feat.positionAlongLine (l/2).firstPoint
            cPoints.append(theP)
            theP=feat.firstPoint
            endPoints.append(theP)
            theP=feat.lastPoint
            endPoints.append(theP)

        arcpy.AddMessage('Computing minimum spanning tree')
        m=len(cPoints)
        nodes=[str(i) for i in range(m)]
        p=list(itt.combinations(range(m), 2))
        edges=[]
        for f,t in p:
            p1=cPoints[f]
            p2=cPoints[t]
            dX=p2.X-p1.X;dY=p2.Y-p1.Y
            lenV=sqrt(dX*dX+dY*dY)
            edges.append((str(f),str(t),lenV))
        MST=prim(nodes,edges)

        mLine=[]
        for edge in MST:
            p1=cPoints[int(edge[0])]
            p2=cPoints[int(edge[1])]
            mLine.append([p1,p2])
        pLine=arcpy.Polyline(arcpy.Array(mLine),SR)

        # create buffer and compute chainage
        buf=pLine.buffer(lMin/2)
        outLine=buf.boundary()
        chainage=[]
        for p in endPoints:
            measure=outLine.measureOnLine(p)
            chainage.append([measure,p])
        chainage.sort(key=lambda x: x[0])

        # built polygon
        pGon=arcpy.Array()
        for pair in chainage:
            pGon.add(pair[1])
        pGon=arcpy.Polygon(pGon,SR)
        curT = arcpy.da.InsertCursor(PGONS,"SHAPE@")
        curT.insertRow((pGon,))
        del curT
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()

Sé que es una bicicleta, pero es mi propio y me gusta

Answer 2

He puesto esta solución para QGIS aquí porque es software libre y fácil de implementar. He considerado que sólo el derecho "rama" de la polilínea capa vectorial; como se puede observar en la imagen siguiente (12 de características en la tabla de atributos):

El código (algoritmo en una línea de python comprensión de lista), para ejecutar en la Consola de Python de QGIS, es:

layer = iface.activeLayer()

features = layer.getFeatures()

features = [feature for feature in features]

n = len(features)

geom = [feature.geometry().asPolyline() for feature in features ]

#multi lines as closed shapes
multi_lines = [[geom[i][0], geom[i][1], geom[i+1][1], geom[i+1][0], geom[i][0]]
               for i in range(n-1)]

#multi polygons
mult_pol = [[] for i in range(n-1)]

for i in range(n-1):
    mult_pol[i].append(multi_lines[i])

#creating a memory layer for multi polygon
crs = layer.crs()
epsg = crs.postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "polygon",
                           "memory")

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

mem_layer.startEditing()

#Set features
feature = [QgsFeature() for i in range(n-1)]

for i in range(n-1):
    #set geometry
    feature[i].setGeometry(QgsGeometry.fromPolygon(mult_pol[i]))
    #set attributes values
    feature[i].setAttributes([i])
    mem_layer.addFeature(feature[i], True)

#stop editing and save changes
mem_layer.commitChanges()

Después de ejecutar el código:

se realizó un polígono capa de memoria (con 11 funciones en su tabla de atributos). Funciona muy bien.

Answer 3

Puede seleccionar los extremos que van a participar en un polígono, crear un TIN a partir de sólo los puntos. Convertir la LATA a los polígonos, disolver los polígonos. El truco para la automatización de este proceso es decidir que apunta a contribuir para cada polígono. Si tienes líneas con direcciones válidas, y esas líneas, todos comparten algún atributo común podría escribir una consulta para la exportación decir, al final vértices usando la línea de vértices a los puntos, a continuación, seleccione por atributo aquellos puntos que tienen en común el valor del atributo.
Mejor sería extraer/seleccionar los puntos, lea las x , y valores mediante un cursor, el uso de la x, y valores para escribir un nuevo polígono. No puedo ver una foto adjunta en su puesto, pero si el punto de orden de asuntos, a continuación, una vez que usted tiene la x, y los valores almacenados en una lista de Python, ordenarlas. http://resources.arcgis.com/EN/HELP/MAIN/10.1/index.html#//002z0000001v000000

Answer 4

La expansión en @iant comentario, el más cercano a la geometría de la instantánea es la forma alfa (alfa casco) de los extremos. Afortunadamente, muchos bien recibido hilos ya han sido respondidas en el SIG SE. Por ejemplo:

• Crear Polígono Complejas desde el Punto de Capa con el único Límite de Puntos en ArcGIS

• ¿Cuáles son Definición, Algoritmos y Soluciones Prácticas para Casco Cóncavo?

Para solucionar su problema, el primer uso de la Característica De Punto para extraer los puntos finales. A continuación, utilice la herramienta de python desde este Enlace para calcular el casco cóncavo.