7

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

  • 1
    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. – Michael Stimson Apr 22 '14 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 '14 at 2:08
14
+50

Here is the 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: Just pipeline

The results: Pipelines and too close

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": enter image description here

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.

6

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
  * @copyright  : Licensed under a Creative Commons Attribution-Share Alike 2.5 Australia License. (http://creativecommons.org/licenses/by-sa/2.5/au/)
*/
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.

5

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. – user30184 May 30 '14 at 20:29
3

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.

http://resources.arcgis.com/en/help/main/10.2/index.html#//000800000019000000 http://resources.arcgis.com/en/help/main/10.2/index.html#//000800000044000000

0

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

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.