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.

[Leo_KrisPalao]

Hi Clint and other users,

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

Thanks for any help,
-Leo

[ClintDow · Accepted solution]

I must have missed the notification for your first reply.

Those results don't seem possible based on the code, I'm a bit confused. The counter value shouldn't be able to be higher than the rastersCount.

I simplified the loop to make sure something wasn't missing to:

rastersCount = 360
counter = 0
dek = 1

while counter < rastersCount:
 for i in range(counter,(counter+10)):
  print i
 counter += 10
 dek += 1

And as expected it prints from 0 - 359. I unfortunately don't have a lot of time today, but I would maybe start by manually setting the rastersCount to 360 and add a few print statements into the loops to observe the values changing with each iteration to begin troubleshooting.

If possible when posting code please wrap it in CODE tags (click the # symbol in the post formatting to add them) so it retains your whitespace, its easier to figure out where loops start and end that way.

In regards to the question about the line: outCellStatistics.save("E:/Dekads/Dek.tif" + outRasterName)
The .tif extension is in the wrong spot, the string should be "E:/Dekads/Dek" + outRasterName +".tif" so the extension is appended to the end of the string rather than in the middle.

[Leo_KrisPalao]

Hi Clint,

It's been a while since my last inquiry and apologies if I have not answered your last reply. the code that you suggested is perfectly running. I tested it using 100 rasters.

However, I want to expand the code a little bit. I have a total of 366 rasters representing day of year. My objective is how to run the script that it will compute every 10 days from raster 1 to raster 360, but on the last step, it will only compute for the remaining 6 rasters. I know this can be done but I cannot imagine how to proceed.

Any help will be much appreciated. In the meantime I will try to research for some examples in the internet.

Thanks,

-Leo

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

# workspace
env.workspace = 'path/to/input/rasters'
rasters = arcpy.ListRasters()
out_ws = 'path/to/output/directory'

rastersCount = len(rasters)
counter = 0
dek = 1

# compute dekad rasters
while counter < rastersCount:
    dekad = []
    for i in range(counter,(counter+10)):
        print i
        dekad.append(rasters)
        outCellStatistics = CellStatistics(dekad, "MEAN", "NODATA")
        outRasterName = out_ws + 'tmax_dek_{:03d}.tif'.format (dek)
        #outRasterName = "00" + str(dek) + ".tif"
        outCellStatistics.save(outRasterName)
    counter += 10
    dek += 1

[Luke_Pinner]

Another way of looping is to use the range function with a step of 7 and list slice notation:

rasterList = arcpy.ListRasters("*", "tif")
rasterCount = len(rasterList)
#loop from 0 to rasterCount incrementing by 7 each iteration
for i in range(0, rasterCount, 7): 
    # Slice 7 elements from the list
    outCellStatistics = CellStatistics(rasterList, "MEAN", "NODATA") 
    outCellStatistics.save("C:/sapyexamples/output/week{}".format(i/7)) #e.g. week0 is first week. For week1 as first, use i/7+1

[Leo_KrisPalao]

Hi Luke, thanks. This code is much simpler and very much readable.

Thanks again for the help.

-Leo