我有一个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
您的示例数据(为了清晰起见，更改了层名(
library(terra)
a1 <- rast(ncol = 10, nrow = 10, vals=rep(5,100), names="A")
a2 <- rast(ncol = 10, nrow = 10, vals=rep(10,100), names="A")
b1 <- rast(ncol = 10, nrow = 10, vals=rep(5,100), names="B")
b2 <- rast(ncol = 10, nrow = 10, vals=rep(10,100), names="B")
z <- c(a1, a2, b1, b2)
你的"；terra"；CCD_ 1是不可能的，因为它不处理返回多个数字的函数。但是，您可以使用当前版本的
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
如果使用函数中的部分，则可以执行以下操作：
首先定义哪些层是哪个组：
> 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，也许你的分组不同。。。