1

I have a script that does some processing but the main part are two things

  1. run a multiple ring buffer using arcpy.MultipleRingBuffer_analysis and four rings
  2. do some zonal statistics using arcpy.sa.ZonalStatisticsAsTable on the four rings

I do this repetitively (through a loop) so I want to dump the results to a table called Final Results. It makes sense to me to have the results organized in five columns (PatientID, Distance1, Distance2, Distance3, Distance4). I making my loop off a search cursor and keeping the PatientID.

I thought to create a pivot table using arcpy.PivotTable_management where the result becomes one row and five columns (the pivot variable and four buffers) and then push that to Final Results. However, it appears I would still need to create a search cursor for my pivot table, grab the first (and only) row, then create an insert cursor to put the data into my Final Resultstable.

Alternatively, I suppose I could just go to the table from #2 above, create a search cursor, and assign four variables (D1, D2, D3, D4) to the zonal stats results with some sort of switch-case statement, then create an insert cursor that concatenates D1, D2, D3, D4 and uses insert row.

Both seem rather inefficient.

Is there a smarter way to just push my results from zonal analysis into a table that I use?

1

You cannot use multiple rings because it produces overlapping circles, no good for zonal statistics. Try this one:

import arcpy, traceback, os, sys
import itertools as it
from arcpy import env
env.overwriteoutput=True
env.workspace=r'in_memory'

points = r'D:\Scratch\points.shp'
pFields=['Shape@',"DISTANCE1","DISTANCE2","DISTANCE3","DISTANCE4"]
distances=[5,10,15,20]
raster=r'D:\Scratch\dem_ext'
victRaster=r'in_memory\rtem'

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)
    pgon=arcpy.Geometry()
    with arcpy.da.UpdateCursor(points,pFields) as cursor:
        m=0
        for row in cursor:
            point=row[0]
            for i in range(4):
                buf=point.buffer(distances[i])
                anExtent=buf.extent
                envelope='%f %f %f %f' %(anExtent.XMin, anExtent.YMin, anExtent.XMax, anExtent.YMax,)
                arcpy.Clip_management(raster,envelope,victRaster,buf,"-3.402823e+038","ClippingGeometry")
                extract = arcpy.RasterToNumPyArray(victRaster,"","","",-9999)
                chain = it.chain.from_iterable(extract)
                reduced=filter(lambda x: x not in [-9999],chain)
                mean=sum(reduced)/float(len(reduced))
                row[1+i]=mean
            m+=1
            cursor.updateRow(row)
            arcpy.AddMessage(m)
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()       

it took 2.5 minutes to process 1000 points. Results stored in parent table, you have to modify inputs, e.g. points, fields, raster. Script assumes that fields exist.

  • Interesting script but I am not so sure it will work for me. I tested it and it took . Distances are 100m, 200m, 300m, and 400m and the cell size of the raster is 1m. – user918967 May 12 '15 at 6:00
  • Hi Felix, interesting script and approach but I am not so sure it will work for me. I tested it and it took 17 minutes to process one point and it bombed with an ERROR 999998: Unexpected Error in the RasterToNumPyArray line. Distances are 100m, 200m, 300m, and 400m and the cell size of the raster is 1m. BTW, I did have to change the "in_memory" to my file system (e.g. the C drive). – user918967 May 12 '15 at 6:12
  • Replace numpyarray by raster statistics – FelixIP May 12 '15 at 9:43
  • I think it was another issue. The points were in a different projection than the raster layer. I projected the points and it ran smoking fast! – user918967 May 12 '15 at 14:57
  • one last question. the raster layer is just 0's and 1's. If I delete the "/float(len(reduced))" I thought it would sum up the raster layer and tell me how many square meters are in the buffer. Instead it just gives me the number of square meters of the buffer (e.g. the size of the buffer). How do I alter it to give me the the area of 1s? – user918967 May 12 '15 at 14:59

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.