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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

[DanPatterson_Retired]

There is a disconnect between your error message and the script your posted.

Where is this line?

outZonalStatistics = ZonalStatisticsAsTable......

check this link

the code samples show that you indeed need an assignment And you should create your zone file as an integer raster so that you have explicit control over its cell size and extent

[DuminduJayasekera]

Thanks Dan. I have updated the code (Please see the code above) but it is still not doing the job I want. Any idea?

[DanPatterson_Retired]

Dumindu, I think you are going to have to throw some print statements in there to ensure that the parameters are indeed correct.  And the big one, is the zonal statistics is looking for a polygon file or its raster equivalent

[DuminduJayasekera]

OK. No error message this time but i cannot see the output. 

[DanPatterson_Retired]

did it get added to the table of contents?  If not, manually add it... better still manually add it to a new dataframe incase the file doesn't have a defined coordinate system and therefore won't play nice with all the other data.  If there is nothing to see, then it made it through the process but didn't work... and I would like to see your 'zone' file to make sure it wasn't a point file.

[DuminduJayasekera]

I just checked the shape file and I have used a point shape file. I have corrected to the polygon shape file. (Please see the updated code) But, the code above does not give an error and shows it completed the run. But, I cannot see the output. What am I doing wrong here?

[DanPatterson_Retired]

Did you do as I suggested and add the result to a new dataframe?  if nothing appeared, use File Explorer and report the file size and save it to the following location and name... there is little point in repeating what isn't working... try a new location, new raster type and see if anything is created... if it works, then figure out how to get what where you want it.

outZonalStatistics.save("C:/temp/x.tif")

[DuminduJayasekera]

I ran the code upto inPointFeatures line in Python window in ArcGIS and gave this error below. Seems to be the code is correct. Cannot understand why.

Runtime error  Traceback (most recent call last):   File "<string>", line 13, in <module>   File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 6477, in MakeFeatureLayer     raise e ExecuteError: ERROR 000622: Failed to execute (Make Feature Layer). Parameters are not valid. ERROR 000628: Cannot set input into parameter in_features.‍‍

[DanPatterson_Retired]

one of these two is wrong.... do the process manually, go to the Results window and copy the python syntax.... ps... your output filename PL..... shouldn't be that long, you might be asking for trouble when you get into long filenames, pathnames and older versions of python

"Public_Land_Survey_System", r"C:/Users/Desktop/KS_All.gdb/PLSS_KS_All_WeeklyPr_PRISM_800m_2016_new")