4

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'
  • I don't understand why you don't process the ZonalStatisticsAsTable as it is, setting polygon shapefile as zone data and the raster as input, getting as output a simple table which countains a row for every polygon with statistics. I had a look in your script, it's not well organized, moreover sshednutret90 variable is not defined. – Nonpenso Jul 11 '14 at 8:09
  • I would, but when I do that zonal stats are not calculated for every polygon (only 486 out of 3897 total). The hope was that this code would cycle through every row and make sure that zonal stats are calculated for each zone. thanks for the indication re: sshednutret90 by the way, I'll fix it. – Fred Jul 11 '14 at 14:36
  • I never used Zonal with overlapped polygons. As workaround I can suggest you to create a single shapefile for every polygon thru a loop, then call the ZonalStatistics. If you have patience I can write down a script on Monday, trying with overlapped polygons. I'm curious to undestand if it is a bug. – Nonpenso Jul 11 '14 at 14:51
  • hm the single shapefile for each polygon may be a path to follow, I'm just hoping there's a faster/less memory intensive way of doing it using a cursor. It's kind of you to offer, I'm going to try to solve this before the weekend but will be sure to post back with how it goes! – Fred Jul 11 '14 at 15:04
  • @Nonpenso, so the code now runs, but returns empty tables. I'm giving it a rest for the day and will work on it with a clear mind later. – Fred Jul 11 '14 at 20:17
3

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
  • this is really great, @Nonpenso! thanks for the help. I'm getting an error telling me that rec is undefined, may I ask what the difference between rec and recs is? could one not just write recs for both? it's running well until the 990th FID and then returning the error. – Fred Jul 14 '14 at 19:40
  • rec and recs (like row and rows) are the cursor variables to iterate in the table. recs can be defined as a list of rows, while rec is used to iterate for every row. If rec is not defined, it means UpdateCursor failed. Please check if temp_tab and/or 'FID_' field exist. – Nonpenso Jul 15 '14 at 7:37
  • Have a look on the 991st polygon. Probably it cannot create a ZonalStatistics table, maybe because it is too small or out of the raster extent. Anyway without or with empty temp_tab the loop with rec returns an error. – Nonpenso Jul 15 '14 at 9:40
  • temp_tab and FID_ exist indeed. the 991st polygon is the smallest in the set, but there are others that small that have been run through the code successfully before it reaches 991. You're right, it is outside of the extent of the raster. Currently working on adding some code to put FIDs outside of the extent in a different table, I'll post it when i've got something that seems to work. thanks again for all the help @nonpenso – Fred Jul 15 '14 at 14:48
  • @Fred, very good. Don't forget to vote for the answer (it has 0 votes => unuseful) and eventually close the question when the script definitely works. Thanks – Nonpenso Jul 15 '14 at 15:08
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/

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.