Czy tworzysz wielokąt łączący punkty końcowe wielu linii za pomocą ArcPy?

10

Próbuję wymyślić, jak utworzyć wielokąt, który łączy wszystkie punkty końcowe pliku kształtu zawierającego zestaw polilin z pythonscript w ArcGIS, mam problem z tym, ponieważ kolejność węzłów w wielokącie jest ważna. Chcę uzyskać szary wielokąt na obrazie z zielonych linii

Chcę połączyć punkty końcowe zielonych linii, aby utworzyć szary wielokąt bez konieczności robienia tego ręcznie

Amanda
źródło
czy twoje linie mają jakiś atrybut do nadania kolejności?
Ian Turton
najpierw potrzebna jest kolejność zdefiniowana jako @iant, a następnie reguła, czy łączy się punkt końcowy z następnym punktem startowym, czy robi to w inny sposób
Matej,
3
w przeciwnym razie może jakiś kadłub alfa w punktach końcowych?
Ian Turton
Linia ma do pewnego stopnia atrybuty, które nadają im porządek. Mają numer identyfikacyjny, ale w powyższym przykładzie prawy oddział ma numer 1-7, lewy 15-21, a po połączeniu identyfikatory to 22-27
Amanda
1
Możesz się bardzo zbliżyć, a) tworząc NIP, używając linii, b) przekształcając NIP w trójkąty c) wybierając trójkąty, które dzielą granicę z liniami. Na górze masz tylko 1 wielokąt do usunięcia
FelixIP

Odpowiedzi:

11

KROKI:

Oblicz punkty środkowe sekcji: wprowadź opis zdjęcia tutaj

Zbudowali swoje euklidesowe minimalne drzewo opinające, rozpuść je i oblicz bufor, odległość równa połowie najkrótszej sekcji: wprowadź opis zdjęcia tutaj

Utwórz punkty końcowe przekroju i oblicz ich pikietaż (odległość wzdłuż linii) na granicy bufora (wersja buforowa zamkniętej polilinii): wprowadź opis zdjęcia tutaj

Posortuj punkty końcowe w porządku rosnącym za pomocą pola pikowania. Punkty poniżej oznaczone ich FID:

wprowadź opis zdjęcia tutaj

Utwórz wielokąt z uporządkowanego zestawu punktów: wprowadź opis zdjęcia tutaj

Scenariusz:

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

Wiem, że to rower, ale jest mój i podoba mi się

FelixIP
źródło
2

Zamieszczam tutaj to rozwiązanie dla QGIS, ponieważ jest to bezpłatne oprogramowanie i łatwe do wdrożenia. Rozważałem tylko właściwą „gałąź” warstwy wektorowej polilinii; jak widać na następnym obrazku (12 cech w tabeli atrybutów):

wprowadź opis zdjęcia tutaj

Kod (algorytm w jednowierszowym rozumieniu listy python) do działania w konsoli Pythona w QGIS to:

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

Po uruchomieniu kodu:

wprowadź opis zdjęcia tutaj

utworzono wieloboczną warstwę pamięci (z 11 cechami w tabeli atrybutów). Działa ładnie.

Xunilk
źródło
1

Możesz wybrać punkty końcowe, które będą uczestniczyć w wielokącie, utworzyć NIP tylko z tych punktów. Przekształć NIP w wielokąty, rozpuść wielokąty. Sztuką automatyzacji tego procesu jest decydowanie, które punkty mają przyczynić się do każdego wielokąta. Jeśli masz linie z prawidłowymi kierunkami, a wszystkie te linie mają wspólny atrybut, możesz napisać zapytanie do wyeksportowania, powiedz, że wierzchołki końcowe używają wierzchołków linii do punktów, a następnie wybierz według atrybutu te punkty, które mają wspólną wartość atrybutu.
Lepiej byłoby wyodrębnić / wybrać punkty, odczytać wartości x, y za pomocą kursora, użyć wartości x, y do napisania nowego wielokąta. Nie widzę załączonego zdjęcia w twoim poście, ale jeśli kolejność punktów ma znaczenie, to kiedy już masz wartości x, y zapisane na liście Pythona, posortuj je. http://resources.arcgis.com/EN/HELP/MAIN/10.1/index.html#//002z0000001v000000

GBG
źródło
1

Rozszerzając komentarz @iant, najbliższą geometrią do twojego migawki jest kształt alfa (kadłub alfa) punktów końcowych. Na szczęście wiele dobrze odebranych wątków zostało już odebranych na GIS SE. Na przykład:

Aby rozwiązać problem, najpierw użyj opcji Funkcja do punktu, aby wyodrębnić punkty końcowe. Następnie użyj narzędzia Python z tego łącza, aby obliczyć kadłub wklęsły.

Farid Cheraghi
źródło
Twój pierwszy link wydaje się być uszkodzony.
PolyGeo