CODEWITHSUNDEEP
rraster
Kongsol
Kongsol

Reputation: 37

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.

Upvotes: 2

Views: 1388

Answers (1)

DJack
DJack

Reputation: 4940

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)

Upvotes: 2

Related Questions