Calculating zonal statistics of raster data in multiple overlapping zones and combining them into one table

Question

The aim here is to calculate the max, min, mean, range, sum, count and std of raster data contained in areas delimited by polygons (much like what the Zonal Statistics as a Table tool from Spatial Analyst does).

The difference being that many of my zones overlap (they are hydrological servicesheds, based on watersheds among other factors), and that I would like to have as output a table with the basic stats for the data within each polygon (one row per polygon). Double counting is not an issue.

Below is the code I have for the moment, which defines a function through which zonal stats are run and the results put into a table.

This is run using a cursor which goes through every row in the dataset which contains every serviceshed (one on each row, 3897 in total). The code does not work, and I was hoping to get guidance from the GIS SE community! any input, be it on method or the code, is greatly welcomed!

I'm using ArcMap 10.2 and have access to an advanced license.

import arcpy
arcpy.env.workspace='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014'
import arcpy.sa
import arcpy.da
import os
arcpy.env.overwriteOutput = True
arcpy.CreateTable_management('C:\\Users\\xxx\\Desktop\\ESRI_summer2014\\zonaldata.gdb', 'sshed_nxport90')
sattable = 'C:\\Users\\xxx\\Desktop\\ESRI_summer2014\\zonaldata.gdb\\sattable'
def zonalstatsfun(row,start,end):
    zone_field="FID"
    value = 'C:\\Users\\xxx\\Desktop\\ESRI_summer2014\\invest_workspace\\nut_ret_1990\\output\\n_export_1990.tif' #raster dataset holding values to be extracted
    arcpy.sa.ZonalStatisticsAsTable(row, zone_field, value, sattable,"","ALL")
    arcpy.Append_management(stattable,sshednutret90)
    arcpy.DeleteRows_management(sattable) #cleanup
    print 'done serviceshed FID {0}'.format(int(in_field)) #to keep track while it is running

fc='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\servicesheds_2014-06-24\\servicesheds_v0.shp'
cursor=arcpy.da.SearchCursor(fc, "FID")
for row in cursor:
    zonalstatsfun(row,0,3896)

MORE RECENT VERSION OF CODE BELOW*

#the aim of this code is to extract data from a raster for each polygon in a feature class where polygons overlap
#this will permit the user to aggregate raster data to overlapping subwatershed or serviceshed level

import arcpy
arcpy.env.workspace='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014'
import arcpy.sa
import arcpy.da
import os
arcpy.env.overwriteOutput = True

sattable = 'C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\zonaldata.gdb\\sattable'
#sshed_nxport90=arcpy.CreateTable_management('C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\zonaldata.gdb','sshed_nxport90', sattable)
def zonalstatsfun(row,start,end):
    value = 'C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\invest_workspace\\nut_ret_1990\\output\\n_export_1990.tif' #XXXX = raster dataset holding values to be extracted
    table='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\zonaldata.gdb\\sshed_nxport90'
    FID=("FID")
    arcpy.sa.ZonalStatisticsAsTable(fc, FID, value, sattable,"DATA","ALL")
    arcpy.Append_management(sattable,table,'TEST')
    ##arcpy.DeleteRows_management(sattable) 
    print 'appended'

fc='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\servicesheds_2014-06-24\\servicesheds_v0.shp'
cursor=arcpy.da.SearchCursor(fc, "FID")
for row in cursor:
    zonalstatsfun(row,0,5) #til total 2896

****FURTHER UPDATE****** The follow code gives good results, but repeats itself indefinitely and it seems the (row, start, end) statement is useless. even with values (row, 0,0) it will repeat until I force quit ArcMap.

Can anyone find a mistake in the code I have missed?

import arcpy
arcpy.env.workspace='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014'
import arcpy.sa
import arcpy.da
import os
arcpy.env.overwriteOutput = True

def zonalstatsfun(FID,start,end):
    value='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\invest_workspace\\sedimentretention1990\\output\\rkls_1990.tif' #XXXX = raster dataset holding values to be extracted
    table='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\zonaldata.gdb\\sshed_rkls90'
    FID="FID"
    arcpy.sa.ZonalStatisticsAsTable(fc, FID, value, statable,"DATA","ALL")
    arcpy.Append_management(statable,table,'NO_TEST')
    print 'appended'
