GSIF (version 0.5-5.1)

spline.krige: Kriging combined with splines

Description

Combines kriging and spline interpolation to speed up the kriging with minimal loss in precision, whilst reducing generation of artifacts. Spline interpolation is implemented via the SAGA GIS function "Multilevel B-Spline Interpolation" (SAGA GIS needs to be installed separately).

Usage

spline.krige(formula, locations, newdata, newlocs = NULL, model, 
    te = as.vector(newdata@bbox), file.name, silent = FALSE, 
    t_cellsize = newdata@grid@cellsize[1], optN = 20, quant.nndist = .5, 
    nmax = 30, predictOnly = FALSE, resample = TRUE, saga.env, 
    saga.lib=c("grid_spline","grid_tools"), saga.module=c(4,0), …)

Arguments

formula

formula that defines the dependent variable as a linear model of independent variables; usually in the form z~1

locations

object of class SpatialPoints; sampling locations

newdata

object of class SpatialPixels*; spatial domain of interest

newlocs

object of class SpatialPoints*; prediction locations produced using the resample.grid function (if missing it will be generated using the resample.grid function)

model

variogram model of dependent variable (or its residuals); see gstat::krige

te

numeric; a vector in the form c(xmin,ymin,xmax,ymax); sets bounding box of the kriging predictions

file.name

character; optional output file name pattern (without any file extension)

silent

logical; specifies whether to print out the progress

t_cellsize

numeric; target cell size (output grid)

optN

integer; optimal number of prediction locations per sampling location e.g. 1 sampling location is used to predict values for 20 new pixels

quant.nndist

numeric; threshold probability to determine the search radius (sigma)

nmax

integer; the number of nearest observations that should be used for kriging

predictOnly

logical; specifies whether to generate only predictions (var1.pred column)

resample

logical; specifies whether to down or upscale SAGA GIS grids to match the grid system of newdata

saga.env

list; path to location of the SAGA binaries (extracted using rsaga.env())

saga.lib

character; names of the SAGA libraries used

saga.module

integer; corresponding module numbers

…

other optional arguments that can be passed to function gstat::krige

Value

Returns an object of class "SpatialGridDataFrame", or an output file name.

References

• Hengl T., (2006) Finding the right pixel size. Computers and Geosciences, 32(9): 1283-1298.

• SAGA GIS (http://sourceforge.net/projects/saga-gis/)

• SpatStat package (https://cran.r-project.org/package=spatstat)

Examples

# NOT RUN {
library(plotKML)
library(spatstat)
library(RSAGA)
library(gstat)
library(raster)
data(eberg)
data(eberg_grid)
data(eberg_grid25)
library(sp)
coordinates(eberg) <- ~X+Y
proj4string(eberg) <- CRS("+init=epsg:31467")
m <- vgm(psill=320, model="Exp", range=1200, nugget=160)
plot(variogram(SNDMHT_A~1, eberg[!is.na(eberg$SNDMHT_A),]), m)
## prediction locations:
gridded(eberg_grid) <- ~x+y
proj4string(eberg_grid) <- CRS("+init=epsg:31467")
gridded(eberg_grid25) <- ~x+y
proj4string(eberg_grid25) <- CRS("+init=epsg:31467")
## prepare prediction locations for spline.krige:
grd <- resample.grid(locations=eberg["SNDMHT_A"], t_cellsize=25,
   newdata=eberg_grid25, optN=5, quant.nndist=.9)
## plot resampled grid:
plot(raster(grd$density))
plot(grd$newlocs)
points(eberg, pch=19, col="red", cex=.7)
env <- rsaga.env()
if(exists("env") & env$version=="2.1.0"){
 ## compare processing time:
 system.time( SND.sok <- spline.krige(locations=eberg["SNDMHT_A"], 
      t_cellsize=25, newdata=eberg_grid25, 
      newlocs=grd$newlocs, model=m, nmax=30) )
 system.time( SND.ok <- krige(SNDMHT_A~1, 
      eberg[!is.na(eberg$SNDMHT_A),], 
      newdata=eberg_grid, m, 
      debug.level = -1, nmax=30) )
 system.time( SND.ok25 <- krige(SNDMHT_A~1, 
      eberg[!is.na(eberg$SNDMHT_A),], 
      newdata=eberg_grid25, m, 
      debug.level = -1, nmax=30) )  
 ## compare outputs visually:
 par(mfrow=c(1,3))
 plot(raster(SND.sok[1]), main="spline.krige (25 m)")
 plot(raster(SND.ok25[1]), main="krige (25 m)")
 plot(raster(SND.ok[1]), main="krige (100 m)") 
}
# }
# NOT RUN {
## conclusion: spline.krige produces less artifacts, 
## and is at order of magnitude faster than simple 'krige'
# }