Создание полигонов, соединяющих конечные точки нескольких линий с помощью ArcPy?

10

Я пытаюсь выяснить, как создать многоугольник, который соединяет все конечные точки шейп-файла, содержащего набор полилиний с pythonscript в ArcGIS, у меня возникают проблемы с этим, так как порядок узлов в многоугольнике важен. Я хочу добиться серого многоугольника на картинке из зеленых линий

Я хочу соединить конечные точки зеленых линий, чтобы создать серый многоугольник без необходимости делать это вручную

Аманда
источник
у ваших строк есть какой-то атрибут, чтобы дать порядок?
Ян Тертон
во-первых, вам нужен порядок, определенный как @iant ask, затем вам нужно правило, соединяющее конечную точку со следующей начальной точкой или выполняющее ее любым другим способом
Matej
3
если не получится, может быть, какой-то альфа-корпус на конечных точках?
Ян Тертон
Строка в некоторой степени имеет атрибуты, чтобы дать им порядок. У них есть идентификационный номер, но для приведенного выше примера правая ветвь имеет ID 1-7, левая 15-21, а после подключения ID 22-27
Аманда
1
Вы можете подойти очень близко, а) создав TIN, используя линии, б) преобразовав TIN в треугольники в) выбрав треугольники, которые разделяют границу с линиями. У вас будет только 1 полигон для удаления вверху
FelixIP

Ответы:

11

ШАГИ:

Вычислить центральные точки секций: введите описание изображения здесь

Построить их евклидово минимальное остовное дерево, распустить его и вычислить буфер, расстояние, равное половине длины кратчайшего участка: введите описание изображения здесь

Создайте конечные точки сечения и вычислите их цепочку (расстояние вдоль линии) на границе буфера (закрытая полилиния, версия буфера): введите описание изображения здесь

Сортируйте конечные точки в порядке возрастания, используя поле цепочки. Точки ниже помечены их FID:

введите описание изображения здесь

Создайте многоугольник из упорядоченного набора точек: введите описание изображения здесь

Автор сценария:

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

Я знаю, что это велосипед, но он мой, и мне это нравится

FelixIP
источник
2

Я публикую это решение для QGIS здесь, потому что оно бесплатное и легко внедряемое. Я рассмотрел только правильную «ветвь» векторного слоя полилинии; как это можно увидеть на следующем изображении (12 объектов в таблице атрибутов):

введите описание изображения здесь

Код (алгоритм в понимании однострочного списка Python) для запуска на консоли Python QGIS:

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

После запуска кода:

введите описание изображения здесь

был создан слой памяти многоугольника (с 11 объектами в таблице атрибутов). Это работает хорошо.

xunilk
источник
1

Вы можете выбрать конечные точки, которые будут участвовать в многоугольнике, создать TIN только из этих точек. Конвертировать TIN в полигоны, растворить полигоны. Хитрость в автоматизации этого процесса состоит в том, чтобы решить, какие точки внести в каждый полигон. Если у вас есть линии с допустимыми направлениями, и у всех этих линий есть общий атрибут, вы можете написать запрос на экспорт, скажем, конечные вершины, используя вершины линий для точек, а затем выбрать по атрибуту те точки, которые имеют значение общего атрибута.
Лучше было бы извлечь / выбрать точки, прочитать значения x, y с помощью курсора, использовать значения x, y, чтобы написать новый многоугольник. Я не вижу прикрепленное изображение в вашем посте, но если порядок точек имеет значение, тогда, когда у вас есть значения x, y, сохраненные в списке Python, сортируйте их. http://resources.arcgis.com/EN/HELP/MAIN/10.1/index.html#//002z0000001v000000

GBG
источник
1

Если развернуть комментарий @iant, ближайшая к вашему снимку геометрия - это альфа-форма (альфа-оболочка) конечных точек. К счастью, на многие хорошо полученные темы уже ответили на GIS SE. Например:

Чтобы решить вашу проблему, сначала используйте Feature To Point, чтобы извлечь конечные точки. Затем используйте инструмент Python из этой ссылки, чтобы рассчитать вогнутый корпус.

Фарид Чераги
источник
Ваша первая ссылка, кажется, не работает.
PolyGeo