R：双线性插值以填补 R 中的空白
R: Bilinear interpolation to fill gaps in R
我有一个网格，其中包含我想使用插值法填充的间隙 (NA)。我的网格在 x 和 y 维度上显示自相关，所以我想尝试双线性插值。我发现的大多数解决方案都集中在 'upsampling'（为了增加 samples/size 网格的数量而进行插值），但我不会 want/need 来更改网格大小。我只想使用插值法填充 NA。其他可能的解决方案 do not seem to handle NAs for the input grid of values (the 'z matrix'), or are rather than bilinear interpoloation, or simply have no answer。
我发现使用栅格包，我可以输入包含 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 提出了一个很好的问题，为什么不在栅格包中使用移动 window 平均值和 focal() 命令。原因是我想尝试双线性插值，我不认为移动 window 平均值总是给出与双线性插值相同的答案。但是，这在我发布的示例中并不清楚（在该示例中移动 window 和双线性插值 do 给出相同的答案），所以我将在一个新示例中进行演示以下。请注意，对于下例 (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)
我是不是漏掉了什么？也许 focal() 的 weights 参数有一些聪明的东西可以产生双线性插值解决方案？
让我们使用等距离的单元格来避免由于 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)
我看到你解雇了 "neighborhood-based solutions rather than bilinear interpoloation"。但问题是为什么。在这种情况下，您可能需要基于邻域的解决方案。
更新。再一次，如果单元格不是近似正方形，则双线性将更可取。
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 值（例如 land/ocean），但在这种情况下可能不行。
`
提取物（r，xyFromCell（r，25））
#NA
提取物（r，xyFromCell（r，25），方法='bilinear'）
#[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]
我有一个网格，其中包含我想使用插值法填充的间隙 (NA)。我的网格在 x 和 y 维度上显示自相关，所以我想尝试双线性插值。我发现的大多数解决方案都集中在 'upsampling'（为了增加 samples/size 网格的数量而进行插值），但我不会 want/need 来更改网格大小。我只想使用插值法填充 NA。其他可能的解决方案 do not seem to handle NAs for the input grid of values (the 'z matrix'), or are
我发现使用栅格包，我可以输入包含 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 提出了一个很好的问题，为什么不在栅格包中使用移动 window 平均值和 focal() 命令。原因是我想尝试双线性插值，我不认为移动 window 平均值总是给出与双线性插值相同的答案。但是，这在我发布的示例中并不清楚（在该示例中移动 window 和双线性插值 do 给出相同的答案），所以我将在一个新示例中进行演示以下。请注意，对于下例 (here is a handy calculator for tests).
，双线性插值解应为 8library(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)
我是不是漏掉了什么？也许 focal() 的 weights 参数有一些聪明的东西可以产生双线性插值解决方案？
让我们使用等距离的单元格来避免由于 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)
我看到你解雇了 "neighborhood-based solutions rather than bilinear interpoloation"。但问题是为什么。在这种情况下，您可能需要基于邻域的解决方案。
更新。再一次，如果单元格不是近似正方形，则双线性将更可取。
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 值（例如 land/ocean），但在这种情况下可能不行。
`
提取物（r，xyFromCell（r，25））
#NA
提取物（r，xyFromCell（r，25），方法='bilinear'）
#[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]