Kriging and Zonal Statistics for multiple rasters using Python

[DavidBetts]

Branched from: Calculate zonal mean for multiple rasters using Python 

I found a solution to the same problem.  This solution will allow you to process all of the raster files in a selected folder instead of naming the raster files individually.  I had 936 rasters I needed to analyze, so the ability to loop through all of the available rasters was a necessity for me. 

The script below is my modificaiton of the discussion in the following post: arcpy - Python Script Zonal Stats as Table Loop Question - Geographic Information Systems Stack Exch... 

The output from this script is an individual .dbf file (and associated files) for each of the rasters that you process.

import arcpy, os, arcinfo
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

# Select the folder containing raster files.  This script will use ALL of
# the raster files in the selected folder. 
env.workspace =
   "C:\\Users\\davebetts\\GIS\\project\\rasters"

# Select the shapefile containing the polygons to use as boundaries
# for zonal statistics
watershedFeat =
   "C:\\Users\\davebetts\\GIS\\project\\zones_polygons.shp"

# Select output folder for saving the output - zonal tables (.dbf files)
outDir =
   "C:\\Users\\davebetts\\GIS\\project\\rasters\\tables\\"

# Something goes wrong with this script during use, perhaps with the
# temporary files.  No error messages are given.  The "print" statements
# inserted within the script keep track of where to restart.  Replace the '0'
# in the "for" statement with the most recently printed integer (printing of
# the variable 'ndx').

x = arcpy.ListRasters()
# the "0" can be replaced by the most recent result of
# "print ndx" in order to restart where the code stopped
for raster in arcpy.ListRasters()[0:]:
   print raster 
    ndx = x.index(raster)
    print ndx
    outTable = outDir + raster + ".dbf"
    arcpy.gp.ZonalStatisticsAsTable_sa(watershedFeat,
      "FID", # Select an attribute in the shape file to identify polygons
      raster,
      outTable,
      "NODATA","MEAN")

arcpy.CheckInExtension("Spatial")
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

!!! This script is still stopping several times partway through the total number of rasters I'm processing.  I have to restart the process using the most recent results from the "prnt ndx" replacing the '0' in line 30. !!!

The script needs improvement.  I would love to get rid of the pauses and having to restart. 

I also want to combine this script with the scripts that I use before and after into a single script.  The steps I follow and a link to the location of these scripts are as follows:

Kriging and zonal statistics  

1. Kriging of gridded climate data over multiple time steps.  Individual rasters are produced for each time step.
2. Zonal statistics of ALL rasters in the selected folder (produced by Step 1)
3.  Combine all .dbf files (zonal tables) in working directory (produced by Step 2) into a single table.
   • Columns = zones
   • Rows = time steps

[DuminduJayasekera]

Thanks Xander for sharing this. I have tried the code below but it does not gives me the zonal mean statistics for all the S_R_T in the County_Grid shape file. But it does provide for some of the grids. I have no clue why it does not produce the zonal mean for all the S_R_Ts. 

I appreciate if If some one can help me to figure this out.

Thanks a lot again.

import arcpy, os, arcinfo
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

env.workspace = "H:/outRaster/Rasters"

# Select the shapefile containing the polygons to use as boundaries
# for zonal statistics
watershedFeat = "H:/outRaster/Grid/County_Grid.shp"

# Select output folder for saving the output - zonal tables (.dbf files)
outDir = "H:/Table/"

x = arcpy.ListRasters()

for raster in arcpy.ListRasters(): #[0:]

raster_name = os.path.basename(raster).rstrip(os.path.splitext(raster)[1])
outTable = outDir + raster_name + "_TBL.dbf"
arcpy.gp.ZonalStatisticsAsTable(watershedFeat,"S_R_T", raster, outTable,"NODATA","MEAN")

arcpy.CheckInExtension("Spatial")

[DavidBetts]

Is there any pattern to which grids are left out?  Are there gaps in the sequence along "S_R_T"? or does the script move through the series and stop before finishing?  Does the attribute "S_R_T" have a unique value for each of the polygons in "County_Grid.shp"?

[DuminduJayasekera]

Thanks for your reply David. No there is no pattern. I checked it by doing manually suing the "zonal statistics as Table" tool but still produced the same output identical to the output by the script. Then I stored the shape file (County_Grid.shp) and the raster in the geodatabase and ran the script but no luck. The grids the shapefile are in 1 mile x 1 mile. Instead, as a trial, I used a county shapfile across the state of Kansas, USA and used the FIPS codes as the ID and it worked by producing zonal means for all the counties for the KS state. But, I really want to make it run for the County_Grid.shp in which grids are in 1 mile x 1 mile resolution to obtain the zonal means for all S_R_T. Any idea?

[DavidBetts]

What is the resolution of your raster data?  Sometimes when the polygons overlap only portions of the underlying raster cells, the zonal statistics isn't able to calculate an average value.  That might produce the odd pattern that doesn't look like a pattern, because it's only every so often that there is enough of an offset between County_Grid.shp and your raster cells that some of the 1x1mile grid cells only contains fractions of the underlying raster cells.  If that is the case, you might have to combine your grids into  2x2mile, 3x3mile or greater.  That's my best guess so far without seeing your raster data, or additional information.