在 R 中用样条插值 rasterstack 时间序列

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) 一起使用，它 returns 一个一般错误：

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

如果我运行它使用f.int.spline(x)它returns一个与内存使用相关的错误：

> 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)

在 R 中用样条插值 rasterstack 时间序列

Interpolating rasterstack time series with spline in R

r

time-series

interpolation

spline

r-raster