Использование критерия для анализа стека растров в R

Question

Получил себе растровый стек со следующим описанием:

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

Каждый RasterLayer представляет собой массив аномалий геопотенциальной высоты со значениями в диапазоне от -300 до 300. Я получил 14 975 временных слоев, что означает каждый день с 1.1.1979 по 31.12.2019.

Одиночная точка сетки называется заблокированной, если аномалия достигает не менее +100 м в течение не менее 10 дней. Итак, я хотел бы написать функцию, которая возвращает запуски, заданные этим условием. Это означает, что я хочу вычислить каждую точку сетки в отдельности через временные слои.

Я не очень хорошо пишу функции, и эта проблема кажется мне слишком сложной для решения, но, возможно, кто-нибудь сможет мне помочь. Большое спасибо!

Answer 1

Начните с упрощения. Создайте вектор

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

Давайте найдем серии длиной 5 значений больше 10. Сначала найдем эти значения и установим для остальных 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

Теперь суммируйте, чтобы получить длину последовательностей. Это немного сложно, нам нужен cumsum, который сбрасывается на 0 (подробнее здесь)

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

И теперь мы можем посчитать количество запусков не менее 5

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

Как функция, которая могла бы быть

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

f(x, 10, 5)

Давайте жестко запрограммируем параметры, которые вы хотите использовать

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

И используйте эту функцию с растровыми данными.

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)