使用{terra}？计算一组栅格堆栈图层的统计数据

bt1cpqcv 于 2023-01-18 发布在其他

关注(0)|答案(2)|浏览(172)

我有一个4层的栅格堆栈。其中两层来自模型1，两层来自模型2。我需要计算每个模型的中值、第5百分位数和第95百分位数。有没有办法一步完成？即，不写出两个中间的栅格堆栈，然后再将它们连接在一起。我的尝试如下，但它没有按组执行函数。

library("terra")   
# Create some toy data
a <- rast(ncol = 10, nrow = 10, vals=rep(5,100), names=1)
b <- rast(ncol = 10, nrow = 10, vals=rep(10,100), names=1)
c <- rast(ncol = 10, nrow = 10, vals=rep(5,100), names=2)
d <- rast(ncol = 10, nrow = 10, vals=rep(10,100), names=2)
z <- c(a, b, c, d)

# Try to write a function to do the work
app(z,
    function(x) {
      c(median(x), quantile(x, c(0.05, 0.95)))
      },
     filename = "grouped_stats.tif)

我想要的结果是一个6层的光栅堆栈。

class       : SpatRaster
dimensions  : 10, 10, 6  (nrow, ncol, nlyr)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
sources     : memory  (3 layers)
              memory  (3 layers)
names       : median_1, q5_1, q95_1, median_2, pc5_2, pc95_2
min values  :      7.5,  5.0,  10.0,      7.5,   5.0,   10.0
max values  :      7.5,  5.0,  10.0,      7.5,   5.0,   10.0

有什么主意吗？谢谢。

作用力1

受@spacedman的启发，我写了这个函数，但它并没有让我达到目的。把它放在这里是为了给其他人提供可能的灵感。

grouped_stats <- function(x) {
  layers_names <- unique(names(x))
  cell_output <- NA
  for (each_layer in layers_names) {
     cell_output <- rbind(cell_output,
                    c(median(x[[each_layer]], na.rm = TRUE),
                      quantile(x[[each_layer]], 0.05, 0.95)))
     names(cell_output) <- glue("{each_layer}_{c('median','pc5','pc95')}")
  }
  cell_output
}

g <- app(z, fun = grouped_stats)

作用力2

我想越来越近了，但还没到。

my_stats_function <- function(x) {c(median(x), quantile(0.05, 0.95))}

app(z, 
    function(x){
      unlist(tapply(x, layer_names, my_stats_function))
      })

class       : SpatRaster 
dimensions  : 10, 10, 4  (nrow, ncol, nlyr)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source      : memory
names       :   11, 1.95%,   21, 2.95%
min values  : 7.50,  0.05, 7.50,  0.05
max values  : 7.50,  0.05, 7.50,  0.05

努力3

我想我快到了。：-）

my_stats_function <- function(x) {c(median(x), quantile(x, c(0.05, 0.95)))}

app(z, 
    function(x){
      unlist(tapply(x, layer_names, my_stats_function))
      })

class       : SpatRaster
dimensions  : 10, 10, 6  (nrow, ncol, nlyr)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source      : memory
names       : 11, 1.5%, 1.95%, 21, 2.5%, 2.95%
min values  :  5,    5,     5,  5,    5,     5
max values  :  5,    5,     5,  5,    5,     5

来源：https://stackoverflow.com/questions/74012639/calculate-stats-on-a-group-of-raster-stack-layers-using-terra

2条答案

按热度按时间

pgvzfuti1#

示例数据（为清楚起见，更改了图层名称）

library(terra)
a <- rast(ncol = 10, nrow = 10, vals=rep(5,100), names="A")
b <- rast(ncol = 10, nrow = 10, vals=rep(10,100), names="A")
c <- rast(ncol = 10, nrow = 10, vals=rep(5,100), names="B")
d <- rast(ncol = 10, nrow = 10, vals=rep(10,100), names="B")
z <- c(a, b, c, d)

您的“terra”版本不可能tapp，因为它不能处理返回多个数字的函数。

f <- function(x) quantile(x, c(0.5, 0.05, 0.95))
x <- tapp(z, names(z), f)
x
#class       : SpatRaster 
#dimensions  : 10, 10, 6  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#names       : A_5%, A_50%, A_95%, B_5%, B_50%, B_95%
#min values  : 5.25,   7.5,  9.75, 5.25,   7.5,  9.75 
#max values  : 5.25,   7.5,  9.75, 5.25,   7.5,  9.75

我认为这是最干净的解决方案，但在本例中，因为有quantile<SpatRaster>方法，所以它 * 可能 * 执行起来更快

probs <- c(0.05, 0.5, 0.95) 
ids <- names(z)
uids <- unique(ids)
x <- lapply(uids, function(i) quantile(z[[ids == i]], probs))
rast(x)
#class       : SpatRaster 
#dimensions  : 10, 10, 6  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#names       : q0.05, q0.5, q0.95, q0.05, q0.5, q0.95 
#min values  :  5.25,  7.5,  9.75,  5.25,  7.5,  9.75 
#max values  :  5.25,  7.5,  9.75,  5.25,  7.5,  9.75

但我只会在光栅非常大的情况下才看。
我认为在所有情况下，使用

function(x) quantile(x, c(0.5, 0.05, 0.95))

代替

function(x) c(median(x), quantile(x, c(0.05, 0.95)))

下面是使用循环和app的效率较低的解决方案

f <- function(x) quantile(x, c(0.05, 0.5, 0.95)) 
ids <- names(z)
uids <- unique(ids)
x <- lapply(uids, function(i) app(z[[ids == i]], f))
rast(x)
#class       : SpatRaster 
#dimensions  : 10, 10, 6  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#sources     : memory  (3 layers) 
#              memory  (3 layers) 
#names       :   5%, 50%,  95%,   5%, 50%,  95% 
#min values  : 5.25, 7.5, 9.75, 5.25, 7.5, 9.75 
#max values  : 5.25, 7.5, 9.75, 5.25, 7.5, 9.75

赞(0）回复(0）举报 2023-01-18

uinbv5nw2#

如果使用函数中的零件，则可以执行以下操作：
首先定义哪些层属于哪个组：

> p1 = c(1,3)
> p2 = c(2,4)

然后将x设为子集并进行处理：

> app(z, function(x){c(median(x[p1]), quantile(x[p1], c(0.05, 0.95)), median(x[p2], ), quantile(x[p2], c(0.05, 0.95)))})
class       : SpatRaster                  
dimensions  : 10, 10, 6  (nrow, ncol, nlyr)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : memory 
names       : lyr.1, lyr.2, lyr.3, lyr.4, lyr.5, lyr.6 
min values  :     5,     5,     5,    10,    10,    10 
max values  :     5,     5,     5,    10,    10,    10

我不知道你是怎么得到你的样本中值为7.5的，所以也许你用了mean，也许你的分组是不同的...

赞(0）回复(0）举报 2023-01-18

我来回答

使用{terra}？计算一组栅格堆栈图层的统计数据

作用力1

作用力2

努力3

2条答案

相关问题

热门标签

最新问答