这是一种方法
library(raster)
r <- raster(nrows=10, ncols=10);
values(r) <- NA
x <- sapply(1:30, function(...) r)
x[[1]] <- setValues(r, runif(ncell(r)))
x[[16]] <- setValues(r, runif(ncell(r))) + 10
x[[30]] <- setValues(r, runif(ncell(r))) + 20
s <- stack(x)
z <- approxNA(s)
plot(z)
plot(1:30, z[1])
这是另一种方法
library(raster)
r <- raster(nrows=10, ncols=10);
x1 <- setValues(r, runif(ncell(r)))
x16 <- setValues(r, runif(ncell(r))) + 10
x30 <- setValues(r, runif(ncell(r))) + 20
s <- stack(x1, x16, x30)
x <- calc(s, fun=function(y) approx(c(1,16,30), y, 1:30)$y)
但如果三层中存在 NA 值,则此操作将会失败。您需要调整功能fun
来处理这个问题(这是一个example).