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
STEPS:
Compute sections centre points:
Built their Euclidean minimum spanning tree, dissolve it and compute buffer, distance equal half of shortest section length:
Create section end points and compute their chainage (distance along line) on the boundary of the buffer (closed polyline version of buffer):
Sort end points in ascending order using chainage field. Points below labelled by their FID:
Create polygon from ordered set of points:
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
I post this solution for QGIS here because it is free software and easy to implement. I considered only the right "branch" of polyline vector layer; as it can be observed at the next image (12 features at attributes table):
The code (algorithm in a one line python list comprehension), for running at the Python Console of QGIS, is:
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()
After running the code:
it was produced a polygon memory layer (with 11 features at its attributes table). It works nicely.
You could select the endpoints that will participate in a polygon, create a TIN from only those points. Convert the TIN to polygons, dissolve the polygons.
The trick to automating this process is deciding which points to contribute to each polygon. If you have lines with valid directions, and those lines all share some common attribute you could write a query to export say the end vertices using line vertices to points, then select by attribute those points that have the common attribute value.
Better would be to extract/select the points, read the x , y values using a cursor, use the x, y values to write a new polygon.
I cannot see an attached picture in your post but if point order matters then once you have the x, y values stored in a Python list, sort them.
http://resources.arcgis.com/EN/HELP/MAIN/10.1/index.html#//002z0000001v000000
Expanding on @iant comment, the closest geometry to your snapshot is the alpha shape (alpha hull) of the endpoints. Fortunately many well received threads have been already answered on GIS SE. For example:
Create Complex Polygon from Point Layer using only Boundary Points in ArcGIS
What are Definition, Algorithms and Practical Solutions for Concave Hull?
To solve your problem, first use Feature To Point to extract the end points. Then use the python tool from this Link to calculate the concave hull.
