9

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

  • do your lines have some attribute to give the ordering? – Ian Turton Oct 13 '15 at 12:17
  • first, you need the ordering defined as @iant asked, then you need the rule whether connecting endpoint to the next startpoint or doing it any other way – Matej Oct 13 '15 at 13:01
  • 3
    failing that maybe some sort of alpha hull on the end points? – Ian Turton Oct 13 '15 at 13:12
  • The line does to some degree have attributes to give them order. They have a ID number, but for the example above the rightbranch have ID 1-7, the left 15- 21 and after they're connected the IDs are 22-27 – Amanda Oct 13 '15 at 14:11
  • 1
    You can get very close by a) creating TIN, using lines, b) converting TIN to triangles c) selecting triangles that share a boundary with lines. You'll have only 1 polygon to delete at the top – FelixIP Oct 13 '15 at 22:27
9

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

3

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

enter image description here

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:

enter image description here

it was produced a polygon memory layer (with 11 features at its attributes table). It works nicely.

1

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

1

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:

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.