Using a criterion to analyze rasterstack in R

Question

I got myself a rasterstack with following description:

 class      : RasterStack 
 dimensions : 221, 121, 26741, 14975  (nrow, ncol, ncell, nlayers)
 resolution : 0.25, 0.25  (x, y)
 extent     : 14.875, 45.125, 24.875, 80.125  (xmin, xmax, ymin, ymax)
 crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
 names      :       layer.1,       layer.2,       layer.3,       layer.4,       layer.5,       
 layer.6,       layer.7,       layer.8,       layer.9,      layer.10,      layer.11,      layer.12,      
 layer.13,      layer.14,      layer.15, ... 
 min values : -291.44314575, -329.02883911, -388.52404785, -377.29467773, -358.85354614, 
 -342.63455200, -268.73141479, -176.84980774, -316.51959229, -267.97073364, -280.14392090, 
 -313.80172729, -287.72854614, -288.26785278, -279.92047119, ... 
 max values :   139.2527466,   101.1818314,   122.2153473,   101.3840714,   163.6441345,   
 197.5134430,   215.1998291,   200.4686127,   159.4233856,   108.7927933,   148.2487488,   
 167.5716858,   198.3386993,   296.5519714,   276.6212463, ...

Every RasterLayer is a geopotential height anomaly array with values ranging from some -300 to 300. I got 14 975 time layers, which mean every day since 1.1.1979 till 31.12.2019.

A single grid point is called blocked, if anomaly reaches at least +100 m for at least 10 days. So I would like to write a function, which returns runs given by this condition. It means that I wanna calculate every single grid point at its own through time layers

I am not very good at writing functions and this problem seems too tough for me to solve, but maybe someone could help me. Thank you very much!

Robert Hijmans · Answer 1 · 2020-10-31T20:48:47.930

Start by simplifying. Create a vector

x <- c(1:20, 1:20)
x
# [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20  1  2  3  4  5
#[26]  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

Let's find runs of length 5 of values larger than 10. First find these values, and set the others to NA

y <- x > 10
y
# [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
#[13]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE
#[25] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
#[37]  TRUE  TRUE  TRUE  TRUE

Now cumulate to get the length of the sequences. That is a bit tricky, we need a cumsum that resets at 0 (more discussion here)

n <- ave(y, cumsum(y == 0), FUN = cumsum)
n
# [1]  0  0  0  0  0  0  0  0  0  0  1  2  3  4  5  6  7  8  9 10  0  0  0  0  0
#[26]  0  0  0  0  0  1  2  3  4  5  6  7  8  9 10

And now we can count the number of runs of at least 5

sum(n==5)
#[1] 2

As a function that could be

f <- function(x, a, b) {
    y <- x > a
    n <- ave(y, cumsum(y == 0), FUN = cumsum)
    sum(n == b)
}

f(x, 10, 5)

Let's hard-code the parameters that you want to use

ff <- function(x) {
    y <- x >= 100
    n <- ave(y, cumsum(y == 0), FUN = cumsum)
    sum(n == 10)
}

And use this function with raster data.

library(raster)
b <- brick(ncol=10, nrow=10, nl=100)
set.seed(2)     
values(b) = sample(90:120, replace=T, ncell(b)*nlayers(b))

r <- calc(b, ff)
r 
#class      : RasterLayer 
#dimensions : 10, 10, 100  (nrow, ncol, ncell)
#resolution : 36, 18  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer 
#values     : 0, 3  (min, max)

Using a criterion to analyze rasterstack in R

1 Answers1