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 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?