4

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.

4

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

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.

enter image description here

  • But i have 13.149 rasters. I can't even load them to arcmap without having it crashing down. Is there not a way to accomplish this without importing them? – Federico Gómez Jan 30 '18 at 15:36
  • 1
    If arc crashes when you load so many rasters, try creating a python script to gather all the layers and sum them using the arcpy cell statistics function. – mikeLdub Jan 30 '18 at 16:03
  • @mikeLdub but how can i make that? can you please give me an example of how gather all the layers in the script? – Federico Gómez Jan 30 '18 at 16:06
  • I agree with @mikeLdub. If ArcMap won't accept that many rasters as input for Cell Statistics, a script is going to be the way to go. – lambertj Jan 30 '18 at 16:20
  • 1
    I've in the past had to add thousands of raster together, I found running cell statistics in batches of 200 would work then add the batches. Throwing all 13,000+ rasters at the tool will blow it up. – Hornbydd Jan 30 '18 at 16:30
0

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)
  • there is a mistake in this code: TypeError: object of type 'int' has no len() – Kais Tounsi Sep 15 at 14:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.