R：为什么矩阵3d线性插值不正确？

Question

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

我试图通过采取在我的矩阵中所有已知值，并通过距离值乘以猜测所有这些值（使得更远的一点是，它的影响力就越小）。这是内插结果是什么样子：

正如你所看到的，这种方法不是很有效，它不影响NA s到已知值最接近的，但后来他们迅速收敛到平均值。我认为这是因为它采用了整个范围，这个范围有许多起伏......而不仅仅是距离它最近的点。

显然，矩阵运算并不是我的专业......我需要改变以正确执行线性插值？

下面的代码：

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) 

注：

1. 到akima::interp第一呼叫执行双线性内插。有关用法和详细信息，请参阅帮助页面?akima::interp。

   • 的一个关键点在于，输入x，y和z为已知数据点不必是在x-y网格。在这种情况下，这些是ControlValues的列。
   • 的akima::interp输出是一个列表，其z组分是内插的值超过其x和y坐标由输入xo和yo，分别被定义在网格的矩阵。在这种情况下，这些只是Cont
   • 的行和列索引正如在帮助页面的凸包外点

   z值返回NA说。

   在这种情况下，对应于yo=1:9的输出的第一九列将是NA秒。

2. 到akima::interp（实际上akima::interp.old）第二呼叫执行数据外推来填充在由第一呼叫离开NA秒。有关此用法的详细信息，请参见this SO quation/answer。

上述方法给出以下结果

执行双线性内插是使用interp.surface函数在fields包的另一种方法。提到这种方法是因为实现是一个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[i],y[i])是要内插的坐标。由于ControlValues中的数据点位于网格上，因此这些情况下的这些要求也得到满足。有关用法和详细信息，请参阅帮助页面?interp.surface。

注：

1. sort(unique(ControlValues[,1]))和sort(unique(ControlValues[,2]))简单地给出了x和y坐标已知数据点的网格
2. 该列表中的z组件是简单地用于再成形为已知数据点的z值在已知数据点的网格上的矩阵
3. 要内插的坐标矩阵由make.surface.grid使用为x和y坐标Conf的行和列索引，分别
4. 坐标进行内插的是位于已知点的网格的外部将导致NA
5. interp.surface一个内插的值返回对应于坐标进行内插z值的矢量。这是然后在坐标格内插，它通过ncol(Cont)

最后的尺寸为nrow(Cont) rehaped一个矩阵，很容易验证两种方法产生相同的结果

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