# Finding close and parallel lines using ArcPy? [closed]

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.

• Finding close lines is a bit of a pain, you can find distance from a point though. Look at "Feature Vertices to Point" help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… (both ends and mid) and "Point Distance" help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… which will give you the distance of a point to ALL other lines in tolerance, just get rid of the ones that match by FID/OID. When comparing your similar azimuth don't forget that 2pi is close to 0. Apr 22, 2014 at 21:30
• Would you be able to edit your Question to mention the version of ArcGIS for Desktop that you are targeting, please? Also, so that your Question can standalone, I think you should edit in the lines of existing code that you are currently just linking to.
– PolyGeo
Jun 3, 2014 at 2:08

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 # Input, feature class
SearchDistance = sys.argv # Input, distance unit or number
OutFeatureClass = sys.argv # 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
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)
# break the input lines into segments
sCur = arcpy.SearchCursor(InFeatureClass)
iCur = arcpy.InsertCursor(SegmentsFC,SR)
lArray = arcpy.Array()

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

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

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

arcpy.CalculateField_management("ND","NearDist.TAngle","!Segments.Angle!","PYTHON")
arcpy.RemoveJoin_management("ND")

# calculate the difference 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.MakeFeatureLayer_management(SegmentsFC,"SegFC")
arcpy.CopyFeatures_management("SegFC",OutFeatureClass)

arcpy.Delete(TempDB)


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

I realize that you are not asking for an SQL solution for your problem, but consider for a second how quickly this problem could be solved using SQL. ArcGIS has the capability to integrate with Microsoft SQL server, and there is nothing saying you can't move the whole procedure to PostGIS. The following code will absolutely solve your problem, given that the links have unique start and end node IDs, and that you have set appropriate spatial and identity indexes. This solution relies on a function called Bearing, detailed here:

/**
* @function   : Bearing
* @precis     : Returns a bearing between two point coordinates
* @version    : 1.0
* @usage      : FUNCTION Bearing(@p_dE1  float,
*                                @p_dN1 float,
*                                @p_dE2 float,
*                                @p_dN2 float )
*                RETURNS GEOMETRY
*               eg select cogo.Bearing(0,0,45,45) * (180/PI()) as Bearing;
* @param      : p_dE1     : X Ordinate of start point of bearing
* @paramtype  : p_dE1     : FLOAT
* @param      : p_dN1     : Y Ordinate of start point of bearing
* @paramtype  : p_dN1     : FLOAT
* @param      : p_dE2     : X Ordinate of end point of bearing
* @paramtype  : p_dE2     : FLOAT
* @param      : p_dN2     : Y Ordinate of end point of bearing
* @paramtype  : p_dN2     : FLOAT
* @return     : bearing   : Bearing between point 1 and 2 from 0-360 (in radians)
* @rtnType    : bearing   : Float
* @note       : Does not throw exceptions
* @note       : Assumes planar projection eg UTM.
* @history    : Simon Greener  - Feb 2005 - Original coding.
* @history    : Simon Greener  - May 2011 - Converted to SQL Server
*/
Create Function [GIS].[Bearing](@p_dE1 Float, @p_dN1 Float,
@p_dE2 Float, @p_dN2 Float)
Returns Float
AS
Begin
Declare
@dBearing Float,
@dEast    Float,
@dNorth   Float;
BEGIN
If (@p_dE1 IS NULL OR
@p_dN1 IS NULL OR
@p_dE2 IS NULL OR
@p_dE1 IS NULL )
Return NULL;

If ( (@p_dE1 = @p_dE2) AND
(@p_dN1 = @p_dN2) )
Return NULL;

SET @dEast  = @p_dE2 - @p_dE1;
SET @dNorth = @p_dN2 - @p_dN1;
If ( @dEast = 0 )
Begin
If ( @dNorth < 0 )
SET @dBearing = PI();
Else
SET @dBearing = 0;
End
Else
SET @dBearing = -aTan(@dNorth / @dEast) + PI() / 2.0;

If ( @dEast < 0 )
SET @dBearing = @dBearing + PI();

