Zonal Statistics Crashes Randomly

[ColinLindeman]

Aloha,

I am having issues with the Spatial Analyst tool ZonalStatisticsAsTable in that it appears to crash the script for no reason.

The input zones (66 individual zones) is a raster file of the same extent, projection and cell size as the value raster. Some runs the script will iterate anywhere from 1 to 32 times (it needs to go the full 72 iterations) before crashing.

I also notice that the output table of the last zonal stats isn't complete when it crashes, it gets about 15 records (of 66) into the statistics table.

Thanks for any feedback and suggestions!

-Colin


############################################################

import string, datetime, arcpy

# set scratchworkspace... for desktop use, not ags use...
arcpy.env.scratchWorkspace = "C:/temp/scratch/"
# set overwrite outputs to true
arcpy.env.overwriteOutput = 1

# Variables...
netdcf_source_folder = "C:/mapservice_data/wrf/net_cdf_files/"
workspace_folder = "%SCRATCHWORKSPACE%/"
workspace_folder_filegdb = workspace_folder + "scratch.gdb/" 
sde_workspace_d2p2 = "C:/mapservice_files/db_connections/connection.sde/"
InMemory_netcdf_raster = "InMemory_netCDF_File"
table_template = "C:/mapservice_data/WRF/wrf.gdb/wrf_table_template_rain"
in_zone_data = "C:/mapservice_data/WRF/wrf.gdb/wrf_zones"
model_runtime = '"' + str(datetime.datetime.now()) + '"'

# Variables based on the user input variable...
input_netcdf_file = netdcf_source_folder + "RAIN_vn10km.nc"
target_table = workspace_folder_filegdb + "wrf_rain_v2"
out_netcdf_table_view = "wrf_rain_table_view"
sde_table = sde_workspace_d2p2 + "D2P2.wrf_rain_v2"
time_table = []

# Check out the spatial analyst extension...
arcpy.CheckOutExtension("Spatial")

try:
 # Create a 'catch all' table for the statistics, this will be copied to SDE
 # Copy Rows...CopyRows_management (in_rows, out_table, {config_keyword})
 arcpy.CopyRows_management(table_template, target_table)
except:
 print "copy rows fail", arcpy.GetMessages()

try:
 # Create table of time slices from netcdf file...
 # MakeNetCDFTableView_md (in_netcdf_file, variable, out_table_view, {row_dimension}, {dimension_values}, {value_selection_method})
 arcpy.MakeNetCDFTableView_md(input_netcdf_file, "rain", out_netcdf_table_view, "time", "", "")
except:
 print "time slice table fail", arcpy.GetMessages()
try:
 # SearchCursor to get time values...
 # SearchCursor (dataset, {where_clause}, {spatial_reference}, {fields}, {sort_fields})
 rows = arcpy.SearchCursor(out_netcdf_table_view, "", "", "time", "")
 currentTime = ""
 for row in rows:
  if currentTime != row.time:
   currentTime = row.time
  time_table.append(row.time)
except:
 print "time_table fail", arcpy.GetMessages()

