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

# Call the source file for the function

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

1 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))
  • TY, I'll give it a go. Feb 2, 2019 at 18:00

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.