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:
- Change the processing extent to 400m around the point
- Find out if the point is within 400m of the edge of an urban area
- If more than one tile involved, stitch tiles together
- Extract the raster data in the 100, 200,300,400m (the data is binary thus only 1s, 0s)
- Grab the count of 1s in the extracted raster
- 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 + " = 1 " if arcpy.Exists(outFC): arcpy.AddMessage(" Begin Deleting " + outFC + " ...") arcpy.Delete_management(outFC) arcpy.CopyFeatures_management(inFC, outFC) arcpy.AddField_management(outFC, "D_100", "DOUBLE") arcpy.AddField_management(outFC, "D_200", "DOUBLE") arcpy.AddField_management(outFC, "D_300", "DOUBLE") arcpy.AddField_management(outFC, "D_400", "DOUBLE") with arcpy.da.UpdateCursor(outFC,pFields) as cursor: m=0 for primaryRow in cursor: arcpy.AddMessage(" Processing ID= " + str(primaryRow) + " ...") x = primaryRow y = primaryRow shape = primaryRow #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.AddMessage(" starting bufferanalysis ...") 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 for r in arcpy.da.SearchCursor("UrbanTiles", ("SHAPE@AREA")) if not r 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)+".tif" myInRaster=os.path.join(TSpace,myTileFullName) myRasterList.append(myInRaster) arcpy.AddMessage(" Your tiles are " + str(myRasterList)) #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: arcpy.AddMessage("processing raster:"+ str(raster)) 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)) primaryRow[4+i]=myCount cursor.updateRow(primaryRow) i += 1 else: primaryRow=-9999 primaryRow=-9999 primaryRow=-9999 primaryRow=-9999 cursor.updateRow(primaryRow)