2

I have several large polyline feature classes and I need to find how many times each polyline intersects itself (where it intersects is of no importance). I wrote the following code, which works fine, the problem is it's been running for several hours and hasn't finished going through the first polyline. Basically, it takes each segment along the polyline and compares it with every other segment to see if they intersect. The first polyline in my sample data has 1488 vertices, (some have a lot more, some have fewer). That's 1487! comparisons for just one line. Does anyone know of a quicker way to do this? Some of my feature classes have over 3000 polylines!

import arcpy


filepath = r'C:\projects\TriCounty.gdb\Florida_GDL_Albers'

arcpy.env.workspace = filepath

These two functions determine if the line segments intersect by looking at orientation

def ccw(a, b, c):
    return (c.Y-a.Y)*(b.X-a.X) > (b.Y-a.Y)*(c.X-a.X)


def intersect(a, b, c, d):
    return ccw(a, c, d) != ccw(b, c, d) and ccw(a, b, c) != ccw(a, b, d)

Determine spatial reference of the input file

sr = arcpy.Describe("LinesTest").spatialReference

Create a new feature class to store info

arcpy.CreateFeatureclass_management(filepath,'Bikemap_SelfInt', "POLYLINE", spatial_reference=sr)
arcpy.AddField_management('Bikemap_SelfInt', 'intCount', "DOUBLE")
arcpy.AddField_management('Bikemap_SelfInt', 'Trip_ID', "DOUBLE")
print 'Table successfully created.'

Prepare insert cursor to insert results into the newly created feature class.

Fields = ['SHAPE@', 'intCount', 'Trip_ID']
insCur = arcpy.da.InsertCursor("Bikemap_SelfInt", Fields)

Search cursor to go through the polylines of the input file. I know there is only one set of points so I didn't have to iterate through parts.

with arcpy.da.SearchCursor("LinesTest", ["SHAPE@", "trip_id"]) as cursor:

    for row in cursor:
        segments = row[0]
        trip_id = row[1]
        intCount = 0
        lenSegs = len(row[0][0])
        print lenSegs

        for i in range(0, lenSegs - 1, 1):
            P = segments[0][i]
            Q = segments[0][i+1]
            for j in range(1, lenSegs-1, 1):
                R = segments[0][j]
                S = segments[0][j+1]

                if intersect(P, Q, R, S):
                    intCount += 1
                    print intCount
        insCur.insertRow([segments, intCount, trip_id])
del insCur

Like I said, this works great, but I don't have months (years?) to let it run! Any ideas on how to speed up this process?

  • You could toss the segments into an in_memory feature class, then iterate with a spatial search looking for interior intersect with id != self. Then you'd eliminate the O(N^3) algorithm. – Vince Mar 30 '20 at 12:14
  • FYI for anyone else with this or a similar question. BERA's method finds all the points where a polyline feature self-intersects. Using the "Export topology errors" tool in the data management toolbox, you can get a point (line, polygon) feature class of errors. In my case I have a trip_id I want associated with the points of self-intersection. I did this by joining my original feature layer to the newly created "error" point layer. – Mathwiz Mar 30 '20 at 12:31
3

Dont know what is fastest but I would try Topology:

Create a file geodatabase - feature dataset - Import your data into the feature dataset then Topology:

enter image description here

Create a new topology with rule "Must not self intersect":

enter image description here

Validate: enter image description here

  • I'll give it a go, thanks. – Mathwiz Mar 30 '20 at 8:27
  • I think this could work for me if there was a way to export that error inspector table (and get the number of errors per line in the feature class). Do you know if there is a way to display the results for only one line? For example, in my test file I have 14 polylines, and there are 4687 errors, but can I display the errors for an individual line and not all 14? – Mathwiz Mar 30 '20 at 9:07
  • Probably, I have not tried it. – BERA Mar 30 '20 at 12:11

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.