# Is there a faster way to see if a polyline intersects itself?

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)
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
trip_id = row
intCount = 0
lenSegs = len(row)
print lenSegs

for i in range(0, lenSegs - 1, 1):
P = segments[i]
Q = segments[i+1]
for j in range(1, lenSegs-1, 1):
R = segments[j]
S = segments[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

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: Create a new topology with rule "Must not self intersect": • 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