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 a feature class with 790,000 points (and it is constantly increasing). I need to select some of them (approx. 350,000), and then buffer that selection.

Currently, the process runs like this (pseudocode):

arcpy.MakeFeatureLayer_management(fcIn, tempLayer, '"CUST" LIKE \'%G%\'')
arcpy.Buffer_analysis(tempLayer, fcOut, "150 Feet", "FULL", "ROUND", "ALL", "")

The Buffer step is prone to failure due to memory errors, which isn't terribly surprising. To get around that, I would like to further subset the selection, creating a bunch of buffers and then merging them all at the end. I know that this will take at least as long as (and probably longer than) a single Buffer operation, but avoiding a complete failure is a higher priority.

The question: how do I structure a loop that iterates selecting 10,000 (or so) points at a time and then buffers them? I am struggling to turn the iteration concept into code that will work with feature layers.

(And of course, if this is a terrible, unworkable idea, please let me know. I'm not very committed to it yet.)

share|improve this question
    
You could create a new layer that contains only the points you want to buffer. Then select points based on FID and run your buffer on the selection. Deselect the points and move on to the next group. –  khafen Sep 3 at 18:39
    
Does the MakeFeatureLayer function generate new FIDs when it subsets? –  Erica Sep 3 at 18:46
    
@Erica no it does not –  Jason Scheirer Sep 3 at 19:34
    
Will buffers overlap, and if so do the buffers need to be individual or can they be dissolved? I note that buffer uses tiling when processing but only if using dissolve. –  Chris W Sep 3 at 19:37
1  
That seems strange to me. If you're using dissolve, it should be tiling, and so I wouldn't expect it to run out of memory. resources.arcgis.com/en/help/main/10.1/index.html#//… –  Chris W Sep 3 at 20:15

3 Answers 3

up vote 3 down vote accepted

As a shapefile FIDs' are contiguous and 0 based, you can use that to your advantage:

import sys, os, arcpy

InFC = sys.argv[1] # must be a shape file
OutFC = sys.argv[2]

# change as appropriate
DivisorSize = 1000
BufferDistance = 100
TempDir = os.environ.get("Temp")

MaxFeat = int(arcpy.GetCount_management(InFC).getOutput(0))
StepRange = range(0,MaxFeat,DivisorSize) # range of FIDs to process
Steps = range(len(StepRange))            # the steps to iterate over

# add the last FID if not an even number
if MaxFeat not in StepRange:
    StepRange.append(MaxFeat)
else:
    # shorten the steps
    Steps = range(len(StepRange)-1)            # the steps to iterate over
# steps must be one less than the length of StepRange
# as the values are StepRange[n] to StepRange[n + 1]

TempFCs = list() # keep track of the temp features created

for ThisStep in Steps:
    arcpy.AddMessage("Buffer %d of %d" % (ThisStep,len(Steps))) # 'still alive' message

    # create a layer with just the features in range
    DefQ = "FID >= %d and FID < %d" % (StepRange[ThisStep],StepRange[ThisStep + 1])
    arcpy.MakeFeatureLayer_management(InFC,"Layer",DefQ)

    # buffer to a temp feature class
    ThisTemp = "%s\\bTemp_%d_to_%d.shp" % (TempDir,StepRange[ThisStep],StepRange[ThisStep + 1])
    arcpy.Buffer_analysis("Layer",ThisTemp,BufferDistance,dissolve_option="ALL")

    TempFCs.append(ThisTemp)         # add to the list to keep track
    arcpy.Delete_management("Layer") # clean up to reuse layer name

arcpy.AddMessage("Finished bufffering, now merge")
arcpy.Merge_management(TempFCs,OutFC)

arcpy.AddMessage("Now cleaning up")
for ThisTemp in TempFCs:
    arcpy.Delete_management(ThisTemp) # delete the temp feature class

Creates a feature layer based on a definition query of a range of FIDs' in the shapefile then buffers to temp, after buffering merge and delete temp files. You should have no problem dissolving the geometries after the merge if that's what you require; as it stands each range will be a single multipolygon.

share|improve this answer
    
This is essentially what we have ended up working with after a handful of meetings to discuss the methods. Essentially, it was preferred because it was easier for people to understand than creating a grid, and therefore in future if some other developer has to reverse engineer the script they'll know what's going on... –  Erica Sep 8 at 12:29

I agree with @Chris W that it is unusual that the buffer operation is running out of memory. Regardless, there are a few routes to take. First, run the buffer tool as a 64-bit process. Second, try the attached script. The general workflow:

  1. Create fishnet of the same extent as the input points
  2. Iterate through the fishnet features and select the points within the feature
  3. Buffer the selected point features
  4. Merge the buffer shapefiles

A few notes here. It is faster to reduce the number of fishnet rows/columns and the subsequent number of output buffer shapefiles. As it is currently configured, the script splits the points into 4 equal pieces and then performs the buffer/merge.


import arcpy, os

arcpy.env.overwriteOutput = True

shp = r'C:\temp\myshp.shp'
outws = r'C:\temp\outws'

def createFN(x):
    """Create the fishnet in_memory"""
    arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(32615)
    desc = arcpy.Describe(x)
    xmax = desc.extent.XMax
    xmin = desc.extent.XMin
    ymin = desc.extent.YMin
    ymax = desc.extent.YMax

    coord = "%s %s" % (xmin, ymin)
    yaxis = "%s %s" % (xmin, ymax)
    corner = "%s %s" %(xmax, ymax)

    arcpy.CreateFishnet_management (r'in_memory\fishnet', coord, yaxis, 0, 0, 2, 2, corner, "NO_LABELS", "", "POLYGON")

createFN(shp)

def selectNBuffer(x, y):
    """Iterate through fishnet features and select/buffer points"""
    if not os.path.exists(y): os.makedirs(y)
    arcpy.MakeFeatureLayer_management(x, "shpLyr")
    mergeList = []
    with arcpy.da.SearchCursor(r'in_memory\fishnet', ["OID@", "SHAPE@"]) as cursor:
        for row in cursor:
            arcpy.SelectLayerByLocation_management("shpLyr", "INTERSECT", row[1])
            if arcpy.GetCount_management("shpLyr").getOutput(0) > 0:
                arcpy.Buffer_analysis("shpLyr", os.path.join(y, "buffer" + str(row[0]) + ".shp"), 1)
                mergeList.append(os.path.join(y, "buffer" + str(row[0]) + ".shp"))
                arcpy.Merge_management(mergeList, r'C:\temp\mergeOut.shp')
    arcpy.Delete_management(y)
share|improve this answer
2  
I like that @Aaron, if you use a shapefile grid you could even run multiple instances on a subselection - handy! One thing you would need to watch is that the "shpLyr" actually has records after the select layer by location operation - it is possible that there are empty grid squares; some tools don't like empty inputs though I can't remember if buffer is one of them. –  Michael Miles-Stimson Sep 3 at 21:40
    
@Michael Good call! I added some logic to the script that will fix the issue you brought up. –  Aaron Sep 3 at 21:49

the quick and dirty solution is to use a temporary feature class (using copy with your layer or select). You will then have new FID's and you can iteratively select bunches of 1000 points for the buffer tool.

A better solution (but a bit longer to write), is to use arcpy.da.searchcursor for reading your points, and arcpy.da.insertcursor for writing the buffer (you buffer the points one by one, see geometry object for details). You can set a condition to delete your insertcursor and create a new one every 1000 points.

share|improve this answer

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.