Using raster stack function time in R?

Question

I successfully ran a function over a smaller test raster stack (48 layers) in about 30 minutes but cannot run it over my large stack (1462 layers). The last one I tried over the large stack ran from 9PM to 7AM and had not finished processing. I stopped it and tried a simpler function and I have not been able to get it to finish yet.

Large stack information: layers = 1462, total elements = 1275602310, 16.4 mb, ncols per layer = 1405, nrows per layer = 621.

# Import libraries ----
library(tidyverse)
library(dplyr)
library(magrittr)
library(raster)

# Call the source file for the function
source("analysis_functions.R")

# Identify file path for files to stack
t_dsn <- "filepath_here"
p_dsn <- "filepath_here"

# List files to stack
t_list <- list.files(t_dsn, pattern = "tif$", full.names = TRUE)
p_list <- list.files(t_dsn, pattern = "tif$", full.names = TRUE)

# Create stacks
t<- stack(t_list)
p<- stack(t_list)

# Calculate a pixel-wise mean value through the stack
tm <-calc(t, mean)
pm <-calc(p, mean)

# Calculate a pixel-wise Nash Sutcliffe Efficiency value through the stack
nse_t <- myNSE(sim = t, obs = p, obs_mean = mp)

#Lastly, here is the function as defined in the source file
myNSE <- function(sim, obs, obs_mean){
  nse <- 1 - ((sum((obs-sim)^2)) / (sum((obs-obs_mean)^2)))
}

The test stack and the large stack are comprised of the same exact files, just a smaller set.

The mean calculation for the test stack takes about 2 minutes.

The mean calculation for the large stack takes about 30 minutes.

The NSE calculation for the test stack takes about 2 minutes.

The NSE calculation in the large stack hangs for more than 9 hours and I have not been able to complete.

Lastly, I have successfully run all commands until I get to the NSE one. This is the one that seems to lag.

Any suggestions?

Answer 1

Couple of problems I can see. One, I think you've got a typo - I'm assuming pm and mp are the same object?

Two, your NSE function assigns the calculation to a variable but then doesn't return it. Get rid of nse <-.

Since you've mentioned that all of the rasters involved can stack on each other, you can rewrite your NSE function to take advantage of that. So for instance,

all_dat <- stack(pm, t, p)

Now each pixel in all_dat is a vector of numbers, call it cell. cell[1] is the mean of t, cell[2:(nlayers(t) + 1)] are your obs data and the rest are your model predictions. So the NSE function written to operate on that stack cell-wise should look something like

cell_NSE <- function(cell) {
  vec_center <- (length(cell) - 1) / 2
  obs_mean <- cell[1]
  obs <- cell[c(2:(vec_center + 1))] 
  sim <- cell[(vec_center + 2):length(cell)]

  1 - ((sum((obs - sim)^2)) / (sum((obs - obs_mean)^2)))
}

and your operation would be

out <- calc(all_dat, cell_NSE)

To speed it up further, clusterR may be of use:

beginCluster(max(1, (parallel::detectCores() - 1))) # always leave a spare
out_fast <- clusterR(all_dat, fun  = calc, args = list(fun = cell_NSE))
endCluster()