我有一些大型栅格数据,代表时间序列(每个图层是一个时间点)。我希望通过在层之间进行插值来增加时间分辨率(在下面的示例中,使用
接下来,我们将
zoo::na.spline
进行插值,尽管如果有助于解决问题,我也可以使用其他插值方法)。我不能直接在栅格数据上使用na.spline
- 它不是为此设计的,并且只会给出Error in as.array.default(X) : 'dimnames' applied to non-array
。
目前,我正在使用一系列循环来完成此操作(第1个循环用于生成具有更多时间层的栅格数据,第2个循环用于将原始栅格数据插入到其中的正确层中,第3个循环用于对每个单元格使用na.spline)。然而,在大型栅格数据上这种方法非常缓慢,似乎是一种相当低效的方法。
一个最小的可重现示例。
首先让我们制作一个表示低时间分辨率数据的初始栅格数据集。请注意,左上角的单元格始终为NA,因为任何解决方案都必须对仅包含NA的单元格进行鲁棒性检查:
library(raster); library(rasterVis); library(zoo)
r1 = raster(matrix(c(NA, rep(1,8)), 3, 3))
r2 = raster(matrix(c(NA, rep(2,8)), 3, 3))
r3 = raster(matrix(c(NA, rep(3,8)), 3, 3))
b1 = brick(list(r1,r2,r3))
levelplot(b1)
现在,让我们创建一个空的栅格砖块,以便插入值:b2 = brick(lapply(1:9, function(x) raster(matrix(NA, 3, 3))))
接下来,我们将
b1
的图层插入到b2
的相应图层中。请注意,我甚至在这里使用了一个循环,因为将选择的图层子分配到光栅砖中不起作用。old.layers = c(1,4,7)
for (i in 1:nlayers(b1)) {
b2[[old.layers[i]]] = b1[[i]]
}
levelplot(b2, layout=c(9,1))
最后,我们循环遍历每个单元格进行时间序列插值,并将结果插入回b2
中:for (cell in 1:ncell(b2)) {
if (all(is.na(as.vector(b2[cell])))) {
# na.spline wil fail if all are NA
new.values = as.numeric(rep(NA, dim(b2)[3]))
} else {
new.values = na.spline(as.vector(b2[cell]), na.rm = F)
}
b2[cell][] = new.values
}
levelplot(b2, layout=c(9,1))
这个方法可以实现,但似乎非常低效和不优雅。有没有更快(最好也更优雅)的方法来做到这一点?
b2 <- calc(b1, fun)
。我在?calc中看到,“如果fun返回多个数字,则返回RasterBrick”,但出于某种原因,我没有得到这种行为。 - dwwcalc
可能返回一个砖块而不是单层,所以这是一个非常有用的学习练习。 - dww