Interpolierende Rasterstack-Zeitreihe mit Spline in R

Question

Ich versuche, die Lücken in einer Zeitreihe von ndvi Bildern mit Spline zu füllen.

Ich habe einen Raster-Stack mit den ndvi-Bildern erstellt, die ich habe und einige Schichten mit nur NA für die Zeitschritte, die ich nicht habe. Die Rasterstapel ist als GeoTiff hier: raster stack download

ich die folgende Funktion geschrieben haben, die fehlenden Werte zu interpolieren:

f.int.spline <- function(x) { # x is a raster stack or brick 
       v=as.vector(x)  # transform x in a vector for easier manipulation 
       z=which(v %in% NA) # find which pixel values are to be interpolated 
       # interpolation with spline 
       interp <- spline(x=c(1:NROW(v)), y = v, 
       xout = z,  # missing values to be interpolated 
       method = "natural") 

x[z] <-interp$y # including the missing values in the original vector 

} 

Die Funktion funktioniert, wenn ich es mit einem Pixel (zB x[ 50, 200 ]), aber wenn ich es mit calc(x, f.int.spline) laufen gibt es einen allgemeinen Fehler:

> error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
    cannot use this function 

wenn ich es laufen mit f.int.spline(x) es einen Fehler zu Speicherauslastung bezogen zurückgibt:

> Error in as.vector(t((matrix(rep(i, rec2), nrow = rec2, byrow = TRUE)) + : 
    error in evaluating the argument 'x' in selecting a method for function 'as.vector': Error in t((matrix(rep(i, rec2), nrow = rec2, byrow = TRUE)) + add) : 
    error in evaluating the argument 'x' in selecting a method for function 't': Error: cannot allocate vector of size 4.9 Gb 

1) Sehen Sie Fehler oder haben Sie Workarounds, um das Problem zu beheben?

2) Ich kann nicht genau verstehen, wie die calc() Funktion für Rasterstapel funktioniert: Nimmt es die Werte jedes Pixels in allen Schichten?

3) Wie von Jeffrey Evans vorgeschlagen Ich suche auch nach anderen Interpolationsfunktionen, die besser für den Job geeignet sind. Irgendeine Idee ?

Answer 1

zuerst eine Funktion erstellen, die auf einem Vektor arbeitet, auch auf einigen Fällen Ecke (das kann oder auch nicht für Sie relevant sein)

f <- function(x) { 
    z <- which(is.na(x)) 
    nz <- length(z) 
    nx <- length(x) 
    if (nz > 0 & nz < nx) { 
     x[z] <- spline(x=1:nx, y=x, xout=z, method="natural")$y 
    } 
    x 
} 

Testen Sie die Funktion

f(c(1,2,3,NA,5,NA,7)) 
##[1] 1 2 3 4 5 6 7 
f(c(NA,NA,5,NA,NA)) 
##[1] 5 5 5 5 5 
f(rep(NA, 8)) 
##[1] NA NA NA NA NA NA NA NA 
f(rep(1, 8)) 
##[1] 1 1 1 1 1 1 1 1 

Dann nutzen calc auf einem RasterStack oder RasterBrick

Beispieldaten

r <- raster(ncols=5, nrows=5) 
r1 <- setValues(r, runif(ncell(r))) 
r2 <- setValues(r, runif(ncell(r))) 
r3 <- setValues(r, runif(ncell(r))) 
r4 <- setValues(r, runif(ncell(r))) 
r5 <- setValues(r, NA) 
r6 <- setValues(r, runif(ncell(r))) 
r1[6:10] <- NA 
r2[5:15] <- NA 
r3[8:25] <- NA 
s <- stack(r1,r2,r3,r4,r5,r6) 
s[1:5] <- NA 

Verwenden Sie die Funktion

x <- calc(s, f) 

Alternativ können Sie approxNA

x <- approxNA(s)