r - R：双线性插值以填补 R 中的空白

Question

我有一个网格，其中包含要使用插值填充的间隙 (NA)。我的网格显示 x 和 y 维度的自相关，所以我想尝试双线性插值。我发现的大多数解决方案都集中在“上采样”（用于增加样本数量/网格大小的插值），但我不想/不需要更改网格大小。我只想使用插值填充 NA。其他潜在的解决方案似乎没有处理输入值网格（“z 矩阵”）的 NA，或者是基于邻域的解决方案而不是双线性插值，或者根本没有答案。

我发现使用 raster 包，我可以输入一个包含 NA 的网格（作为光栅），并使用“resample”命令输出相同大小的网格。然而，结果看起来像最近邻插值而不是双线性插值。

我是否遗漏了一些东西，以至于有一种方法可以使用 raster 包进行双线性插值？或者有没有更好的方法来做双线性插值来填充 NA？

library(raster)

# raster containing gap
r <- raster(nrow=10, ncol=10)
r[] <- 1:ncell(r)
r[25] <- NA

# The s raster is the same size as the r raster
s <- raster(nrow=10, ncol=10)
s <- resample(r, s, method='bilinear')
plot(r)
plot(s)
s[25]
s[35]
# s[25] appears to have been filled with neighbor s[35]

更新

Akima 包似乎是上述光栅方法的一个有前途的替代方案，但如果输入值网格（Z 矩阵）中有 NA，我会遇到麻烦。这是一个与上面示例平行的示例来演示。（再次，我插入到与原始大小相同的网格）。

library(akima)

# Use bilinear interpolation (no NAs in input)
rmat<-matrix(seq(1,100,1), nrow = 10, ncol = 10, byrow = T)
x <- seq(1,10,1)
y <- seq(1,10,1)
smat <- bilinear.grid(x, y, rmat, nx = 10, ny = 10) # works
plot(raster(rmat), main = "original")
plot(raster(smat$z), main = "interpolated")

# Try using bilinear interpolation but with an NA
rmat<-matrix(seq(1,100,1), nrow = 10, ncol = 10, byrow = T)
rmat[3,5] <- NA
x <- seq(1,10,1)
y <- seq(1,10,1)
smat <- bilinear.grid(x, y, rmat, nx = 10, ny = 10) # Error about NAs

更新2

@Robert Hijmans 提出了一个很好的问题，即为什么不在 raster 包中的 focus() 命令中使用移动窗口平均值。原因是我想尝试双线性插值，我不认为移动窗口平均值总是给出与双线性插值相同的答案。但是，这在我发布的示例中并不清楚（在该示例中移动窗口和双线性 interp确实给出了相同的答案），因此我将在下面的新示例中进行演示。请注意，对于下面的示例，双线性插值解应该是 8（这里是一个方便的测试计算器）。

library(raster)
r <- raster(nrow=10, ncol=10)

# Different grid values than earlier examples
values(r) <- c(rep(1:5, 4), rep(4:8, 4), rep(1:5, 4), rep(4:8, 4), rep(1:5, 4))
r[25] <- NA
plot(r)

# See what the mean of the moving window produces
f <- focal(r, w=matrix(1,nrow=3, ncol=3), fun=mean, NAonly=TRUE, na.rm=TRUE) 
f[25] # Moving window gives 5 but bilinear interp gives 8

# Note that this seems to be how the moving window works with equal weights
window_test <- c(r[14:16], r[24:26], r[34:36])
 mean(window_test, na.rm = T)

我在这里错过了什么吗？也许焦点（）的权重参数有一些聪明的东西可以产生双线性插值解决方案？

标签： rinterpolationgeospatial

Answer 1

让我们使用等距离的单元格来避免由于单元格大小与 lon/lat 数据的差异而产生的差异

library(raster)
r <- raster(nrow=10, ncol=10, crs='+proj=utm +zone=1 +datum=WGS84', xmn=0, xmx=1, ymn=0, ymx=1)

对于此示例，您可以使用focal

values(r) <- 1:ncell(r)
r[25] <- NA
f <- focal(r, w=matrix(1,nrow=3, ncol=3), fun=mean, NAonly=TRUE, na.rm=TRUE) 

我看到您驳回了“基于邻域的解决方案而不是双线性插值”。但问题是为什么。在这种情况下，您可能需要基于邻域的解决方案。

更新。再说一次，如果单元格不是近似正方形，双线性将是更可取的。

values(r) <- c(rep(1:5, 4), rep(4:8, 4), rep(1:5, 4), rep(4:8, 4), rep(1:5, 4))
r[25] <- NA

双线性插值的问题通常使用 4 个连续的单元格，但在这种情况下，如果您想要单元格中心的值，则适当的单元格将是单元格本身的值，因为到该单元格的距离为零，并且因此，这就是插值结束的地方。例如，对于单元格 23

extract(r, xyFromCell(r, 23))
#6
extract(r, xyFromCell(r, 23), method='bilinear')
#[1] 6

在这种情况下，焦点单元格为 NA，因此您得到焦点单元格和另外 3 个单元格的平均值。问题是哪三个？这是任意的，但要使其工作，NA 单元必须获得一个值。该raster算法将 NA 单元格下方的值分配给该单元格（此处也是 8）。我认为这可以很好地处理边缘（例如陆地/海洋）的 NA 值，但在这种情况下可能不是。` extract(r, xyFromCell(r, 25)) #NA extract(r, xyFromCell(r, 25), method='bilinear') #[1] 8

这也是resample给

resample(r, r)[25]
# 8

这也是在线计算器的建议吗？

这对微小的变化非常敏感

extract(r, xyFromCell(r, 25)+0.0001, method='bilinear')
#[1] 4.998997

在这种情况下，我真正想要的是 rook-neighbors 的平均值

mean(r[adjacent(r, 25, pairs=FALSE)])
[1] 6

或者，更一般地，局部反距离加权平均值。您可以通过设置具有焦点的权重矩阵来计算

# compute weights matrix
a <- sort(adjacent(r, 25, 8, pairs=F, include=TRUE))
axy <- xyFromCell(r, a)
d <- pointDistance(axy, xyFromCell(r, 25), lonlat=F)
w <- matrix(d, 3, 3)
w[2,2] <- 0
w <- w / sum(w)

# A simpler approach could be: 
# w <- matrix(c(0,.25,0,.25,0,.25,0,.25,0), 3, 3)


foc <- focal(r, w, na.rm=TRUE, NAonly=TRUE)
foc[25]

在这个例子中，这很好；但如果焦点区域中有多个 NA 值（因为权重之和不再为 1），则将不正确。我们可以通过计算权重之和来纠正这一点

x <- as.integer(r/r)
sum_weights <- focal(x, w, na.rm=TRUE, NAonly=TRUE)

fw <- foc/sum_weights
done <- cover(r, fw)
done[25]