Search code examples
rgisraster

Cumulative sum of multiple rasters in R


I have 200 rasters with air temperatures for each day stored in a directory and I would like to do a cumulative sum of temperatures for each day, e.g. day1 = day1; day2 = day1 + day2; day3 = day1 + day2 + day3, etc. I tried to do it in raster package, but the cumsum function calculates the cumulative sum of all cells in each raster and not the cumulative sum of individual rasters (at least it seems to me from my results). I tried to do it like this:

library(raster)

setwd("C:/Air_temperatures/AT") # set working directory

AT.files <- list.files(pattern="AT") # list files with name AT
All.AT <- unique(AT.files) # unique files
  for (i in 1:200) {
    AT.i <- All.AT[i] # make i
    AT.raster.i<-raster(AT.i) # make rasters from files
    AT.sum.i <- cumsum(AT.raster.i) # ???calculate cumulative sum???
    writeRaster(AT.sum.i,
               paste(AT.i, "_cumsum.tif", sep = ""),
              datatype='FLT4S', overwrite=TRUE) # write new files and paste new name
  } 

This loop worked when I tried for example to add a constant to each raster and so on, but I have no idea how to calculate the cumulative sum.


Solution

  • The Reduce function might help here. This function successively combines objects based on a specified function. Since there is a + method for Raster* objects, you can use this as the function and obtain a cumulative sum of the raster values.

    First, create the 200 raster objects in a list:

    theATs <- lapply(All.AT, raster)
    

    Then, use the Reduce function with the accumulate argument set to TRUE

    res <- Reduce("+", theATs, accumulate = TRUE)
    

    Then write the results to files:

    lapply(seq_along(res), function(x) {
        writeRaster(res[[x]],
    # Apparently, using .tif files threw an error
    #     paste(All.AT[x], "_cumsum.tif", sep = ""),
    # Alternatively, omit file extension to write default type (usually .grd)
          paste(All.AT[x], "_cumsum", sep = ""),
          datatype = 'FLT4S', overwrite = TRUE)
        })
    

    The rasters you create here are stored in memory, so if the original files are very large, you may run into problems.