Monday, 3 February 2020

geometry - arcpy: get list of adjacent lines per line


I am working with polylines with hundreds of thousands of features, and would like to write a script that returns the OIDs of all lines adjacent (intersecting with) each line. I have been experimenting with the Geometry object and the shp.touches() method, but it seems as though a cursor with a nested cursor on the same feature is the only way to use it, and this takes much too long. The main goal is to perform an "Eliminate" task on a polyline, where short lines are merged to adjacent longer lines that share some unique ID. Collecting information in a dictionary has always been my technique of choice because of its speed, but dealing with geometries complicates this significantly. Unfortunately, Eliminate only works on polygons, and is quite slow itself. I have already tried to record shape geometries in a dictionary (dict[row.OBJECTID] = [row.Shape_Length, row.Shape, row.uniqueID]) and attempted a nested iteration, but the geometry objects don't seem to be interacting in the same way they do within a cursor iteration (shp.touches() doesn't work). Has anyone come across any techniques I could use for this problem?



Answer



I try to avoid answering my own questions, but I came up with an arcpy solution. Kudos to kenbuja, who provided a good ArcObjects solution, but the question asked for an arcpy solution. To address PolyGeo, this particular solution does not require creating a topology, but simply accessing the Shape object. A topology-based solution may also work, but I've had trouble automating the corrections within the topology framework. I would love to see a topology solution though. Here's my cleanLineGeom function, which assumes that adjacent lines share a startPoint / endPoint combination.


This function was written for stream networks, hence the params "streamID" and "segID", but can be used for any polyline. "StreamID" refers to the line ID, while "segID" refers to the ID for smaller lines that share a line ID. The function will only merge short (length < lineClusterTolerance) line segments to longer (length > lineClusterTolerance) line segments that share a "streamID". Lastly, the output is a new feature class, so the input is not altered.


Notice -- there will be some short lines left untouched in the process. These are lines that had no longer lines adjacent. These lines are exceptions, and can simply be dissolved normally after the function is run.


def cleanLineGeom(inLine,outName,streamID,segID,lineClusterTolerance):

lyrs = []
inLineName = arcpy.Describe(inLine).name
oidFieldName = arcpy.Describe(inLine).oidFieldName
# Separate short and long lines into different layers, then select all longs that touch shorts
shortLines = arcpy.MakeFeatureLayer_management(inLine,'shortLines','"SHAPE_LENGTH"<='+str(lineClusterTolerance))
lyrs.append(shortLines)
longLines = arcpy.MakeFeatureLayer_management(inLine,'longLines','"SHAPE_LENGTH">'+str(lineClusterTolerance))
lyrs.append(longLines)
arcpy.SelectLayerByLocation_management(longLines,"BOUNDARY_TOUCHES",shortLines,'',"NEW_SELECTION")


# Make a dictionary relating shortLine streamID/segID pairs to their origin- and endpoint coordinate pairs
shortDict = {}
rows = arcpy.SearchCursor(shortLines)
for row in rows:
shp = row.Shape
shortDict[(row.getValue(streamID),row.getValue(segID))] = [(shp.firstPoint.X,shp.firstPoint.Y),(shp.lastPoint.X,shp.lastPoint.Y)]
del rows

# Make a dictionary relating longLine origin- and endpoint coordinate pairs to segIDs
longDict = {}

rows = arcpy.SearchCursor(longLines)
for row in rows:
shp = row.Shape
firstCoords = (shp.firstPoint.X,shp.firstPoint.Y)
lastCoords = (shp.lastPoint.X,shp.lastPoint.Y)
longDict[firstCoords] = (row.getValue(streamID),row.getValue(segID))
longDict[lastCoords] = (row.getValue(streamID),row.getValue(segID))
del rows

# Create new dictionary relating shortLine segIDs to longLine segIDs that share a point

dissolveDict = {}
# If a shortLine's coordinate pair matches an entry in longDict,
# and the longLine's streamID matches, add their segIDs to dissolveDict
for ids, coordPairs in shortDict.iteritems():
for coords in [coordPairs[0],coordPairs[1]]:
if coords in longDict.iterkeys():
if longDict[coords][0] == ids[0]:
dissolveDict[ids[1]] = longDict[coords][1]

# Give all longLines a 'dissolve' value equal to their segID

arcpy.AddField_management(inLine,'dissolve','LONG')
arcpy.SelectLayerByAttribute_management(longLines,"CLEAR_SELECTION")
arcpy.CalculateField_management(longLines,'dissolve','[{0}]'.format(segID),'VB')

# If shortLine in dissolveDict, give it a 'dissolve' value equal to the dissolveDict value
# Else give it its own segID
urows = arcpy.UpdateCursor(shortLines,'','',segID+';dissolve')
for urow in urows:
if dissolveDict.get(urow.getValue(segID)):
urow.dissolve = dissolveDict[urow.getValue(segID)]

else:
urow.dissolve = urow.getValue(segID)
urows.updateRow(urow)
del urows

arcpy.Dissolve_management(inLine,inLineName+'D',"dissolve",'','SINGLE_PART')
cleaned = fieldFunc.fasterJoinAll(inLineName+'D',"dissolve",inLine,segID,outName)
arcpy.DeleteField_management(cleaned,'dissolve')
arcpy.CalculateField_management(cleaned,segID,'[{0}]'.format(oidFieldName),'VB')


#Clean-up
for lyr in lyrs:
arcpy.Delete_management(lyr)

return cleaned

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...