Wednesday, 13 September 2017

Creating polygon connecting endpoints of multiple lines using ArcPy?


I'm trying to figure out how to create a polygon that connects all the endpoints of a shapefile containing a set of polyilnes with pythonscript in ArcGIS, I'm having trouble doing this as the order of the nodes in the polygon is important. I want to achieve the grey polygon in the picture from the green lines


I want to connect the endpoints of the green lines to create the grey polygon without having to do it manually



Answer




STEPS:


Compute sections centre points: enter image description here


Built their Euclidean minimum spanning tree, dissolve it and compute buffer, distance equal half of shortest section length: enter image description here


Create section end points and compute their chainage (distance along line) on the boundary of the buffer (closed polyline version of buffer): enter image description here


Sort end points in ascending order using chainage field. Points below labelled by their FID:


enter image description here


Create polygon from ordered set of points: enter image description here


Script:


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

I know it is a bicycle, but it’s my own and I like it


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...