try:    
 # Create a list of index numbers in the netcdf file to process...
 x = 0
 for time in time_table:
  if x < 10:
   y = "0" + str(x)#prefixes a "0" to numbers less than 10 to help in sorting.
  else:
   y = str(x)
  output_netcdf_raster = workspace_folder + "rain_" + y + ".img"
  zonalst_table = workspace_folder_filegdb + "rain_" + y + "_zonalst"
  time_slice_number = time
  
  try:
   # Make NetCDF Raster Layer...MakeNetCDFRasterLayer_md (in_netcdf_file, variable, x_dimension, y_dimension, out_raster_layer, {band_dimension}, {dimension_values}, {value_selection_method})
   arcpy.MakeNetCDFRasterLayer_md(input_netcdf_file, "rain", "longitude", "latitude", InMemory_netcdf_raster, "", time_slice_number, "BY_VALUE")
  except:
   print "Make NetCDF Raster Layer", time_slice_number, arcpy.GetMessages()

  try:
   # Copy Raster for zonal stats input...CopyRaster_management (in_raster, out_rasterdataset, {config_keyword}, {background_value}, {nodata_value}, {onebit_to_eightbit}, {colormap_to_RGB}, {pixel_type})
   arcpy.CopyRaster_management(InMemory_netcdf_raster, output_netcdf_raster,"","","","","","")
  except:
   print "Copy Raster", time_slice_number, arcpy.GetMessages()

  try:
   # Zonal Stats as table...ZonalStatisticsAsTable (in_zone_data, zone_field, in_value_raster, out_table, {ignore_nodata}, {statistics_type})
   arcpy.sa.ZonalStatisticsAsTable(in_zone_data, "D_NAME", output_netcdf_raster, zonalst_table, "DATA", "SUM")
  except:
   print "zonal stats as table", time_slice_number, arcpy.GetMessages()

  try:
   # Process: Add Field...AddField_management (in_table, field_name, field_type, {field_precision}, {field_scale}, {field_length}, {field_alias}, {field_is_nullable}, {field_is_required}, {field_domain})
   arcpy.AddField_management(zonalst_table, "INDEXSLICE", "SHORT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
  except:
   print "add field indexslice", index_slice_number, arcpy.GetMessages()
  try:
   # Process: Calculate Field...CalculateField_management (in_table, field, expression, {expression_type}, {code_block})
   expression = y
   arcpy.CalculateField_management(zonalst_table, "INDEXSLICE", expression, "PYTHON")
  except:
   print "calc field indexslice", index_slice_number, arcpy.GetMessages()

  try:
   arcpy.AddField_management(zonalst_table, "RUNTIME", "TEXT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
  except:
   print "add field runtime", time_slice_number, arcpy.GetMessages()
  try:
   arcpy.CalculateField_management(zonalst_table, "RUNTIME", model_runtime, "PYTHON")
  except:
   print "calc field runtime", time_slice_number, arcpy.GetMessages()

##  try:
##   arcpy.AddField_management(zonalst_table, "TIMESLICE", "DATE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
##  except:
##   print "add field timeslice", time_slice_number, arcpy.GetMessages()
##  try:
##   arcpy.CalculateField_management(zonalst_table, "TIMESLICE", time_slice_number, "PYTHON")
##  except:
##   print "calc field timeslice", time_slice_number, arcpy.GetMessages()
##
##  try:
##   arcpy.AddField_management(zonalst_table, "TIMESLICETXT", "TEXT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
##  except:
##   print "add field timeslicetxt", time_slice_number, arcpy.GetMessages()
##  try:
##   arcpy.CalculateField_management(zonalst_table, "TIMESLICETXT", time_slice_number, "PYTHON")
##  except:
##   print "calc field timeslicetxt", time_slice_number, arcpy.GetMessages()    

  try:
   # Process: Append...Append_management (inputs, target, {schema_type}, {field_mapping}, {subtype})
   arcpy.Append_management(zonalst_table, target_table, "NO_TEST", "", "")
  except:
   print "append", time_slice_number, arcpy.GetMessages()

  print y
  # add 1 to the counter
  x += 1

## #should add schema lock test here...
## TestSchemaLock (dataset)
## arcpy.TestSchemaLock()
 
## try:
##  # Copy the statistics table to SDE, overwriting the old table...
##  # Copy Rows...CopyRows_management (in_rows, out_table, {config_keyword})
##  arcpy.CopyRows_management(target_table, sde_table)
## except:
##  print "copyrows", time_slice_number, arcpy.GetMessages()

except:
 print "general failure", arcpy.GetMessages()

# Check in the spatial analyst extension...
arcpy.CheckInExtension("Spatial")

[ChrisSnyder]

Make sure the zone raster and/or the value raster has a raster attribute table. If not, use the  BuildRasterAttrbibuteTable (or is it called CreateRasterAttributeTable? I forget).

Generally, if the raster has more than 500 (I think) values, the raster attribute table is not automatically created.

[ColinLindeman]

I wanted to let others know that i have found a workaround to my problem, although I would rather it work how it was intended to.

There is a bug report for this and its nimbus id is NIM063814, and here is a link
http://resources.arcgis.com/content/nimbus-bug?bugID=TklNMDYzODE0

I had used a os.spawnv to perform the "Zonal Statistics as Table" for each individual raster, which kept the memory isolated and ran without any crashing.

I will share the scripts here soon or send me an email if they aren't up yet.

thanks
-Colin


edit:
here are some links that got me on the track to use os.spawnv
http://forums.esri.com/Thread.asp?c=93&f=1729&t=177007
http://forums.esri.com/Thread.asp?c=93&f=1729&t=173610&mc=33

[ChrisSnyder]

In the olden days you used to be able to run the SA tools (like zonal stats) in a loop for weeks with no memory leaks. Whatever ESRI did to some of the SA tools in v9.3 and v10 seemed to:

1. Increase the speed of execution (sometimes by a whole lot). We all appreciate that...
2. Cause many of these algorithms to run out of memory. Yar...

One example is the Fill tool. In years past, the thing might take a week to finish, but it finished. Now the stupid things just bombs with an 'out of memory' error. Hmm...

[NilsBabel]

I wanted to let others know that i have found a workaround to my problem, although I would rather it work how it was intended to.

There is a bug report for this and its nimbus id is NIM063814, and here is a link
http://resources.arcgis.com/content/nimbus-bug?bugID=TklNMDYzODE0

I had used a os.spawnv to perform the "Zonal Statistics as Table" for each individual raster, which kept the memory isolated and ran without any crashing.

I will share the scripts here soon or send me an email if they aren't up yet.

thanks
-Colin


edit:
here are some links that got me on the track to use os.spawnv
http://forums.esri.com/Thread.asp?c=93&f=1729&t=177007
http://forums.esri.com/Thread.asp?c=93&f=1729&t=173610&mc=33

Hey Colin, I think I'm coming across the same exact issue.  I'm running ZonalStatisticsAsTable on a bunch of rasters using ListRasters.  It gets through about 4 or 5 before it crashes.  Can you post some of your code that shows how you used the spawnv?  BTW I came across those forum posts about spawnv also but noticed they're from about 2005.  I can't believe we are still dealing with this issue.

Thanks in advance,

Nils

[DagnachewLegesse]

I am having this same issue with ArcGIS 10 SP2.
Any solution yet?
What is strange is the script works just fine for about 10-15 files while looping through 1200 RASTER GRIDs of climate data and it crashes with the following. If I rerun it again, it works again but crashes at another file (may be after 40 files). It is totaly random and I don't know why.

Traceback (most recent call last):
  File "C:\DATA2011\scripts\zonalStat_CRU.py", line 36, in <module>
    outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, infile, outTable, "NODATA", "MEAN")
  File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 5837, in ZonalStatisticsAsTable
    statistics_type)
  File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 5829, in wrapper
    statistics_type)
  File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\_base.py", line 474, in <lambda>
    return lambda *args: val(*gp_fixargs(args))
