Making zonal analysis faster?

I created this python script which performs zonal analysis for the 100m, 200m, 300m, and 400m surrounding a list of points. It works but I have two main problems:

A) The raster I am processing has 1m cell size for only urban areas within the state of California. The original data is tiled because it is huge.

B) Right now it runs about 90 seconds per point and my problem is that I have 34,000 points. That means we are talking a month of processing time. I am hoping to somehow speed the process up.

My basic processing steps are:

1. Change the processing extent to 400m around the point
2. Find out if the point is within 400m of the edge of an urban area
3. If more than one tile involved, stitch tiles together
4. Extract the raster data in the 100, 200,300,400m (the data is binary thus only 1s, 0s)
5. Grab the count of 1s in the extracted raster
6. Rinse-repeat for each point.

I fear that that this is really not workable for 34,000 points.

I might be able to speed things up a little by pre-processing to remove Step #2, however we expect only 3-4,000 points to be kicked out. That helps but likely not enough.

I might be able to reduce Step #3 by pre-processing and kicking out all points that need more than one tile and processing them separately.

I am in a bit of a quandary and am hoping that perhaps my code is not very efficient or you GISers could tighten my code (below).

I know lots of people use QGIS - would that be any faster?

``````import arcpy, traceback, os, sys, time
import numpy as np
import itertools as it
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension('Spatial')

arcpy.env.overwriteOutput=True
arcpy.env.workspace = r"in_memory"

WSpace = r'C:\Users\Don\Documents\AirQuality\AirBufferAnalysis.gdb'  #  Where the file geodatabase is with your points in it
TSpace =  r'L:\AirQuality\TreeTiles\Tile_2012'                       #  Where the FRAP tiles reside

inFC           = os.path.join(WSpace, "SingleTBTestPoint")               #  Point file with the geocoded addresses
outFC          = os.path.join(WSpace, "TB_Processed")                #  Final Output Name
UrbanTileFC    = os.path.join(WSpace,"Urban_TreeIndex")              #  This is the combination of Urban Census and Tree Tile index
tempBuffer     = os.path.join(WSpace,"tempBuff")                     #  This is the temporary buffered point
tempIntersect  = os.path.join(WSpace,"tempIntersect")                #  This is the temporary intersection with Urban Tree and buffered point
rtemp = os.path.join(WSpace,"myRastTemp")
arcpy.env.snapRaster = os.path.join(TSpace,"R01_C01_TreeCanopy_CA_2012.tif")

pFields=["SHAPE@X", "SHAPE@Y", "SHAPE@", "id", "D_100","D_200","D_300","D_400"]
distances=[100,200,300,400]                                          #distances are in meters

fullDist=400
bufferUnit = "meters"

rFields=['Value','Count']
SQLWhereClause= " " + rFields[0] + " = 1 "

if arcpy.Exists(outFC):
arcpy.AddMessage("       Begin Deleting " + outFC + " ...")
arcpy.Delete_management(outFC)
arcpy.CopyFeatures_management(inFC, outFC)

with arcpy.da.UpdateCursor(outFC,pFields) as cursor:
m=0

for primaryRow in cursor:
arcpy.AddMessage("  Processing ID= " + str(primaryRow[3]) + " ...")
x = primaryRow[0]
y = primaryRow[1]
shape = primaryRow[2]
#set processing extent to a little bigger than the distance surrounding the point
arcpy.env.extent = arcpy.Extent(x-405, y-405, x+405, y+405)
arcpy.AddMessage("       extent " + str(arcpy.env.extent) + " ...")
#First buffer by 400 meters and set up the raster(s) needed
arcpy.Buffer_analysis(shape, tempBuffer, "%s Meters" % str(fullDist))
arcpy.DefineProjection_management(tempBuffer, arcpy.SpatialReference(3310))
inIntersectFeatures = [tempBuffer, UrbanTileFC]
arcpy.Intersect_analysis(inIntersectFeatures, "myPointIntersect")
arcpy.MakeFeatureLayer_management("myPointIntersect", "UrbanTiles")
#check if the entire 400m surrounding the point is within a feature inside of UrbanTiles
area = sum([r[0] for r in arcpy.da.SearchCursor("UrbanTiles", ("SHAPE@AREA")) if not r[0] is None])
if area > 502650:
myRasterList = []
#find out which raster tiles are within 400m of point
rows = arcpy.da.SearchCursor("UrbanTiles", "TILENAME")
for r in rows:
myTileFullName=str(r[0])+".tif"
myInRaster=os.path.join(TSpace,myTileFullName)
myRasterList.append(myInRaster)
#tempRasterList=[os.path.join(TSpace,"R09_C07_TreeCanopy_CA_2012.tif"),os.path.join(TSpace,"R09_C06_TreeCanopy_CA_2012.tif")]
#arcpy.MosaicToNewRaster_management(tempRasterList, WSpace, "inRaster", arcpy.SpatialReference(3310), "1_BIT", "1", "1", "LAST","FIRST")
i=0
#stitch together the raster tiles to then extract data from
for raster in myRasterList:
out1 = Con(IsNull(raster),0,raster)
if i==0:
out2 = out1
i += 1
else:
out2 = out2 + Con(IsNull(raster),0,raster)
i += 1
#out2.save(os.path.join(WSpace,'myOut2'))
i=0
for d in distances:
rast_d = ExtractByCircle(out2, arcpy.Point(x,y),d,"INSIDE")
if arcpy.Exists(rtemp):
arcpy.Delete_management(rtemp)
rast_d.save(rtemp)
for myCount in arcpy.da.SearchCursor(rtemp, rFields, where_clause=SQLWhereClause):
arcpy.AddMessage("       For distance " + str(d) + " count is " + str(myCount[1]))
primaryRow[4+i]=myCount[1]
cursor.updateRow(primaryRow)
i += 1
else:
primaryRow[4]=-9999
primaryRow[5]=-9999
primaryRow[6]=-9999
primaryRow[7]=-9999
cursor.updateRow(primaryRow)
``````
• have you tried using a profiling tool to identify which parts of your script are slow? That would allow you to focus on improving the relevant section(s) Commented Jan 27, 2016 at 2:16
• I did not know there was one, I will try and report back! Commented Jan 27, 2016 at 16:58