2016-12-09 6 views
0

Ich versuche, die krige Funktion im gstat Paket von R zu verwenden, um einige räumliche Ozeantiefendaten in R zu interpolieren. Ich finde für mehr als ungefähr ~ 1000 Punkte Funktion beginnt unzumutbare Mengen an Zeit zu beenden (dh Stunden bis Tage zu nie beendet). Ist das normal oder mache ich etwas falsch? Ich bin besonders besorgt, weil mein letztendliches Ziel die räumlich-zeitliche Kriggerung eines sehr großen Datensatzes (> 30.000 Datenpunkte) ist, und ich bin besorgt, dass es angesichts dieser Laufzeiten einfach nicht machbar sein wird.Wie langsam ist zu langsam beim Kriging mit Gstat in R

Ich führe gstat-1.1-3 und R-3.3.2. Ich betreibe Sie den Code unten ist:

library(sp); library(raster); library(gstat) 
v.utm # SpatialPointsDataFrame with >30,000 points 

# Remove points with identical positons 
zd = zerodist(v.utm) 
nzd = v.utm[-zd[,1],] # Layer with no identical positions 

# Make a raster layer covering point layer 
resolution=1e4 
e = extent(as.matrix([email protected]))+resolution 
r = raster(e,resolution=resolution) 
proj4string(r) = proj4string(v.utm) 

# r is a 181x157 raster 

# Fit variogram 
fv = fit.variogram(variogram(AVGDEPTH~1, nzd),model=vgm(6000,"Exp",1,5e5,1)) 

# Krige on random sample of 500 points - works fine 
size=500 
ss=nzd[sample.int(nrow(nzd),size),] 
depth.krig = krige(AVGDEPTH~1,ss,as(r,"SpatialPixelsDataFrame"), 
       model=depth.fit) 

# Krige on random sample of 5000 points - never seems to end 
size=5000 
ss=nzd[sample.int(nrow(nzd),size),] 
depth.krig = krige(AVGDEPTH~1,ss,as(r,"SpatialPixelsDataFrame"), 
       model=depth.fit) 

Antwort

1

Die Komplexität der Choleski Zersetzung (oder ähnlich) ist O (n^3), was bedeutet, dass, wenn Sie die Anzahl der Punkte mit 10 multiplizieren, die Zeit wird es dauern steigt mit dem Faktor 1000. Es gibt zwei Möglichkeiten, dieses Problem, zumindest so weit wie gstat betrifft:

  1. installieren eine optimierte Version von BLAS (zB OpenBLAS oder MKL) - das löst nicht die O (n^3) Problem, kann aber maximal einen Faktor n mit der Anzahl der verfügbaren Kerne beschleunigen.
  2. Vermeiden Sie die vollständige Zerlegung der Kovarianzmatrix durch die Auswahl lokale Nachbarschaften (Argumente maxdist und/oder nmax)
+0

Danke für die hilfreiche Antwort. Ich habe nmax auf einen vernünftigen Betrag gesetzt und das hat die Dinge sehr schnell gemacht. Ich frage mich, ob es etwas Ähnliches zu beschleunigen VariogramST gibt? Ich habe versucht, die Twindow und Cutoff Argumente zu setzen, aber Variogramm Berechnung ist immer noch sehr langsam. Der genaue Code, den ich verwende, ist: vst = variogramST (BOTTEMP ~ AVGDEPTH + lat, Daten = v, cutoff = 7e5, twindow = 5 * 365,25, tlags = (0:15) * 365,25, tunit = "Tage ") – user3004015

+0

Variogramm-Berechnung ist eine O (n^2) -Operation. –

+0

Ich denke, meine Frage ist, ob die Argumente "cutoff" und "twindow" die Anzahl der Berechnungen im selben Sinne reduzieren, wie die Argumente nmax und maxdist für 'krige'? Außerdem glaube ich nicht, dass ich das 'twindow' Argument korrekt in dem von mir geposteten Code verwendet habe. Können Sie bestätigen, dass die Einheiten von "twindow" die Anzahl der Beobachtungen und keine Zeiteinheit sind? Vielen Dank! – user3004015