1

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)
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[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.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[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)
            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[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)
  • 1
    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) – Stephen Lead Jan 27 '16 at 2:16
  • I did not know there was one, I will try and report back! – user918967 Jan 27 '16 at 16:58
1

Rather than running the zonal statistics multiple times on every point (which involves stitching the rasters together, then throwing them away each time) it may make more sense to run the operation using the entire point layer, using something like this:

  • buffer the point layer by 400m
  • determine which rasters fall within the buffer
  • stitch those rasters together into a single raster
  • perform the zonal statistics on the single raster
  • I feel that one huge zonal stats of 34,000 points against a 1m of CA will completely crush my computer. I may be able to do the buffer of the 34,000 points first, then zoom into each one and punch out the zonal stats... – user918967 Jan 27 '16 at 17:00
  • My gut feeling (without having tested this) is that it's the opposite - 34k points in one featureclass isn't exceptional, and ZonalStatistics should be able to handle that more easily than 34k individual runs. Give it a test over the weekend? – Stephen Lead Jan 27 '16 at 21:57

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.