R:双线性插值以填补 R 中的空白

Question

我有一个包含要使用插值填充的间隙 (NA) 的网格。我的网格在 x 和 y 维度显示自相关，所以我想尝试双线性插值。我发现的大多数解决方案都集中在“上采样”(用于增加样本数量/网格大小的插值)，但我不想/不需要更改网格大小。我只想使用插值填充 NA。其他潜在的解决方案 do not seem to handle NAs for the input grid of values (the 'z matrix') ，或者是 neighborhood-based solutions 而不是双线性插值，或者只是 have no answer 。

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

我是否遗漏了一些可以使用光栅包进行双线性插值的方法？或者有没有更好的方法来简单地填充 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 提出了一个很好的问题，即为什么不在光栅包中使用带有 focus() 命令的移动窗口平均值。原因是我想尝试双线性插值，我不认为移动窗口平均值总是给出与双线性插值相同的答案。但是，这在我发布的示例中并不清楚(在该示例中移动窗口和双线性 interp 确实给出了相同的答案)，因此我将在下面的新示例中进行演示。请注意，以下示例 ( here is a handy calculator for tests ) 的双线性插值解应为 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)

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

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 值，但在这种情况下可能不是。
`
提取(r，xyFromCell(r，25))
#NA
提取(r，xyFromCell(r，25)，方法=“双线性”)
#[1] 8

那也是什么 resample给

resample(r, r)[25]
# 8

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

这对小的变化非常敏感

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

在这种情况下我真正想要的是车邻居的平均值

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]