nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Building python script to compute weekly average from daily rainfall rasters

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.

example_1.png

Vind meer berichten getagd met

Accepted answers

ClintDow

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.

All comments

Leo_KrisPalao

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

Thanks again for the help.

-Leo

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

ClintDow

rastersCount = 360
counter = 0
dek = 1

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

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

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

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

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

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

dkwiens, thanks for the conceptual advice.

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

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

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.

Doe mee aan het gesprek

Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.

Snelle Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.