spline.krige: Kriging combined with splines

In GSIF: Global Soil Information Facilities

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.

Note

This function adjusts grid density (prediction locations) in reference to the actual local sampling intensity. High resolution grids are created where sampling density is higher and vice versa (Hengl, 2006). Low resolution grids (due to sparse data) are then downscaled to the target resolution using spline interpolation. This allows for speeding up the kriging with minimal loss in precision, whilst reducing generation of artifacts. Spline interpolation is implemented via the SAGA GIS v2.1 function "Multilevel B-Spline Interpolation" using the default settings. This function is especially suitable for producing predictions for large grids where the sampling locations show high spatial clustering. It is NOT intended for predicting using point samples collected using sampling designs with constant spatial sampling intensity e.g. point samples collected using simple random sampling or grid sampling.

Author(s)

Tomislav Hengl

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)") 
}

## End(Not run)
## conclusion: spline.krige produces less artifacts, 
## and is at order of magnitude faster than simple 'krige'

GSIF documentation built on May 2, 2019, 5:44 p.m.