R: Warum funktioniert diese Matrix 3d lineare Interpolation nicht korrekt?

Question

Ich habe eine Matrix von Werten und Nullen, wobei Null = NA. Die Werte sind um die Matrix verteilt und ich möchte die Werte aller NA Werte interpolieren. Dies sind die Daten:

Ich versuche, alle diese Werte zu erraten, indem er alle bekannten Werte in meiner Matrix nehmen und den Wert durch den Abstand multipliziert wird (so dass die weiter entfernt ist ein Punkt , je weniger Einfluss es hat). Dies ist, was das interpolierte Ergebnis sieht so aus:

Wie Sie sehen können, ist diese Methode nicht sehr wirksam ist, es hat die NA s am nächsten zu den bekannten Werten beeinflussen, aber dann konvergieren sie schnell auf einen Durchschnittswert . Ich denke, das liegt an der Tatsache, dass es das GESAMTE BEREICH nimmt, das viele Höhen und Tiefen hat ... und nicht nur die nächsten Punkte.

Offensichtlich sind Matrixoperationen nicht meine Spezialität ... Was muss ich ändern, um die lineare Interpolation korrekt durchzuführen?

Hier ist der Code:

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

Ein Ansatz zur Interpolation und extrapolieren Daten in R ist das akima Paket zu verwenden. Im Folgenden wird eine bilineare Interpolation gefolgt von einer Extrapolation durchgeführt, wobei die bekannten Datenpunkte im Datenrahmen ControlValues als Eingabe verwendet werden, um die Nullen in Cont zu füllen.

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) 

Anmerkungen:

1. Der erste Aufruf von akima::interp führt die bilineare Interpolation. Weitere Informationen finden Sie auf der Hilfeseite ?akima::interp.

   • Ein wesentlicher Punkt ist, dass die Eingänge x, y und z für die bekannten Datenpunkte nicht auf einem x-y Gitter zu sein brauchen. In diesem Fall sind dies die Spalten ControlValues.
   • Der Ausgang des akima::interp ist eine Liste, deren z Komponente ist eine Matrix von interpolierten Werten über das Gitter, dessen x und y Koordinaten durch die Eingänge xo und yo jeweils definiert sind. In diesem Fall sind dies nur die Zeilen- und Spaltenindizes Cont
   • Wie in der Hilfeseite erklärt

   z-Werte für die Punkte außerhalb der konvexen Hülle als NA zurückgegeben.

   In diesem Fall sind die ersten neun Spalten des Ausgangs auf yo=1:9 entspricht, wird NA s sein.

2. Der zweite Aufruf von akima::interp (eigentlich akima::interp.old) führt die Daten Extrapolation in dem NA s durch den ersten Anruf zu füllen. Einzelheiten zu dieser Verwendung finden Sie unter this SO quation/answer.

Der obige Ansatz ergibt folgendes Ergebnis

Einen anderen Ansatz bilineare Interpolation durchführt, die in der Funktion interp.surfacefields Paket zu verwenden ist. Dieser Ansatz wird erwähnt, da die Implementierung ein R-Skript ist, das durch Eingabe des Funktionsnamens interp.surface in der R-Befehlszeile aufgelistet werden kann.

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 

Hier sind die Anforderungen, die entgegengesetzt zu denen für akima::interp. Konkret müssen die bekannten Datenpunkte auf einem x-y Raster liegen. Die zu interpolierenden Koordinaten müssen jedoch nicht auf einem Gitter liegen und sind stattdessen eine Matrix, die entsprechende Spaltenvektoren von x und y Koordinaten enthält, wobei jedes Tupel (x[i],y[i]) eine x-y Koordinate zum Interpolieren ist. Da die Datenpunkte in ControlValues auf einem Raster liegen, sind diese Anforderungen auch für diesen Fall erfüllt. Informationen zur Verwendung und Details finden Sie auf der Hilfeseite ?interp.surface.

Anmerkungen:

1. sort(unique(ControlValues[,1])) und sort(unique(ControlValues[,2])) einfach geben die x und y Koordinaten für das Raster der bekannten Datenpunkte
2. Die z Komponente in der Liste ist einfach die z Werte für die bekannten Datenpunkte neu geformt, wie eine Matrix über dem Gitter von bekannten Datenpunkten
3. Die Matrix der zu interpolierenden Koordinaten wird erzeugt durch make.surface.grid unter Verwendung als x und y koordiniert die Zeilen- und Spaltenindizes Conf jeweils
4. Koordinatenmeßgerät zu interpolieren, die außerhalb des Gitters von bekannten Punkten liegt in einem interpolierten Wert von NA
5. interp.surface liefert einen Vektor von z Werte entsprechend den Koordinaten führt zu interpolieren, . Diese wird dann an eine Matrix über das Gitter von Koordinaten rehaped zu interpolieren, die Dimensionen nrow(Cont) von ncol(Cont)

Schließlich hat, ist es leicht, zu überprüfen, ob die beiden Ansätze das gleiche Ergebnis

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