ExecuteError: ERROR 999999: Error executing function.
Workspace or data source is read only.
Workspace or data source is read only. [The C:\DATA2011\DATA\CRU3\hru\tmn\ workspace is read only.]
ERROR 010067: Error in executing grid expression. Zonal statistics program failed
Failed to execute (ZonalStatisticsAsTable).

Here is my code:


[PHP]import arcpy
import os
from arcpy import env
from arcpy.sa import *
try:
    #def zonal_stat(indir):
    #Set scratch workspace
    arcpy.env.scratchWorkspace = "c:\\mytempfile.gdb"
    arcpy.env.overwriteOutput = True
    arcpy.CheckOutExtension("Spatial") #chaeckout spatial analyst extension
    # Set environment settings
    basedir = "C:\\DATA2011\\DATA\CRU3\\resampled\\tmn\\"
    indirs = ['tmn\\', 'vap\\', 'cld\\']
    #Set the feature class zone data
    inZoneData = r'C:\DATA2011\Model_500\bnile_hruf.shp'
    zoneField = "GRIDCODE"
    for indir in indirs:
        env.workspace = basedir + indir
        outtemp = env.workspace       
        outpath = outtemp.replace('resampled', 'hru')
        rasterList = arcpy.ListRasters("*", "GRID")
        if not os.path.exists(outpath):
            os.makedirs(outpath)
       for infile in rasterList:
            print "Zonal Stat  ..." + outpath + infile
            outTable = outpath + "\\CRU_" + infile + ".dbf"
            outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, infile, outTable, "NODATA", "MEAN")
except:
     print arcpy.GetMessages(2)[/PHP]

[ChrisBater]

I've found that first manually (or programmatically) converting my zone polygons to rasters improves chances of success. Definitely not a reliable solution though....

[ChrisSnyder]

When run in a loop, the Watershed SA tool will also randomly fail if the input "pour points" are in vector format. If you 1st convert the points to a GRID, the process become stable.

[ChrisBater]

I wanted to let others know that i have found a workaround to my problem, although I would rather it work how it was intended to.

There is a bug report for this and its nimbus id is NIM063814, and here is a link
http://resources.arcgis.com/content/nimbus-bug?bugID=TklNMDYzODE0

I had used a os.spawnv to perform the "Zonal Statistics as Table" for each individual raster, which kept the memory isolated and ran without any crashing.

I will share the scripts here soon or send me an email if they aren't up yet.

thanks
-Colin


edit:
here are some links that got me on the track to use os.spawnv
http://forums.esri.com/Thread.asp?c=93&f=1729&t=177007
http://forums.esri.com/Thread.asp?c=93&f=1729&t=173610&mc=33

I wonder, would writing and then calling external modules accomplish the same thing as os.spawnv in terms of memory management?

See for example: http://docs.python.org/tutorial/modules.html

[DagnachewLegesse]

I wonder, would writing and then calling external modules accomplish the same thing as os.spawnv in terms of memory management?

See for example: http://docs.python.org/tutorial/modules.html

Thanks Chris for the suggestion about using RASTER for the zones. I will give it a try. As to the external module suggestion, as you can see in my code above, I have tried that too (commented line). Will also give that a more try with the raster zone and report back.

have a good day
Dag