0

我得到了一个带有以下描述的光栅堆栈：

 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 m 并持续至少 10 天，则将单个网格点称为阻塞。所以我想写一个函数，它返回这个条件给出的运行。这意味着我想通过时间层自己计算每个网格点

我不太擅长编写函数，这个问题对我来说似乎太难解决了，但也许有人可以帮助我。非常感谢！

4

1 回答 1

0

从简化开始。创建矢量

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

让我们找到大于 10 的长度为 5 的游程。首先找到这些值，然后将其他值设置为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

现在累积以获得序列的长度。这有点棘手，我们需要一个重置为 0 的 cumsum（更多讨论在这里）

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)
于 2020-10-31T20:06:29.887 回答