ecodrive

Zonal Statistics Crashes Randomly

Discussion created by ecodrive on Mar 8, 2011
Latest reply on Oct 13, 2012 by forreststevens
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")

Outcomes