Function to sum each grid cells of raster stack using other rasters as an indicator

## input raster
s <- stack(list.files("~/dailyraster", full.names=TRUE)) # daily raster stack
r_start <- raster("~/stackSumSTART.asc") # this raster contain starting Julian day
r_end <- raster("~/stackSumEND.asc") # this raster contain ending Julian day
noNAcells <- which(!is.na(r[])) # cell numbers which contain values

## dummy raster
x <- r
x[] <- NA 

## loop
for (i in noNAcells) {      
  x[i] <- sum(s[[r_start[i]:r_end[i]]][i])
}

I would like to create a function like stackApply(), but I want it to work on a cell basis.

Above is a for() loop version and it works well, but it takes too much time.

The point is that each cell gets the range of sum() from two raster layers, r_start, r_end in above script.

Now I am struggling to transform this code using apply() family.

Is there any possibility to improve the speed with for() loop? or please give me some tips to write this code in apply()

Any comments will help me, thank you.

1 answer

  • answered 2018-04-17 05:17 DJack

    Your approach

    x <- s$layer.1
    
    system.time(
    for (i in 1:ncell(x)) {
         x[i] <- sum(s[[r_start[i]:r_end[i]]][i], na.rm = T)
       }
    )
    
       user  system elapsed 
      0.708   0.000   0.710
    

    My proposal

    You can add the rasters used as indices at the end of your stack and then use calc to highly speed up the process (~30-50x).

    s2 <- stack(s, r_start, r_end)
    sum_time <- function(x) {sum(x[x[6]:x[7]], na.rm = T)}
    
    system.time(
       output <- calc(s2, fun = sum_time)
    )
    
       user  system elapsed 
      0.016   0.000   0.015
    
    all.equal(x, output)
    
    [1] TRUE
    

    Sample Data

    library(raster)
    
    # Generate rasters of random values
    r1 <- r2 <- r3 <- r4 <- r5 <- r_start <- r_end <- raster(ncol=10, nrow=10)
    
    r1[] <- rnorm(ncell(r1), 1, 0.2)
    r2[] <- rnorm(ncell(r2), 1, 0.2)
    r3[] <- rnorm(ncell(r3), 1, 0.2)
    r4[] <- rnorm(ncell(r4), 1, 0.2)
    r5[] <- rnorm(ncell(r5), 1, 0.2)
    s <- stack(r1,r2,r3,r4,r5)
    
    r_start[] <- sample(1:2, ncell(r_start),replace = T)
    r_end[] <- sample(3:5, ncell(r_end),replace = T)