topic Re: Calculate zonal mean for multiple rasters using Python in Python Questions

Calculate zonal mean for multiple rasters using Python

DuminduJayasekera — Fri, 10 Dec 2021 21:29:49 GMT

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

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

Re: Calculate zonal mean for multiple rasters using Python

DanPatterson_Retired — Tue, 29 Aug 2017 01:31:55 GMT

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

Re: Calculate zonal mean for multiple rasters using Python

DuminduJayasekera — Tue, 29 Aug 2017 02:44:52 GMT

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

Re: Calculate zonal mean for multiple rasters using Python

DanPatterson_Retired — Tue, 29 Aug 2017 02:52:44 GMT

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

Re: Calculate zonal mean for multiple rasters using Python

DuminduJayasekera — Tue, 29 Aug 2017 03:12:33 GMT

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

Re: Calculate zonal mean for multiple rasters using Python

DanPatterson_Retired — Tue, 29 Aug 2017 03:15:29 GMT

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.

Re: Calculate zonal mean for multiple rasters using Python

DuminduJayasekera — Tue, 29 Aug 2017 03:43:52 GMT

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?

Re: Calculate zonal mean for multiple rasters using Python

DanPatterson_Retired — Sun, 12 Dec 2021 15:56:58 GMT

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

Re: Calculate zonal mean for multiple rasters using Python

DuminduJayasekera — Tue, 29 Aug 2017 15:16:46 GMT

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

Re: Calculate zonal mean for multiple rasters using Python

DanPatterson_Retired — Tue, 29 Aug 2017 15:42:25 GMT

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

Re: Calculate zonal mean for multiple rasters using Python

XanderBakker — Thu, 31 Aug 2017 12:45:27 GMT

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.