Tuesday 24 July 2018

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


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'


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...