1

I'm trying to intersect two polyline feature classes and generate a point output using geometry objects. Its intersecting almost 2 million features with 250,000+ features, so running it with arcpy.Intersect_analysis has been running for...probably a week at this point with no progress. So I'm looking for a faster way. I'm getting AttributeError: 'list object has no attribute 'intersect', this occurs on line 13 at intersect = ln_geom.intersect(fn_geom, 1). All the documentation I've found says that is how to do it.

import os
import arcpy

Fishnet = r"C:\Users\AO_Fishnet_Testing"
Workspace = r"C:\Users\Testing.gdb"
Lines = r"C:\Users\Testing.gdb\AO_Testing_Ln"
Output = r"C:\Users\Testing.gdb\Testing_Points"

ln_geom = arcpy.CopyFeatures_management(Lines, arcpy.Geometry())
print "Lines Copied"
fn_geom = arcpy.CopyFeatures_management(Fishnet, arcpy.Geometry())
print "Fishnet Copied"
intersect = ln_geom.intersect(fn_geom, 1)
print "Intersected points"
arcpy.CopyFeatures_management(intersect, Output)
  • 1
    Please edit the question to show which line is generating that error. If you really need to intersect all features, doing so with a binary will be faster than by Python. You can try using an in_memory workspace, though you may be a bit large for that option. If you choose to use Python, you need to set up a nested loop using DA cursors -- this is not a linear operation, which is why your code fails. – Vince Jun 15 '16 at 14:37
  • Pretty sure Python will run out of memory if I try to do it in memory. You are saying to set up a cursor to loop through each feature, do the intersect, then add the geometry to another feature class? – BenStrike Jun 15 '16 at 14:45
  • I'm not an arc-gis user, but the Python error is telling you that there is no intersect method of the ln_geom object. I would inspect that object after it's creation to see what I've got. I've always found that the best way to do something like this is by using FME. Right tool for the right job, and so on. You can download the trial version and get the whole thing done, including the learning curve, in a couple of hours. – marcp Jun 15 '16 at 14:53
  • 1
    ln_geom is a list of polyline geometries, so you are calling a method "intersect" on a list which clearly makes no sense. @Vince way is best, use a cursor stepping through each feature, get its geometry and then call the intersect method as you have done and store the resulting point in another featureclass. – Hornbydd Jun 15 '16 at 14:56
  • 1
    "in_memory" is a special class of workspace. Even millions of features can work in_memory, provided that the attributes are not too wide. You should be careful with ordering of parameters, so that the most efficient comparisons are made (though some ArcGIS commands include feature count and spatial index optimization). – Vince Jun 15 '16 at 16:30
3

To do this efficiently you should read the 250,000 lines into a list with a da cursor (I assumed AO_Testing_Ln had the smaller number of features, but if AO_Fishnet_Testing has fewer lines it should be placed in the list instead). Then you should loop through a cursor on the Fishnet and loop through each line in the lineList to check for intersections. You should first test to determine if the extents intersect using a helper method I use then test for disjoint of the actual line shapes before intersecting, since both of those methods are faster than intersect for eliminating non-intersecting lines. I assumed you wanted to retain the OID of the Lines and Fishnet to get attributes from the lines later or to create a relate. The code below is untested, so errors may occur, but the code should be structured correctly:

import arcpy

Lines = r"C:\Users\Testing.gdb\AO_Testing_Ln"

lineList = [(r[0], r[1], r[1].extent) for r in arcpy.da.SearchCursor(Lines, ["OID@", "SHAPE@"])]

Fishnet = r"C:\Users\AO_Fishnet_Testing"
Output = r"C:\Users\Testing.gdb\Testing_Points"
env.workspace = r"C:\Users\Testing.gdb"

# Create a point feature class with the spatial reference of Lines
arcpy.CreateFeatureclass_management(r"C:\Users\Testing.gdb", "Testing_Points", "POINT", "", "DISABLED", "DISABLED", r"C:\Users\Testing.gdb\AO_Testing_Ln")

# Add a lineid field so that later attributes can be joined using Lines OIDs
arcpy.AddField_management("Testing_Points", "lineID", "LONG", 9, "", "", "lineID", "NULLABLE", "REQUIRED")

# Add a fishnetid field so that later attributes can be joined using Fishnet OIDs
arcpy.AddField_management("Testing_Points", "fishnetID", "LONG", 9, "", "", "fishnetID", "NULLABLE", "REQUIRED")

# create an insert cursor
fields = ['lineID', 'fishnetID', 'SHAPE@XY']
insertCursor = arcpy.da.InsertCursor(Output, fields)

counter = 0
# loop through the Fishnet lines and only intersect lines with extent overlap and not disjoint
with arcpy.da.SearchCursor(Fishnet, ["OID@", "SHAPE@"]) as cursor:
    for row in cursor:
        for i in range(0,len(lineList)):
            if rectIntersect(row[1].extent, lineList[i][2]):
                if not row[1].disjoint(lineList[i][1]):
                    pts = row[1].intersect(lineList[i][1], 1)
                        if pts.type == "point":
                            pts = [pts]
                        for p in pts:
                            insertCursor.insertRow((linelist[i][0], row[0],[p.X, p.Y]))
                            del p
                    del pts
        counter += 1
        if counter%10000 == 0: # adjust modulus amount to monitor progress
            print("Completed processing of {0} Fishnet Lines".format(counter))

del insertCursor

# Helper method for determining if two feature extents do or do not intersect
def rectIntersect(R1, R2):
    if ((R1.XMax < R2.XMin) or (R1.XMin > R2.XMax) or (R1.YMax < R2.YMin) or (R1.YMin > R2.YMax)):
        return False
    else:
        return True
  • I honestly don't need to do anything other than generating the intersections as a point, so I can split the lines on those. Basically what I'm doing is trying to count the number of lines within each grid of the fishnet and create a raster from it. But I want total features, not just features with their centroid in the cell, so I'm splitting the lines where they intersect the grid. I will play around with this though and let you know if I get it working. – BenStrike Jun 15 '16 at 17:15
  • The real problem with even my method is that if you process the entire set of features at once you have to make 500,000,000,000 comparisons. You would be better to take an encompassing polygon that is cut into 10 pieces then use Select By Location on the lines for each cut polygon, and process about 200,000 fishnets against 25,000 lines 10 separate times and merge the results. That reduces the number of feature comparisons to 50,000,000,000, which is 450,000,000,000 fewer comparisons. If you did 100 polygons the comparisons drop to about 5,000,000,000, or 495,000,000,000 fewer comparisons. – Richard Fairhurst Jun 15 '16 at 20:13
  • So if you want to take 10 to 100 times longer to get the same results, process everything all at once rather than 10 or 100 separate spatial selections on the two feature sets and running intersect the corresponding number of times (which could be scripted). – Richard Fairhurst Jun 15 '16 at 20:18
  • Will do, I'm trying to figure out what's the max number I can do before the search cursor gives me a memory error, which it does for around 700k lines. Then I'll just divide it up by that and process it that way. – BenStrike Jun 16 '16 at 12:46
  • I added a del p line and a del pts line to the code to see if that affects memory. Or is the memory error happening before it gets to the loop where the lines are intersected? Did the counter actually get to 700k? Or is this number of lines based on what you can load into the python list? Anyway, I would like to know more about the memory usage of the code to help isolate any possible memory leaks or develop coding techniques to manage memory with large data sets. – Richard Fairhurst Jun 16 '16 at 14:32

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.