Membuat poligon yang menghubungkan titik akhir dari beberapa baris menggunakan ArcPy?

10

Saya mencoba mencari cara untuk membuat poligon yang menghubungkan semua titik akhir dari sebuah shapefile yang berisi satu set poligon dengan skrip pythons di ArcGIS, saya mengalami kesulitan melakukan ini karena urutan node dalam poligon adalah penting. Saya ingin mencapai poligon abu-abu pada gambar dari garis hijau

Saya ingin menghubungkan titik akhir dari garis hijau untuk membuat poligon abu-abu tanpa harus melakukannya secara manual

Amanda
sumber
apakah baris Anda memiliki beberapa atribut untuk dipesan?
Ian Turton
pertama, Anda perlu urutan didefinisikan sebagai @iant bertanya, maka Anda memerlukan aturan apakah menghubungkan titik akhir ke titik awal berikutnya atau melakukannya dengan cara lain
Matej
3
gagal bahwa mungkin semacam hull alpha pada titik akhir?
Ian Turton
Garis ini pada tingkat tertentu memiliki atribut untuk memberi mereka perintah. Mereka memiliki nomor ID, tetapi untuk contoh di atas cabang kanan memiliki ID 1-7, 15-21 kiri dan setelah mereka terhubung ID adalah 22-27
Amanda
1
Anda dapat menjadi sangat dekat dengan a) membuat TIN, menggunakan garis, b) mengubah TIN menjadi segitiga c) memilih segitiga yang berbagi batas dengan garis. Anda hanya akan memiliki 1 poligon untuk dihapus di atas
FelixIP

Jawaban:

11

LANGKAH:

Hitung titik pusat bagian: masukkan deskripsi gambar di sini

Membangun pohon spanning minimum Euclidean mereka, melarutkannya dan menghitung buffer, jarak yang sama dengan setengah panjang bagian terpendek: masukkan deskripsi gambar di sini

Buat titik akhir bagian dan hitung rantainya (jarak sepanjang garis) pada batas buffer (versi buffer polyline tertutup): masukkan deskripsi gambar di sini

Urutkan titik akhir dalam urutan menaik menggunakan bidang rantai. Poin di bawah ini dilabeli oleh FID mereka:

masukkan deskripsi gambar di sini

Buat poligon dari set poin yang dipesan: masukkan deskripsi gambar di sini

Naskah:

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

Saya tahu ini sepeda, tapi itu milik saya dan saya menyukainya

FelixIP
sumber
2

Saya memposting solusi ini untuk QGIS di sini karena ini adalah perangkat lunak gratis dan mudah diterapkan. Saya dianggap hanya "cabang" yang tepat dari layer vektor polyline; karena dapat diamati pada gambar berikutnya (12 fitur di tabel atribut):

masukkan deskripsi gambar di sini

Kode (algoritme dalam pemahaman daftar python satu baris), untuk berjalan di Konsol Python QGIS, adalah:

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

Setelah menjalankan kode:

masukkan deskripsi gambar di sini

itu diproduksi lapisan memori poligon (dengan 11 fitur di tabel atributnya). Ini bekerja dengan baik.

xunilk
sumber
1

Anda dapat memilih titik akhir yang akan berpartisipasi dalam poligon, membuat TIN hanya dari titik-titik tersebut. Konversikan TIN menjadi poligon, larutkan poligon. Trik untuk mengotomatiskan proses ini adalah menentukan titik mana yang berkontribusi pada setiap poligon. Jika Anda memiliki garis dengan arah yang valid, dan garis-garis itu semua memiliki atribut yang sama, Anda bisa menulis kueri untuk mengekspor, katakan simpul akhir menggunakan simpul garis ke titik, lalu pilih dengan atribut titik-titik yang memiliki nilai atribut umum.
Lebih baik mengekstrak / memilih titik, membaca nilai x, y menggunakan kursor, menggunakan nilai x, y untuk menulis poligon baru. Saya tidak bisa melihat gambar terlampir di posting Anda, tetapi jika urutan poin penting maka setelah Anda memiliki nilai x, y disimpan dalam daftar Python, urutkan mereka. http://resources.arcgis.com/EN/HELP/MAIN/10.1/index.html#//002z0000001v000000

GBG
sumber
1

Memperluas pada komentar @iant, geometri terdekat dengan foto Anda adalah bentuk alpha (alpha hull) dari titik akhir. Untungnya banyak utas yang diterima dengan baik telah dijawab di GIS SE. Sebagai contoh:

Untuk mengatasi masalah Anda, pertama-tama gunakan Feature To Point untuk mengekstrak poin akhir. Kemudian gunakan alat python dari Tautan ini untuk menghitung lambung cekung.

Farid Cheraghi
sumber
Tautan pertama Anda tampaknya rusak.
PolyGeo