1

I have a raster brick b of global fAPAR data which is 35 years’ worth of daily data (i.e. large).

The Aim

Is there a way I can loop through this brick and run a function on 365 layers (days) at a time? The end goal is to run a function to identify the layer (date) in which the value of each gridcell goes above a threshold of 20% of the max value for that gridcell for each year.

So the code would need to:

  1. Work out max value of each gridcell in a year
  2. Return the layer (date) on which each gridcell crosses 20% of that max value

Code

The following code finds the max for each year, but doesn't work on my data due to its size (daily, global, 0.5 deg. data). The key may be to work with 365 layers at a time, not the whole brick (if possible).

# Example 10-year brick:
b <- brick(nrow=108,ncol=132,nl=3650)
values(b) <- runif(ncell(b)*nlayers(b))

i <- rep(1:nlayers(b), each=365)
# if incomplete last set of days:  
i <- i[1:nlayers(b)]
x <- stackApply(b, i, max)

My goal is quite similar to this and this. I also wonder whether I need to divide my data up into individual years and run something on multiple files, given the size issues.

Konrad Rudolph
  • 530,221
  • 131
  • 937
  • 1,214
Ndharwood
  • 123
  • 3
  • 11

1 Answers1

2

This is probably because raster assess the memory needs wrong. You can set lower memory limits with rasterOptions.

You could loop over years (I note that you think they all have 365 days) with i

yrs <- unique(i)
out <- list()
for (j in yrs) {
  bb <- b[[which(i==j)]]
  out[[j]] <- max(bb) 
}
out <- stack(out)

You can also try to do this with terra (the new package replacing raster)

library(terra)
x <- rast(b)
z <- tapp(x, i, max)

this is not well defined, I think:

Return the layer (date) on which each gridcell crosses 20% of that max value

As there can be many such dates. Below, out1 gets them all and out1 gets the first occurrence:

library(terra)
x <- rast(nrow=108,ncol=132, nl=3650)
set.seed(1)
values(x) <- runif(size(x))
i <- rep(1:nlayers(b), each=365)
i <- i[1:nlayers(b)]

zz <- tapp(x, i, function(i) 0.2 * max(i))

yrs <- unique(i)
out1 <- out2 <- list()
for (j in yrs) {
    xx <- x[[which(i==j)]]
    rr <- zz[[j]] < xx
    out1[[j]] <- rr
    out2[[j]] <- which.lyr(rr)
}

out1 <- rast(out1)
out2 <- rast(out2)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Excellent - thanks Robert, both work on my data. Did you have any ideas for part 2 under 'Aim'? Would I be able to use quantile(0.2) or similar? – Ndharwood Sep 10 '21 at 15:27
  • I have expanded the answer – Robert Hijmans Sep 13 '21 at 07:49
  • Thanks a million, Robert! Is there a reason I'm getting 'Error: cannot allocate vector of size..' when I run `zz <- tapp(x, i, function(i) 0.2 * max(i))` on my dataset? The `tapp(x, i, max)` code you suggested above ran on my data, so it's weird that adding a function in causes it to run out of memory. (R version 4.1.0, terra v 1.3-22) – Ndharwood Sep 16 '21 at 13:10