在R中用样条插值栅格堆栈时间序列

Question

我正在尝试使用样条线填补 ndvi 图像时间序列中的空白。

我已经用我拥有的 ndvi 图像和一些只有NA 的图层创建了一个光栅堆栈，用于我没有的时间步长。 光栅堆栈在此处作为 Geotiff 提供：raster stack download

我编写了以下函数来插入缺失值：

f.int.spline <- function(x) {    # x is a raster stack or brick
              v=as.vector(x)     # transform x in a vector for easier manipulation
              z=which(v %in% NA) # find which pixel values are to be interpolated
              # interpolation with spline
              interp <- spline(x=c(1:NROW(v)), y = v,  
                 xout = z,       # missing values to be interpolated
                 method = "natural") 

x[z] <-interp$y # including the missing values in the original vector

}

如果我将它与一个像素一起使用（例如x[ 50, 200 ]），该函数可以工作，但如果我使用calc(x, f.int.spline) 运行它，它会返回一个一般错误：

> error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

如果我使用f.int.spline(x) 运行它，它会返回与内存使用相关的错误：

> Error in as.vector(t((matrix(rep(i, rec2), nrow = rec2, byrow = TRUE)) +  : 
  error in evaluating the argument 'x' in selecting a method for function 'as.vector': Error in t((matrix(rep(i, rec2), nrow = rec2, byrow = TRUE)) + add) : 
  error in evaluating the argument 'x' in selecting a method for function 't': Error: cannot allocate vector of size 4.9 Gb

1) 您是否发现任何缺陷或有任何解决方法来使其发挥作用？

2) 我无法准确理解 calc() 函数如何用于光栅堆栈：它是否采用所有层中每个像素的值？

3) 正如 Jeffrey Evans 所建议的那样，我也在寻找更适合这项工作的其他插值函数。有什么想法吗？

Answer 1

首先创建一个适用于向量的函数，包括一些极端情况（可能与您相关也可能不相关）

 f <- function(x) {  
    z <- which(is.na(x))
    nz <- length(z)
    nx <- length(x)
    if (nz > 0 & nz < nx) { 
        x[z] <- spline(x=1:nx, y=x, xout=z, method="natural")$y
    }
    x
}

测试功能

f(c(1,2,3,NA,5,NA,7))
##[1] 1 2 3 4 5 6 7
f(c(NA,NA,5,NA,NA))
##[1] 5 5 5 5 5
f(rep(NA, 8))
##[1] NA NA NA NA NA NA NA NA
f(rep(1, 8))
##[1] 1 1 1 1 1 1 1 1

然后在 RasterStack 或 RasterBrick 上使用 calc

示例数据

r <- raster(ncols=5, nrows=5)
r1 <- setValues(r, runif(ncell(r)))
r2 <- setValues(r, runif(ncell(r)))
r3 <- setValues(r, runif(ncell(r)))
r4 <- setValues(r, runif(ncell(r)))
r5 <- setValues(r, NA)
r6 <- setValues(r, runif(ncell(r)))
r1[6:10] <- NA
r2[5:15] <- NA
r3[8:25] <- NA
s <- stack(r1,r2,r3,r4,r5,r6)
s[1:5] <- NA

使用函数

x <- calc(s, f)

或者，您可以使用approxNA

x <- approxNA(s)