Finding close and parallel lines using ArcPy?

I need help with script for QA of water network. I need to find all lines that are almost parallel and close to each other less then 2 meters.

The script has 3 stages:

1. add azimuth field and calculate azimuth.



2. iterate and find close lines.


3. If lines have almost the same azimuth - select them.

I've used these lines here for first stage and need help with the loop that will go thru all the line and find close pairs.

Answer

Bearing does not apply to a whole line, rather to a segment (2-point line), so the first step is to break down the polylines to two point lines. With the segments apply the angle, but simplified to only the first two quadrants of the circle (0-180 degrees) by flipping the lines if the Y coordinate is lower at the end - lines go both ways. Then generate a near table, compare the angles, join and export.

I started with some data I had lying around:

The results:

The too close and parallel lines are red. Unmarked but parallel lines are outside of the distance (I used 2 metres) they are about ten metres apart.

It also finds "dog legs":

There is no reason for those two vertices to be there.

Now, the code!

import sys, os, arcpy, math, time

InFeatureClass = sys.argv[1] # Input, feature class
SearchDistance = sys.argv[2] # Input, distance unit or number
OutFeatureClass = sys.argv[3] # Output, feature class

TempDir = os.environ.get("Temp")

# assumed values
DBname = "Temp_" + time.strftime("%Y-%m-%d_%H%M%S") + ".gdb"
AngleTolerance = 10
FlipAngle = 180 - AngleTolerance

TempDB = TempDir + "\\" + DBname
SegmentsFC = TempDB + "\\Segments"
NearTable = TempDB + "\\NearTable"
NearDist = TempDB + "\\NearDist"
NearAngle = TempDB + "\\NearDistAngle"


def GetAngle(FromX,FromY,ToX,ToY):
    # Calculate the angle of the segment
    # simplified in the 0 to pi range
    if FromY < ToY:
        fX = FromX
        fY = FromY
        tX = ToX
        tY = ToY
    else:

        fX = ToX
        fY = ToY
        tX = FromX
        tY = FromY
    dX = tX - fX
    dY = tY - fY
    RadAngle = math.atan2(dY,dX)
    DegAngle = (RadAngle * 180) / math.pi
    return DegAngle # In degrees, would work in radians but the numbers are easier


# Get some important information about the input
desc = arcpy.Describe(InFeatureClass)
SR = desc.spatialReference
ShapeName = desc.shapeFieldName

# create temp database and feature class
arcpy.AddMessage("Creating temp database and feature class")
arcpy.CreateFileGDB_management(TempDir,DBname)
arcpy.CreateFeatureclass_management(TempDB,"Segments","POLYLINE",None,"DISABLED","DISABLED",SR)
arcpy.AddField_management(SegmentsFC,"Angle","DOUBLE")

# break the input lines into segments
sCur = arcpy.SearchCursor(InFeatureClass)
iCur = arcpy.InsertCursor(SegmentsFC,SR)
lArray = arcpy.Array()

arcpy.AddMessage("Breaking lines into segments")

for fromRow in sCur:
    geom = fromRow.getValue(ShapeName)
    partNum = 0

    for part in geom:
        thisPart = geom.getPart(partNum)
        firstPoint = True
        for pnt in thisPart:
            if pnt != None:
                if firstPoint:
                    # Skipping the first point, a segment needs two
                    # points therefore segments = points - 1
                    firstPoint = False
                    prePoint = pnt

                else:
                    # use the previous point and create a segment
                    lArray.add(prePoint)
                    lArray.add(pnt)

                    # Create a new row buffer and set shape
                    feat = iCur.newRow()
                    feat.shape = lArray

                    # Set the angle of the segment

                    feat.setValue("Angle",GetAngle(prePoint.X,prePoint.Y,pnt.X,pnt.Y))

                    # insert the new feature and clear the array
                    iCur.insertRow(feat)
                    lArray.removeAll()

                    prePoint = pnt # save the pnt for the next segment
        partNum += 1
del sCur
del iCur


# Generate near table of ALL features within search distance
arcpy.AddMessage("Generating near table")
arcpy.GenerateNearTable_analysis(SegmentsFC,SegmentsFC,NearTable,SearchDistance,"NO_LOCATION","NO_ANGLE","ALL")
# reduce the near table to just the non-touching features
arcpy.TableSelect_analysis(NearTable,NearDist,"NEAR_DIST > 0")
# add fields for from feature angle, to feature angle
arcpy.AddField_management(NearDist,"FAngle","DOUBLE")
arcpy.AddField_management(NearDist,"TAngle","DOUBLE")
arcpy.AddField_management(NearDist,"AngDiff","DOUBLE")


# create a join to copy the angles to the from and to angle fields
arcpy.AddMessage("Copying angles")
arcpy.MakeTableView_management(NearDist,"ND")
arcpy.AddJoin_management("ND","IN_FID",SegmentsFC,"OBJECTID")
arcpy.CalculateField_management("ND","NearDist.FAngle","!Segments.Angle!","PYTHON")
arcpy.RemoveJoin_management("ND")

arcpy.AddJoin_management("ND","NEAR_FID",SegmentsFC,"OBJECTID")
arcpy.CalculateField_management("ND","NearDist.TAngle","!Segments.Angle!","PYTHON")

arcpy.RemoveJoin_management("ND")

# calculate the difference in angle
arcpy.AddMessage("Resolving differences in angle")
arcpy.CalculateField_management(NearDist,"AngDiff","abs(!FAngle! - !TAngle!)","PYTHON")
# flip where one is near 180 and the other is near 0
arcpy.MakeTableView_management(NearDist,"NDA","AngDiff > %s" % FlipAngle) 
arcpy.CalculateField_management("NDA","AngDiff","180 - !TAngle!","PYTHON")

# Reduce the near table to similar angles

arcpy.TableSelect_analysis(NearDist,NearAngle,"AngDiff < %s" % str(AngleTolerance))

# join to the table and export to OutFeatureClass
arcpy.AddMessage("Exporting records")
arcpy.MakeFeatureLayer_management(SegmentsFC,"SegFC")
arcpy.AddJoin_management("SegFC","OBJECTID",NearAngle,"IN_FID","KEEP_COMMON")
arcpy.CopyFeatures_management("SegFC",OutFeatureClass)

arcpy.AddMessage("Cleaning up")
arcpy.Delete(TempDB)

This code can be copied and put directly into a python tool and run from a toolbox.