Building python script to compute weekly average from daily rainfall rasters

[BL]

Using:  ArcGIS 10 SP4 (ArcInfo) student evaluation with all extensions.

Goal: I need help finding a solution to integrate the cell statistics tool using Python to create a new weekly rainfall average layer for each week of 490 consecutive days of raster data.

Data: I have 490 raster files(.tif) of rainfall amounts for each day. Each raster was converted from a corresponding shapefile containing i) an even-spaced grid of points with no coincident points and ii) the same amount of grid points.

Example Raster for one day of data:
[ATTACH=CONFIG]13682[/ATTACH]

So far this is my logic for how to do this using Python:

0) Import arcpy, env from arcpy
1) Set workspace directory containing raster data for 490 consecutive days
2) Create list of rasters in workspace (rasterList = arcpy.ListRasters("*", "tif"))
3) Checkout spatial analyst extension

While rasterList is not empty
4) Create 7 variables for the first 7 consecutive days of rasters
5) Call Cell_Statistics function:  outCellStatistics = CellStatistics([raster1,..., raster7], "MEAN", "NODATA")   
6) Save new output raster:  outCellStatistics.save("C:/sapyexamples/output/week1")
7) Populate next consecutive 7 days into variables raster1,...,raster7
😎 Repeat steps 5 to 7 for 70 iterations (70 weeks)

I have tried messing around with the example in ArcGIS 10.0 Desktop Help. I am having an explicit problem with assigning the next 7 raster files into the 7 input variables for the next iteration. Any advice is appreciated.

[DarrenWiens2]

Create a counter = 1 before the loop. Enter loop. Is counter <= 7? Yes. Load raster into first variable. Add 1 to counter. Next.

Reset the counter and run your stats when it gets to 8.

[BL]

This is my code so far. I don't get any syntax or any errors for that matter, but get no output.

# Import system modules
import arcpy, os, sys
from arcpy import env
from arcpy.sa import *

# Set workspace
env.workspace = "D:/dailyrainfall/subset"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Create list of rasters
rasters = arcpy.ListRasters("*", "tif")
        
# Create Weekly rainfall raster
try:
        index = 70
        for index in rasters: 
      
                count = 0
        
                while count <=7:
                        # Set local variables
                        inRaster0 = rasters.Next()
                        count += 1
                        inRaster1 = rasters.Next()
                        count += 1
                        inRaster2 = rasters.Next()
                        count += 1
                        inRaster3 = rasters.Next()
                        count += 1
                        inRaster4 = rasters.Next()
                        count += 1
                        inRaster5 = rasters.Next()
                        count += 1
                        inRaster6 = rasters.Next()
                        count += 1
                index -= 1

                # Execute CellStatistics
                outCellStatistics = CellStatistics([inRaster0, inRaster1, inRaster2, inRaster3, inRaster4, inRaster5, inRaster6], "SUM", "NODATA")

                output = outCellStatistics.name

                # Save the output 
                outCellStatistics.save("D:/dailyrainfall/subset/result/" + output + ".tif")
except:
    # If an error occurred while running a tool, then print the messages.
    print arcpy.GetMessages()

[MarcinGasior]

Here is a code which I tested for 21 rasters - seemed to work well:

try:
    # Create list of rasters
    rasters = arcpy.ListRasters("*", "tif")
    rastersCount = len(rasters)

    counter = 0
    weekNum = 1

    while counter < rastersCount:
        #collect 7 consecutive rasters
        weekRasters = []
        for i in range(counter,(counter+7)):
            weekRasters.append(rasters)

        # Execute CellStatistics (one can use list of rasters directly here)
        outCellStatistics = CellStatistics(weekRasters, "MEAN", "NODATA")
        # Save the output
        outRasterName = "week_"+ str(weekNum)
        outCellStatistics.save("C:/tmp/PrecipRasters/" + outRasterName )

        counter += 7
        weekNum += 1

[BL]

m.gasior, that code you posted solved my problems.

dkwiens, thanks for the conceptual advice.

[Leo_KrisPalao]

Hi m.gasior,

I tried to implement the python script that you suggested but I get an error when I ran it in in PythonWin and in ArcGIS 10. I have 366 rasters in TIFF format and I want to calculate the mean values every 10 days. I tried to modify the script but get an error. This is the modified script that I am trying to use. Could you help me revise my script?

Thanks,
-Leo

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
env.workspace = "E:\PCP_TRMM_1998-2013"
rasters = arcpy.ListRasters("*", "tif")
rastersCount = len(rasters)
counter = 0
weekDek = 1
while counter < rastersCount:
    dekadRasters = []
    for i in range(counter,(counter+10)):
        dekRasters.append(rasters)
        outCellStatistics = CellStatistics(dekRasters, "MEAN", "NODATA")
        outRasterName = "dekad_"+ str(weekDek)
        outCellStatistics.save("E:\PCP_TRMM_1998-2013\Dekads" + outRasterName )
        counter += 10
        weekDek += 1

[ClintDow]

Hi m.gasior,

I tried to implement the python script that you suggested but I get an error when I ran it in in PythonWin and in ArcGIS 10. I have 366 rasters in TIFF format and I want to calculate the mean values every 10 days. I tried to modify the script but get an error. This is the modified script that I am trying to use. Could you help me revise my script?

