R：为什么这个矩阵 3d 线性插值不能正常工作？

Question

我有一个由值和零组成的矩阵，其中零 = NA。这些值散布在矩阵周围，我想做的是对所有 NA 值进行插值。这是数据：

我试图通过获取矩阵中的所有已知值并将该值乘以距离（这样一个点离得越远，它的影响越小）来猜测所有这些值。这是插值结果的样子：

如你所见，这个方法不是很有效，它确实影响了最接近已知值的NA，但随后它们很快收敛到一个平均值。我认为这是因为它采用了整个范围，其中有很多起伏......而不仅仅是最接近它的点。

显然，矩阵运算不是我的专长...我需要更改什么才能正确执行线性插值？

代码如下：

library(dplyr)
library(plotly)

Cont <- structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 1816, 2320, 1406, 2028, 1760, 1932, 1630, 
                    1835, 1873, 1474, 1671, 2073, 1347, 2131, 2038, 1969, 2036, 1602, 
                    1986, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 2311, 1947, 2094, 1947, 2441, 1775, 1461, 1260, 
                    1494, 2022, 1863, 1587, 2082, 1567, 1770, 2065, 1404, 1809, 1972, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                    0, 0, 0, 2314, 1595, 2065, 1870, 2178, 1410, 1994, 1979, 2111, 
                    1531, 1917, 1559, 2109, 1921, 1606, 1469, 1601, 1771, 1771), .Dim = c(19L, 
                                                                                       30L))

  ## First get real control values
  idx <- which(Cont > 0, arr.ind=TRUE)
  V <- Cont[idx]
  ControlValues <- data.frame(idx,V)

  ## Make data.frame of values to fill
  toFill <- which(Cont == 0, arr.ind=TRUE) %>% as.data.frame
  toFill$V <- 0

  ## And now figure out the weighted value of each point
  for (i in 1:nrow(toFill)){
    toFill[i,] -> CurrentPoint

    Xs <- (1/abs(CurrentPoint[,1] - ControlValues[,1])) 
    Xs[is.infinite(Xs)] <- 0
    Xs <- Xs/sum(Xs)/100

    Ys <- (1/abs(CurrentPoint[,2] - ControlValues[,2])) 
    Ys[is.infinite(Ys)] <- 0
    Ys <- Ys/sum(Ys)/100

    ControlValues1 <- data.frame(Xs,Ys)
    toFill[i,3] <- sum(rowMeans(ControlValues1) * ControlValues$V)*100
  }

  ## add back in the controls and reorder
  bind_rows(ControlValues,toFill) -> Both
  Both %>% arrange(row,col) -> Both

  ## and plot the new surface
  NewCont <- matrix(Both$V,max(Both$row),max(Both$col),byrow = T)
  plot_ly(z=NewCont, type="surface",showscale=FALSE)

Answer 1

在 R 中内插和外推数据的一种方法是使用 akima 包。以下执行双线性插值，然后使用数据帧 ControlValues 中的已知数据点作为输入进行外推，以填充 Cont 中的零。

library(akima)
library(plotly)

NewCont <- akima::interp(x=ControlValues[,1], y=ControlValues[,2], z=ControlValues[,3],
                         xo=1:nrow(Cont), yo=1:ncol(Cont), linear=TRUE)$z
NewCont[,1:9] <- akima::interp.old(x=ControlValues[,1], y=ControlValues[,2], 
                                   z=ControlValues[,3], xo=1:nrow(Cont), 
                                   yo=1:9, ncp=2, extrap=TRUE)$z

plot_ly(z=NewCont, type="surface",showscale=FALSE)

备注：

第一次调用 akima::interp 执行双线性插值。有关用法和详细信息，请参阅帮助页面 ?akima::interp。
- 一个关键点是已知数据点的输入 x、y 和 z 不需要在 x-y 网格上。在本例中，这些是 ControlValues 的列。
- akima::interp 的输出是一个列表，其 z 分量是网格上的插值矩阵，其 x 和 y 坐标由输入定义分别为 xo 和 yo。在这种情况下，这些只是 Cont
- 如帮助页面所述
z-values for points outside the convex hull are returned as NA.

在这种情况下，对应于yo=1:9的输出的前九列将是NAs。
第二次调用 akima::interp（实际上是 akima::interp.old）执行数据外推以填充第一次调用留下的 NA。有关此用法的详细信息，请参阅。

上述方法给出了以下结果

执行双线性插值的另一种方法是使用 fields 包中的 interp.surface 函数。提到这种方法是因为实现是一个R脚本，可以通过在R命令行输入函数名interp.surface来列出。

library(fields)

loc <- make.surface.grid(list(x=1:nrow(Cont), y=1:ncol(Cont)))
NewCont2 <- matrix(interp.surface(list(x=sort(unique(ControlValues[,1])),
                                       y=sort(unique(ControlValues[,2])),
                                       z=matrix(ControlValues[,3],
                                                nrow=length(unique(ControlValues[,1])),
                                                ncol=length(unique(ControlValues[,2])))),
                                  loc), nrow=nrow(Cont), ncol=ncol(Cont))
NewCont2[,1:9] <- akima::interp.old(x=ControlValues[,1], y=ControlValues[,2], 
                                    z=ControlValues[,3], xo=1:nrow(Cont), 
                                    yo=1:9, ncp=2, extrap=TRUE)$z

此处的要求与akima::interp相反。具体来说，已知数据点必须位于 x-y 网格上。但是，要插值的坐标不需要在网格上，而是一个包含 x 和 y 坐标的相应列向量的矩阵，其中每个元组 (x[i],y[i]) 是一个 x-y 坐标插值。由于 ControlValues 中的数据点在网格上，因此本例也满足这些要求。有关用法和详细信息，请参阅帮助页面 ?interp.surface。

备注：

sort(unique(ControlValues[,1])) 和 sort(unique(ControlValues[,2])) 只是给出已知数据点网格的 x 和 y 坐标
列表中的 z 分量只是已知数据点的 z 值，在已知数据点的网格上重塑为矩阵
要插值的坐标矩阵由make.surface.grid生成，分别使用x和y坐标Conf的行和列索引
位于已知点网格之外的插值坐标将导致插值 NA
interp.surface returns 对应于要插值的坐标的 z 值的矢量。然后将其转换为要插值的坐标网格上的矩阵，其尺寸为 nrow(Cont) by ncol(Cont)

最后，很容易验证两种方法给出相同的结果

print(max(abs(NewCont - NewCont2)))
##[1] 4.547474e-13

R：为什么这个矩阵 3d 线性插值不能正常工作？

R: Why isn't this matrix 3d linear interpolation working correctly?

interpolation

r

matrix

linear-interpolation