statable=arcpy.CreateTable_management('C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\zonaldata.gdb','statable', 'C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\socioecon_settlements.gdb\\USLE_sshed90')
sshed_rkls90=arcpy.CreateTable_management('C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\zonaldata.gdb','sshed_rkls90', 'C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\zonaldata.gdb\\statable')
fc='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\servicesheds_2014-06-24\\servicesheds_v0.shp'
cursor=arcpy.da.SearchCursor(fc, "FID")
for FID in cursor:
    zonalstatsfun(FID,0,0) #til total 2896

-----------ANOTHER EDIT------------

Here is a new version of the code, which runs MUCH faster and gets the job done. Only one issue with it: It runs infinitely (until ArcMap crashes). I've tried implementing a range in many different ways but none seem to work.

Anybody have an idea of how I could limit the FIDs to either run all once or manually set the range to FID range(0,3896)?

import arcpy
from arcpy import da, sa, os
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True

workdir = r'C:\Users\frederichoffmann\Desktop\ESRI_summer2014'
zonaldb = workdir + r'\zonaldata.gdb'
tab_template = workdir + r'\zonaldata.gdb\template_sum'

statable=arcpy.CreateTable_management(zonaldb,'statable', tab_template)
values=arcpy.CreateTable_management(zonaldb,'sshed_rkls90', statable)
fc='C:\\Users\\frederichoffmann\\Desktop\\ESRI_summer2014\\servicesheds_2014-06-24\\servicesheds_v0.shp'
cursor=arcpy.da.SearchCursor(fc, "FID")

for FID in cursor:
    0<=FID<=5
    value=workdir + r'\invest_workspace\sedimentretention1990\output\rkls_1990.tif' #XXXX = raster dataset holding values to be extracted
    table= zonaldb + r'\sshed_rkls90'
    FID="FID"
    arcpy.sa.ZonalStatisticsAsTable(fc, FID, value, statable,"DATA","SUM")
    arcpy.Append_management(statable,table,'NO_TEST')
    print 'appended'

Answer 1

I tested the ZonalStatistics function with overlapped polygons shapefile as zone data. Shapefile is converted into a raster, so overlapped areas are lost.

As workaround I tried to extract every polygon feature from the shapefile and process the ZonalStatistics. The script does not contain def.

Summing up what it does:

• Define the input variables
• Create the output table
• Iterate for every polygon feature of the zonal shapefile
• Extract a polygon creating a temp shapefile
• Compute ZonalStatistics in a the temp DBF table
• Add the FID of original polygon to temp table
• Append the row to the output table
• Delete temp shape and table

This is the code.

import arcpy
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True

workdir = r'C:\Users\frederichoffmann\Desktop\ESRI_summer2014'
inraster = workdir + r'\invest_workspace\sedimentretention1990\output\rkls_1990.tif'
zonal_shp = workdir + r'\servicesheds_2014-06-24\servicesheds_v0.shp'
tab_template = workdir + r'\socioecon_settlements.gdb\USLE_sshed90'

stat_table = workdir + r'\zonaldata.gdb\statable'
arcpy.CreateTable_management(workdir + r'\zonaldata.gdb','statable', tab_template)

with arcpy.da.SearchCursor(zonal_shp, ['FID']) as rows:
    for row in rows:

        fid = row[0]

        expression = '"FID" = ' + str(fid)
        temp_shp = workdir + r'\tempshp.shp'
        arcpy.Select_analysis(zonal_shp, temp_shp, expression)

        temp_tab =  workdir + r'\temptab.dbf'
        arcpy.sa.ZonalStatisticsAsTable(temp_shp, 'FID', inraster, temp_tab, "DATA", "ALL")

        with arcpy.da.UpdateCursor(temp_tab, ['FID_']) as recs:
            for rec in recs:
                rec[0] = fid
                recs.updateRow(rec)
        del rec, recs

        arcpy.Append_management(temp_tab, stat_table, 'NO_TEST')

        arcpy.Delete_management(temp_tab)
        arcpy.Delete_management(temp_shp)

del row, rows

Answer 2

I found an ESRI supplemental tool called Zonal Statistics As Table 2, which allows me to calculate overlapping polygons. http://blogs.esri.com/esri/arcgis/2013/11/26/new-spatial-analyst-supplemental-tools-v1-3/