Return @dBearing;
End
End;

GO


You would then call on this select procedure (assuming you want objects within 2 meters/feet and consider less than 15 degrees to be almost paralell):

      SELECT AKey, BKey
FROM
(
SELECT  *,CASE
WHEN ABS(ABearing - BBearing) <= 90  THEN ABS(ABearing - BBearing)
WHEN ABS(ABearing - BBearing) <= 180 THEN 180-ABS(ABearing - BBearing)
WHEN ABS(ABearing - BBearing) <= 270 THEN ABS(ABearing - BBearing)-180
WHEN ABS(ABearing - BBearing) <= 360 THEN 360-ABS(ABearing - BBearing)
END AS AngleDiff
FROM
(
SELECT  A.ID AS AKey,
B.ID AS BKey,
A.Bearing AS ABearing,
B.Bearing AS BBearing
FROM    (
SELECT  ID,
us_node_id,
ds_node_id,
GIS.Bearing(
Shape.STStartPoint().STX
,Shape.STStartPoint().STY
,Shape.STEndPoint().STX
,Shape.STEndPoint().STY
)* 180/PI()  as Bearing,
Shape
FROM    YourTableOfPipes
) AS A
INNER JOIN
(
SELECT  ID,
us_node_id,
ds_node_id,
GIS.Bearing(
Shape.STStartPoint().STX
,Shape.STStartPoint().STY
,Shape.STEndPoint().STX
,Shape.STEndPoint().STY
)* 180/PI() as Bearing,
Shape
FROM    YourTableOfPipes
) AS B
ON  A.us_node_id <> B.us_node_id
AND
A.ds_node_id <> B.ds_node_id
AND
A.us_node_id <> B.ds_node_id
AND
A.ds_node_id <> B.us_node_id
AND
A.Shape.STDistance(B.Shape) <=2
) AS Angle
) AS XX
WHERE AngleDiff <= 15


For a few thousand objects, this should take seconds. A couple million objects should take a few minutes. Should you not have a us/ds identifier, simply match on A.ID <> B.ID. You would need the us/ds identifiers to avoid matching conduits that are upstream/downstream of one another (spatially connected), but I get the feeling this sort of match is something that would not be target activity.

Your final objective seems to be the matching of line segments. We've found that open jump with roadmatcher was very usefull for this kind of process, and it is open source.

Within ArcGIS, the tool "collapse dual line to centerline" could also help you to identify those lines (there is a minimum and maximum width parameter).

• The topology plugin in OpenJUMP Plus has also a tool for finding almost same segments. Tools takes distance of the lines and angle between the lines as parameters. It could be worth trying. May 30, 2014 at 20:29

Perhaps you could create a 2m buffer with the dissolve option set to NONE around the lines and find the buffers inside of which more than one line falls. The arcpy tools I am thinking could be used for this are:

arcpy.Buffer_analysis("lines.shp","lines_buffer2m.shp","2 Meters","FULL","ROUND","NONE")
arcpy.TabulateIntersection_analysis("lines_buffer2m.shp","uniqueID","lines.shp","out_table","","Count")


You would need to add a field with unique IDs in the buffer polygons. You would also need to create a field in the lines feature class that has a value of 1 for all features (the field named "Count" in my example.) The buffer act as zones, and the Tabulate Intersection tool will sum on the Count field for every line in each buffer zone.

This is just an idea. It worked on a small sample that I just tried, but it may not work if your water network. It won't really work on multi-part lines or lines that are completely connected. I had to run Planarize lines first to make sure lines were broken at intersections.

I know this post is old, but I think that the best tool for doing what the original post described is the Detect Feature Changes tool that was first released with ArcGIS 10.2. At ArcGIS 10.4 it also detects line direction changes (flipped lines) that are otherwise identical. The tool optionally creates a table to identify 1:1 segment matches or 1:M, M:1 or M:N matches, so that you can quickly identify the set of lines that may need to split or be merged. For pictures of a practical use of the tool see my response to Selecting parallel line segments but not segments that intersect