Making cumulative sum of daily precipitation raster in ArcGIS Desktop?

Question

I have daily precipitation rasters from 1981 to 2016 and i need to get the total sum of each pixel in the raster and the mean of them.

How can i iterate through them to sum all of them (more than 13000 rasters)?

I've tried Cell statistics with the iterate rasters loop in arcmap but it does not work.

Answer 1

Try creating a python script similar to the below... from ListRasters and Cell Statistics

import arcpy

# Set the current workspace
arcpy.env.workspace = "c:/data/yourRasters"

# Get and print a list of GRIDs from the workspace
rasters = arcpy.ListRasters("*", "GRID")

# Run Cell Statistics arcpy tool
outCellStats = CellStatistics(rasters, "SUM", "DATA")
outCellStats.save("c:/data/yourSummedRasters.tif")

Answer 2

I would try using the regular Cell Statistics tool. If you put all your precipitation rasters in the same directory, you can select them all at once as the input for Cell Statistics and calculate SUM and MEAN. No need to use an iterator in ModelBuilder.

Answer 3

You can do it somewhat easily in GDAL. Might need to install the Python bindings for it if you don't have it already. Didn't test this code, as I don't have data ready to test, but it'll probably work. It'll be slow though, as I had no access to your metadata nor your setup to optimize it. I'd first give it a try with, say, 5 rasters, and see if it works. If it does, do with all 13k+.

This code expects all rasters to be in the same folder. It also expects them to be of the same file type (e.g. all GeoTIFFs), same extent, and to align properly. It also makes no evaluation of nodata values, though that should be easy enough to implement if needed. You'll need to change the "name_of_driver" to the GDAL code of your raster type (e.g. "GTiff" for GeoTIFFs). You'll also need to change the path to your rasters folder, and where it reads '.ext' (line 11 and 17) change it to the appropriate filetype extension (e.g. '.tif' for GeoTIFFs).

import os
import numpy as np
from osgeo import gdal
from osgeo.gdalconst import *
from osgeo.gdal_array import *

driver = gdal.GetDriverByName("name_of_driver")
driver.Register()

os.chdir(r'path/to/rasters/folder')
li_rasters = [raster for raster in os.listdir(os.getcwd()) if os.path.splitext(raster)[-1] == '.ext']

raster = gdal.Open(li_rasters[0], GA_ReadOnly)
x, y = raster.RasterXSize, raster.RasterYSize
band = raster.GetRasterBand(1)

final_raster = driver.Create('final_raster.ext', x, y)
final_band = final_raster.GetRasterBand(1)
final_data = np.zeros((x, y), dtype = np.dtype(GDALTypeCodeToNumericTypeCode(band.DataType)))

for i in range(len(li_rasters) - 1):
    rast1 = gdal.Open(li_rasters[i], GA_ReadOnly)
    band1 = rast1.GetRasterBand(1)
    rast2 = gdal.Open(li_rasters[i + 1], GA_ReadOnly)
    band2 = rast2.GetRasterBand(1)
    for j in range(len(x)):
        for k in range(len(y)):
            data1 = band1.ReadAsArray(k, j, 1, 1)
            data2 = band2.ReadAsArray(k, j, 1, 1)
            final_data[k, j] += (data1 + data2)

final_band.WriteArray(final_data)
final_band.FlushCache()
final_band.GetStatistics(0, 1)