根据另一个栅格对 Terra 中的大型栅格堆栈进行重新分类

Question

我正在尝试执行以下操作：

我有一堆每月温度栅格，并据此计算堆栈的第 90 个分位数。

然后，我想计算每个栅格单元等于或高于该分位数的频率。这样，我最终会得到一个栅格，其中包含超过第 90 分位数的频率的总和值。我知道这似乎没有必要，因为我采用的是第 90 分位数，但它会使对一堆投影（未来）栅格重复该过程变得更容易，因此我可以查看超过频率的变化过去与现在之间的门槛。

所以我从以下开始： 然后取第 90 分位数：

然后，我尝试查看 hist_temp_stack 的每一层，并希望对于栅格中的每个像元，该层超过 T90 值的计数：

  hist_extreme_temp <- app(hist_temp_stack, fun = function(x) ifelse(x >= T90[], 1, 0))
 
  hist_extreme_temp_count <- sum(hist_extreme_temp)

导致

There were 50 or more warnings (use warnings() to see the first 50)

>   warnings()
Warning messages:
1: In x >= T90[] :  longer object length is not a multiple of shorter object length

在所有这些警告之后......

>   hist_extreme_temp_count
class       : SpatRaster 
dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +no_defs 
source      : spat_362c1f613e88_13868.tif 
name        :   sum 
min value   :     4 
max value   : 62977 

但是这个最大值没有任何意义。 最大值不应超过 492（hist_temp_stack 中的层数），因为我正在尝试对每个单元格 hist_temp_stack 中超过第 90 分位数 (T90) 的层数求和。

为了避免与 hist_temp_stack 长于 T90 相关的警告，我还执行了以下操作，但这会导致内存错误，提示我需要 8TB 内存...

 T90 <- rep(T90, nlyr(hist_temp_stack))

  hist_extreme_temp <- app(hist_temp_stack, fun = function(x) ifelse(x >= T90[], 1, 0))

Answer 1

示例数据

library(terra)
r <- rast(nrow=5, ncol=5, nlyr=100, xmin=0, xmax=1, ymin=0, ymax=1, crs="local")
r <- init(r, runif)

获取第 90 分位数

q <- quantile(r, .9)

q
#class       : SpatRaster 
#dimensions  : 5, 5, 1  (nrow, ncol, nlyr)
#resolution  : 0.2, 0.2  (x, y)
#extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter) 
#source(s)   : memory
#name        :      q0.9 
#min value   : 0.8204033 
#max value   : 0.9477986 

对值 >= 分位数的单元格数量求和

sum(r >= q)

#class       : SpatRaster 
#dimensions  : 5, 5, 1  (nrow, ncol, nlyr)
#resolution  : 0.2, 0.2  (x, y)
#extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter) 
#source(s)   : memory
#name        : sum 
#min value   :  10 
#max value   :  10