This package provides functionality for automatic interpolation of spatial data. The package was originally developed as the computational back-end of the intamap web service, but is now a stand-alone package as maintenance of the web service has ended.
The normal work flow for working with the
intamap package can best be illustrated
with the following R-script. The procedure starts with reading data and meta data,
then setting up an object which is used in the following functions: preprocess data,
estimate parameters, compute spatial predictions, and post process them
(i.e., write them out):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
library(intamap) # set up intamap object, either manually: obj = list( observations = readOGR("PG:dbname=postgis", "eurdep.data"), predictionLocations = readOGR("PG:dbname=postgis", "eurdep1km.grid"), targetCRS = "+init=epsg:3035", params = getIntamapParams() ) class(obj) = c("idw") # or using createIntamapObject obj = createIntamapObject( observations = readOGR("PG:dbname=postgis", "eurdep.data"), predictionLocations = readOGR("PG:dbname=postgis", "eurdep1km.grid"), targetCRS = "+init=epsg:3035",class = c("idw") ) # run test: checkSetup(obj) # do interpolation steps: obj = preProcess(obj) obj = estimateParameters(obj) # faster obj = spatialPredict(obj) obj = postProcess(obj)
Our idea is that a script following this setup will allow the full statistical analysis required for the R back-end to the automatic interpolation service, and provides the means to extend the current (over-simplistic) code with the full-grown statistical analysis routines developed by INTAMAP partners. Running the package independently under R gives the user more flexibility in the utilization than what is possible through the web-interface.
Let us look into detail what the code parts do:
library(intamap) loads the R code of the
package to the current R session, along with the packages required for
this (sp, gstat, akima, automap, mvtnorm, evd, MASS). It is also recommended
to install and load rgdal on those platforms where this package is available.
All packages need to be available to the
R session, which is possible after downloading them from
the Comprehensive R Network Archives (CRAN) (https://cran.r-project.org)
1 2 3 4 5 6 7
This code sets up a list object called
obj, and assigns a class
(or a group of classes) to it. This list should hold anything we need in the next steps, and the
bare minimum seems to be measured point data (which will be extended to
polygon data) and prediction locations,
and a suggestion what to do with it. Here, the data are read from a
PostGIS data base running on localhost; data base connections over a
network are equally simple to set up. From the data base
inspire1km.grid are read; it
is assumed that these have their SRID (spatial reference identifier) set.
The suggestion what to do with these data is put in the classes,
idw. This will determine which versions of
parameterEstimate etc will be used:
intamap provides methods
for each of the generic functions
Although it would be possible to apply two classes in this case (
dataType in addition to
as the choice of pre- and post-processing steps
tend to be data-dependent, we have tried to limit the number of classes to one for most applications.
The S3 method mechanism (used here) hence requires these versions to
postProcess.idw (and eventually
To see that, we get in an interactive session:
1 2 3 4 5 6 7 8 9 10
> library(intamap) Loading required package: sp Loading required package: gstat Loading required package: rgdal Geospatial Data Abstraction Library extensions to R successfully loaded > methods(estimateParameters)  estimateParameters.automap* estimateParameters.copula*  estimateParameters.default* estimateParameters.idw*  estimateParameters.linearVariogram* estimateParameters.transGaussian*  estimateParameters.yamamoto*
Now if a partner provides additional methods for BayesianKriging, one could integrate them by
class(obj) = "BayesianKriging"
and provide some or all of the functions
postProcess.BayesianKriging, which would be called automatically
when using their generic form (
It is also possible to provide a method that calls another
method. Further, for each generic there is a default method. For
spatialPredict these print an
error message and stop, for the pre- and postprocessing the default
methods may be the only thing needed for the full procedure; if no
preProcess.BayesianKriging is found,
will be used when the generic (
preProcess) is called.
If a method does something, then it adds its result to the object it received, and returns this object. If it doesn't do anything, then it just passes (returns) the object it received.
To make these different methods exchangable, it is needed that they can all make the same assumptions about the contents of the object that they receive when called, and that what they return complies with what the consequent procedures expect. The details about that are given in the descriptions of the respective methods, below.
Because a specific interpolation method implemented may have its peculiar
characteristics, it may have to extend these prescriptions by passing
more information than described below, for example information about
The choice between methods is usually done based on the type of problem
(extreme values present, computation time available etc.). The possibility
for parallel processing of the prediction step is enabled for some of the main methods.
To be able to take advantage of multiple CPUs on a computer, the package
doParallel must be installed, additionally the parameter nclus must be set to
a value larger than 1.
object of class \
value that is the target variable.
object extending class
character; target CRS or missing
formula string for parameter estimation and prediction functions
list of parameters, to be set in
getIntamapParams. These parameters include:
Defining whether anisotropy should be calculated
Defining whether biases should be removed, and in case yes, which ones
Defining which biases to be added in the
This has not yet been implemented.
character; specifies which methods to use to remove bias. See below.
Defining if the predictions should be subject to segmentation. Segmentation has been implemented, but not the use of it.
for local kriging: the number of nearest observations that should be used for a kriging prediction or simulation, where nearest is defined in terms of the space of the spatial locations. By default, 50 observations are used.
The number of grid points to be used if an Averaged Cumulative Distribution Function (ACDF) needs to be computed for unbiased kriging
Number of simulations when needed
Block size; a vector with 1, 2 or 3 values containing the size
of a rectangular in x-, y- and z-dimension respectively
(0 if not set), or a data frame with 1, 2 or 3 columns,
containing the points that discretize the block in the
x-, y- and z-dimension to define irregular blocks relative to
(0,0) or (0,0,0) - see also the details section of
By default, predictions or simulations refer to the support of the data values.
If known - the distribution of the data. Defaults to gaussian, analytical solutions also exists in some cases for logNormal. This setting only affects a limited number of methods, e.g. the block prediciton
If set, the program will attempt conform projections in
calling the function
The number of clusters to use, if applying to a method which can
run processes in parallel. Currently implemented for methods
Used in some functions for giving additional output. See individual functions for more information.
Additional parameters that do not exist in the default parameter set,
particularly parameters necessary for new methods within the
Pebesma, E., Cornford, D., Dubois, G., Heuvelink, G.B.M., Hristopulos, D., Pilz, J., Stohlker, U., Morin, G., Skoien, J.O. INTAMAP: The design and implementation f an interoperable automated interpolation Web Service. Computers and Geosciences 37 (3), 2011.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.