arcpy: get list of adjacent lines per line

Question

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 1

I'm working an a project that does something similar, but I'm using ArcObjects. Its goal is to aggregate polygons that are smaller than a user-designated area with their larger adjacent neighboring polygons.

My logic was to find the smallest polygon and merging it with a larger, adjacent polygon. I chose this method instead of just cycling straight through a feature cursor, since after a polygon is merged with another one, it could still be smaller than the minimum area. This ensures that the resultant polygon can still be returned by the search for the smallest polygon. The SQL code to find the smallest polygon is for the layer "layername" (assuming it's in a geodatabase)

"Shape_Area" = (Select Min("Shape_Area") from layername

The adjacent polygon that it gets merged with is not chosen randomly, but has the closest value for a selected numeric attribute. You could modify this to only merge with a feature that has the same unique attribute.

The only little bug I'm still working on is when there is an island polygon. Currently, the code will exit since otherwise it would just get stuck in an endless loop, with the query filter always returning that island polygon as the smallest feature.

Here's the code for that section

pQFCounter.WhereClause = Chr(34) & "Shape_Area" & Chr(34) & " < " & CDbl(txtMMU.Text)
If pScratchFeatureClass.FeatureCount(Nothing) = 0 Then Exit Sub

pFLayer.FeatureClass = pScratchFeatureClass
pQFilter = New ESRI.ArcGIS.Geodatabase.QueryFilter
pQFilter.WhereClause = Chr(34) & "Shape_Area" & Chr(34) & " = (Select Min(" & Chr(34) & "Shape_Area" & Chr(34) & ") from " & pFLayer.Name & ")"
pFCursor = pScratchFeatureClass.Search(pQFilter, False)
pFeature = pFCursor.NextFeature

'I'm only going to use the first feature returned, since there may be multiple features that have the same minimum area
Do Until pFeature Is Nothing
    pArea = pFeature.Shape

    'avoiding the island polygon problem.
    If pFeature.OID = OldOID And pArea.Area = OldArea Then
        Debug.Print(pFeature.OID)
        Exit Do
    End If
    OldOID = pFeature.OID
    OldArea = pArea.Area

    pSFilter = New ESRI.ArcGIS.Geodatabase.SpatialFilter
    pSFilter.Geometry = pFeature.Shape
    pSFilter.GeometryField = pFClass.ShapeFieldName
    pSFilter.SpatialRel = ESRI.ArcGIS.Geodatabase.esriSpatialRelEnum.esriSpatialRelTouches

    pSelFCursor = pScratchFeatureClass.Search(pSFilter, False)
    pSelFeature = pSelFCursor.NextFeature

    'Finding the feature that has the closest attribute
    If pSelFeature IsNot Nothing Then
        pMergeFeature = pSelFeature
        MinDiff = Math.Abs(pFeature.Value(FieldIndex) - pSelFeature.Value(FieldIndex))

        pSelFeature = pSelFCursor.NextFeature

        Do Until pSelFeature Is Nothing
            Diff = Math.Abs(pFeature.Value(FieldIndex) - pSelFeature.Value(FieldIndex))
            If MinDiff > Diff Then
                pMergeFeature = pSelFeature
                MinDiff = Diff
            End If
            pSelFeature = pSelFCursor.NextFeature
        Loop

        pSelArea = pMergeFeature.Shape
        pTopoOp = CType(pMergeFeature.ShapeCopy, ESRI.ArcGIS.Geometry.ITopologicalOperator4)
        pTopoOp.IsKnownSimple_2 = False
        pTopoOp.Simplify()
        pMergeFeature.Shape = pTopoOp.Union(pFeature.ShapeCopy)
        pTopoOp = CType(pMergeFeature.ShapeCopy, ESRI.ArcGIS.Geometry.ITopologicalOperator4)
        pTopoOp.IsKnownSimple_2 = False
        pTopoOp.Simplify()

        pFeature.Delete()
        pMergeFeature.Store()
    Else
        Debug.Print("Feature " & pFeature.OID & " appears to be an island feature.")
        If Len(Excluded) > 0 Then
            Excluded += ", " & pFeature.OID
        Else
            Excluded = pFeature.OID
        End If
    End If


    System.Runtime.InteropServices.Marshal.FinalReleaseComObject(pSelFCursor)

    If Len(Excluded) > 0 Then
       'this is not giving me the expected cursor that would exclude the island polygons
        pQFSelection = New ESRI.ArcGIS.Geodatabase.QueryFilter
        pQFSelection.WhereClause = pScratchFeatureClass.OIDFieldName & " not in (" & Excluded & ")"
        pSelectionSet = pScratchFeatureClass.Select(pQFSelection, ESRI.ArcGIS.Geodatabase.esriSelectionType.esriSelectionTypeHybrid, ESRI.ArcGIS.Geodatabase.esriSelectionOption.esriSelectionOptionNormal, pWorkSpace)
        pSelectionSet.Search(pQFilter, False, pFCursor)
    Else
        pFCursor = pScratchFeatureClass.Search(pQFilter, False)
    End If
    pFeature = pFCursor.NextFeature
    pArea = pFeature.Shape

    If pArea.Area > CDbl(txtMMU.Text) Then Exit Do

Loop

Part of the code I didn't include copies the original feature class into a geodatabase, ensuring that the Shape_Area field exists and is automatically calculated.

Answer 2

I worked in a project something similar of this, but I used ArcObjects. My goal was connecting to adjacent polylines if one of it's end point is another one's starting point to make two short polylines a single polyline. My process was:

1. Dictionary<PointKey, FeatureDataList> polylineDictionary;

• PointKey is a class that contains a point.
• FeatureDataList is a class that contains List of IFeatures.

• Both classes overrides "Equals" and "GetHashCode" methods.

  Dictionary<PointKey, FeatureDataList> ToPointDictionary;
  Dictionary<PointKey, FeatureDataList> FromPointDictionary;
  public void CreateDictionary(IFeatureLayer featureLayer)
  {
      var featureFunctionality = new FeatureFunctionality();
      List<IFeature> features = GetAllFeatures(featureLayer.FeatureClass);
  
      foreach (var feature in features)
      {
          IPolyline polyline = GetPolylineFromFeature(feature);
          AddFeatureInDictionary(ToPointDictionary, feature, polyline.ToPoint);
          AddFeatureInDictionary(FromPointDictionary, feature, polyline.FromPoint);
      }
  }
  
  
  
  void AddFeatureInDictionary(Dictionary<PointKey, FeatureDataList> polylineDictionary, IFeature feature, IPoint point)
  {
      FeatureDataList featureDataList;
      PointKey key = PointKey.GetKey(point);
      if (!polylineDictionary.ContainsKey(key))
      {
          featureDataList = new FeatureDataList();
          featureDataList.Add(feature);
          polylineDictionary.Add(key, featureDataList);
      }
      else
      {
          featureDataList = polylineDictionary[key];
          featureDataList.Add(feature);
      }
  }
  

By these processes I made two dictionaries. After creating dictionaries, I check for if both dictionaries contain same point and in both dictionaries, that key has only one feature in feature list, then I created a new polyline with those two polylines and delete two short polylines.

Answer 3

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