Calculate zonal mean for multiple rasters using Python

[DuminduJayasekera]

Hi,

I need to calculate the zonal mean (areal precipitation) for grids in a shape file for multiple rasters using a Python script. I have created a script but it gives me an error. 

Any help is appreciated to fix the issue. Please see the code below.

import arcpy
from arcpy import env
from arcpy.sa import *
import os
import arcgisscripting, sys
arcpy.env.overwriteOutput = True
# Set environment settings
env.workspace = "C:/Users/Desktop/KS_All.gdb"
# Make a layer from the feature class 
arcpy.MakeFeatureLayer_management("Public_Land_Survey_System", "lyr") ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

# Write the selected features to a new featureclass
arcpy.CopyFeatures_management("Public_Land_Survey_System", r"C:/Users/Desktop/KS_All.gdb/PLSS_KS_All_WeeklyPr_PRISM_800m_2016_new")

# Set environment settings
arcpy.env.workspace = "H:/PRISM_800m_Daily"

# Set local variables
inPointFeatures = r"C:/Users/Desktop/KS_All.gdb/PLSS_KS_All_WeeklyPr_PRISM_800m_2016_new.shp"

inRasterList = [["Week_1_Avg2016_10.tif", "Oct_PCPwk1"],
["Week_2_Avg2016_10.tif", "Oct_PCPwk2"],
["Week_3_Avg2016_10.tif", "Oct_PCPwk3"],
["Week_4_Avg2016_10.tif", "Oct_PCPwk4"],
["Week_5_Avg2016_10.tif", "Oct_PCPwk5"],
["Week_1_Avg2016_11.tif", "Nov_PCPwk1"],
["Week_2_Avg2016_11.tif", "Nov_PCPwk2"]]

arcpy.CheckOutExtension("Spatial")

zoneField = "S_R_T"
outZonalStatistics = ZonalStatistics(inPointFeatures, zoneField, inRasterList, "MEAN", "DATA")

# Save the output 
outZonalStatistics.save("C:/Users/Desktop/KS_All.gdb/zonalstattblout")

fieldName = "YEAR"
expression = "getClass(!YEAR!)"
codeblock = """def getClass(YEAR):
  if YEAR >= 1998:
        return 2016
  else:
        return -9999"""

# Execute CalculateField 
arcpy.CalculateField_management(inPointFeatures, fieldName, expression, "PYTHON_9.3", codeblock)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[XanderBakker]

Not sure if the code represents the latest version of what you are trying to run, but I noticed that you used the following:

inPointFeatures = r"C:/Users/Desktop/KS_All.gdb/PLSS_KS_All_WeeklyPr_PRISM_800m_2016_new.shp"‍

In this case you are trying the write a shapefile in the "folder" that is a file geodatabase. You should loose the ".shp" at the end of the statement.  

If possible you could share a sample of your data to see if the error can be reproduced and to correct any errors in the code.