Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have donut buffers generated to protect anonymity of household locations. They need to be clipped so that they are only inside the census tract that the original point was in. The code takes about 8 hours on a very fast computer. Any optimizations to the following code would be greatly appreciated:

''' Clipbuff.py
'''
# Import arcpy module
import arcpy, sys, os, datetime

print datetime.datetime.now()
# Set up workspace and scratch workspace
workspace = "C:/Users/percy/projects/CEE/PerturbationProcessDataset/jitters.gdb/UTM/"
arcpy.env.workspace = workspace
scratch_workspace = "C:/users/percy/projects/CEE/PerturbationProcessDataset/scratch/"


# Allow overwriting
arcpy.env.overwriteOutput = True

# Local variables:
points = "OHAS_household"
id_Field = "OBJECTID"
polygons = "Geography_Tract"
buffers = "buff10"


# Open a search cursor on the buffers layer
rows = arcpy.SearchCursor(buffers)
print rows
# Create a for-loop to iterate through buffer IDs (via searchcursor),
# make feature layers of both the buffer and census tract layers,
# select containing census tract,
# clip by selected tract  
x = 1
for row in rows:
    myVal = row.getValue("SAMPN")
    mySelection = '"SAMPN" IN ('+ str(row.getValue("SAMPN")) +")"
    # Abstract the current row value
    print datetime.datetime.now()
    # Make a feature layer with the buffers layer, selecting the OBJECTID of the     current row/record
    arcpy.MakeFeatureLayer_management(buffers, "buffers_layer", mySelection)

    # Make a feature layer with the corresponding points layer
    arcpy.MakeFeatureLayer_management(points, "points_layer", mySelection)

    # Make a feature layer for the census tracts layer   
    arcpy.MakeFeatureLayer_management(polygons, "polygons_layer")

    # Select by location on the census tracts layer
    arcpy.SelectLayerByLocation_management("polygons_layer", "CONTAINS", "points_layer")

    # Clip the queried buffer by its spatially coincident census tract
    arcpy.Clip_analysis("buffers_layer", "polygons_layer", scratch_workspace + "buffclip" + str(x))
    # Add 1 to the file suffix
    x += 1
    # Delete the feature layers    
    arcpy.Delete_management("polygons_layer")
    arcpy.Delete_management("buffers_layer")
    arcpy.Delete_management("points_layer")

# Set the scratch workspace as the current workspace to list the clipped buffer featureclasses
arcpy.env.workspace = scratch_workspace

# Merge all clipped buffers into one layer
mylist = arcpy.ListFeatureClasses("buffclip*")
mytarget = mylist[0]

del mylist[0]

arcpy.Append_management(mylist, mytarget, "NO_TEST")
print datetime.datetime.now()

New version of code:

# Import arcpy module
import arcpy, sys, os, datetime

print datetime.datetime.now()
# Set up workspace and scratch workspace
workspace = "C:/projects/CEE/PerturbationProcessDataset/jitters.gdb/UTM/"
arcpy.env.workspace = workspace
scratch_workspace = "C:/projects/CEE/PerturbationProcessDataset/scratch1.gdb/"


# Allow overwriting
arcpy.env.overwriteOutput = True

# Local variables:
points = "OHAS_household_subset"
id_Field = "OBJECTID"
polygons = "Geography_Tract"
buffers = "buff1"

arcpy.MakeFeatureLayer_management("Geography_Tract", "polygons" )
arcpy.MakeFeatureLayer_management("buff1", "buffers")
arcpy.MakeFeatureLayer_management(points, "points")

# Open a search cursor on the buffers layer
rows = arcpy.SearchCursor(points)
x = 1
#for row in rows:
row = rows.next()
myVal = row.getValue("SAMPN")
mySelection = '"SAMPN" IN ('+ str(row.getValue("SAMPN")) +")"
print mySelection
#print row[0]
print datetime.datetime.now()
arcpy.SelectLayerByAttribute_management("points", "NEW_SELECTION",mySelection)
ptcount = arcpy.GetCount_management("points")
arcpy.SelectLayerByLocation_management("polygons", "CONTAINS", "points")
arcpy.SelectLayerByAttribute_management("buffers", "NEW_SELECTION",mySelection)
polycount = arcpy.GetCount_management("polygons")
buffcount = arcpy.GetCount_management("buffers")
print ptcount, polycount, buffcount
arcpy.Clip_analysis(polygons, buffers, scratch_workspace + "clippy" + str(x))
share|improve this question
    
wow, I have to figure out how to make my post look less ugly! :-) –  Percy Apr 26 '13 at 0:45
    
okay, that's better! –  Percy Apr 26 '13 at 5:36
    
by converting to file geodatabase for the scratch output, my estimate is about five hours for the 5000 buffer, but we need to do this for multiple inner and outer ring buffer combinations... Thanks for any insight! –  Percy Apr 26 '13 at 6:15
add comment

1 Answer

Almost any time you use cursors to manipulate geodatabase rows one by one, you should expect a slower process than if you can find a way to do the same thing but on larger selections. Compare the difference of, say, iterating through a point feature class containing 5000 points with a search cursor and buffering each one individually then merging the results compared to doing the same thing on the entire feature class in one go with a single command. Checking geometry (select by location) can also be slow, especially if you need to do a lot of it on complex geometry and ArcGIS adds a tons of vertices when it creates buffers - so even if the geometry looks simple to us (circular) it can still be effectively complex to an algorithm which just sees a ton of vertices).

So, I would look for ways to cut down the number of database calls and geometry checks. Here are some ideas I'd look at in your scenario:

  1. Do a spatial join of the census areas feature class on the buffers feature class to give the buffers the same id and the census track that contains their centroid. This will save time with repeatedly checking geometry with a 'select-by-location' call later. Selecting by an indexed attribute will usually be a lot quicker than testing topological relationships. In your current code, for every single buffer polygon, you must check ALL the census tracks' geometry. That's a lot of geometry testing! It's better to do it just once.
  2. Revise your code to iterate through the census tracks (of which I assume there are far fewer rows to iterate than the buffers) instead of the buffers (reduces the number of loops using a cursor).
  3. Select groups of buffers which share the same id (from the spatial join in step 1) as the currently selected census tract and clip them all in one go (batch processing of buffers and no geometry/topology testing).
  4. Rinse and repeat

This way you clip groups of buffers instead of doing it one by one, you avoid tons of effectively redundant topology checks and when you merge your files at the end, you have fewer of them to merge.

share|improve this answer
    
Thanks, that's a great idea, I was thinking along the same lines. The buffers already have the census tract id. I'm worried about the implementation of step 3, however... –  Percy Apr 27 '13 at 13:52
    
You are just clipping one selection with another. It should work just the same as if you were doing it from the UI. –  MappaGnosis Apr 28 '13 at 14:28
    
Ah, I see the problem I was worried about... I need the buffers to still retain one record per buffer. I think what you are suggesting in step 3 will result in one multipart buffer polygon per census tract. –  Percy Apr 29 '13 at 15:23
    
You fear is unfounded! You won't get a multipolygon. You will get what you want (one clipped poly per polygon in your input selection). It is just that doing the operation on a selection is quicker than doing the same thing using a cursor. –  MappaGnosis Apr 29 '13 at 18:17
    
You mean I will get back the same number of donuts that I started with? Cool, can't wait to try it! I'll let you know how it goes and mark your answer... Thanks again! –  Percy Apr 30 '13 at 21:19
show 2 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.