Thanks,
-Leo

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
env.workspace = "E:\PCP_TRMM_1998-2013"
rasters = arcpy.ListRasters("*", "tif")
rastersCount = len(rasters)
counter = 0
weekDek = 1
while counter < rastersCount:
 dekadRasters = []
 for i in range(counter,(counter+10)):
  dekRasters.append(rasters)
  
 outCellStatistics = CellStatistics(dekRasters, "MEAN", "NODATA")
 outRasterName = "dekad_"+ str(weekDek)
 outCellStatistics.save("E:\PCP_TRMM_1998-2013\Dekads" + outRasterName )
 counter += 10
 weekDek += 1

Added quote tags for readability.

You assign your list variable the name dekadRasters then try to append to dekRasters, you will need to fix that variable name on lines 12 & 14

Also your path strings have two issues, first the '\' characters need to either be replaced with '/' or use a raw string. Otherwise the backslash will produce unexpected results in the string since it is the escape character. See here for more info on that concept: Python Docs: String Literals

Lastly you might want to revisit your loop, since counter is increasing by 10 each iteration and raster count may not be a multiple of 10 you will get an IndexError when counter+10 exceeds the total number of rasters.

[Leo_KrisPalao]

Hi Clint,

I tied to modify the script but still I get error. I also reduced my rasters to 360 so my counter will not exceed the total number of rasters. Below is my modified script.

Thanks,
-Leo

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
env.workspace = "E:/PCP_TRMM_1998-2013/Dekads"   # character "\" is already replaced with "/"
rasters = arcpy.ListRasters("*", "tif")
rastersCount = len(rasters)
counter = 0
dek = 1
while counter < rastersCount:
    dekRasters = []                      # I already checked the variables. This is already the same with line 12 and 13
    for i in range(counter,(counter+10)):
        dekRasters.append(rasters)
        outCellStatistics = CellStatistics(dekRasters, "MEAN", "NODATA")
        outRasterName = "dekad_" + str(dek)
        outCellStatistics.save("E:/PCP_TRMM_1998-2013/Dekads" + outRasterName )  # character "\" is already replaced with "/"
        counter += 10
        dek += 1

By the way, this is the error I am getting when I run it in PythonWin:

PythonWin 2.6.5 (r265:79096, Mar 19 2010, 21:48:26) [MSC v.1500 32 bit (Intel)] on win32.
Portions Copyright 1994-2008 Mark Hammond - see 'Help/About PythonWin' for further copyright information.
Traceback (most recent call last):
  File "C:\Python26\ArcGIS10.0\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 312, in RunScript
    exec codeObject in __main__.__dict__
  File "E:\Python_scripts\ArcGIS\cellstatistics.py", line 13, in <module>
    outCellStatistics = CellStatistics(dekRasters, "MEAN", "NODATA")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 2780, in CellStatistics
    ignore_nodata)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 2776, in wrapper
    return _wrapLocalFunctionRaster(u"CellStatistics_sa", [function] + Utils.flattenLists(in_rasters_or_constants))
RuntimeError: ERROR 000824: The tool is not licensed.

[ClintDow]

The error 'RuntimeError: ERROR 000824: The tool is not licensed.' Indicates you have an issue with your licenses. CellStatistics uses the spatial analyst license.
Its probably easier to check out the license in the script then to ensure the user has it enabled similar to how tlupher has it in their code snippet.

if arcpy.CheckExtension("Spatial") == "Available":
  arcpy.CheckOutExtension("Spatial")

[Leo_KrisPalao]

Hi Clint,

I ran the code successfully, but I have issues with the results:

1) I only have 360 rasters so I expected my output to be 36 rasters only (360/10). But in the output it generated 40 rasters. I entered in the interactive window the following commands just to double check the results:

>>> Print rastersCount
360; this is correct
>>> Print counter
400; this seems not correct. Would this value only have to be 360?
>>> Print dek
41; this seems not correct. I think this value would only have to be 36
>>> Print outRasterName
040; this is also not correct. I am only expecting to have 036.

2) I compared the results from my manual computation and the automated computation, and the results were different:
Values generated for manual computation (max: 37.63888931; min: 0; mean: 2.4461317333); I think these are the right values
Values generated by the python script (max: 44.71593856811523, min: 0, but a weird mean value of -1.#QNAN)

This is the code that I used:

import arcpy, os, sys
if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.CheckOutExtension("Spatial")
from arcpy import env
from arcpy.sa import *
env.workspace = "E:/PCP_TRMM_1998-2013_dekad"
rasters = arcpy.ListRasters("*", "tif")
rastersCount = len(rasters)
counter = 0
dek = 1
while counter < rastersCount:
    dekad = []
    for i in range(counter,(counter+10)):
        dekad.append(rasters)
        outCellStatistics = CellStatistics(dekad, "MEAN", "NODATA")
        outRasterName = "0" + str(dek)
        outCellStatistics.save("E:/Dekads/Dek" + outRasterName)
        counter += 10
        dek += 1


By the way Clint, if you are kind enough, I would like to save my output into a GeoTIFF file. So I modified my code like this:
          outCellStatistics.save("E:/Dekads/Dek.tif" + outRasterName)

but in the error message it says: File "E:\Python_scripts\ArcGIS\cellstatistics.py", line 17, in <module>
    outCellStatistics.save("E:/Dekads/Dek_.tif" + outRasterName)
RuntimeError: ERROR 010093: Output raster format UNKNOWN is unsupported.

Thanks for your kind assistance.

-Leo