Using Intamap

The INTAMAP System

Examples using the intamap package

For more details, see the section on R-package installation. Once installed the package can be loaded into a session by: > library(intamap) This also loads all the packages that the intamap package is dependent on. Information about the functionality is possible by typing a question mark and the name of the function or data set you are interested in. Additionally there is a documentation object describing the use of the package. > package?intamap
> ?interpolate
 Most of the examples below uses the Meuse data set, which consists of 155 observations of heavy metal concentrations together with soil properties and distance to the river, and 3103 prediction locations in the object meuse.grid, both from the gstat-package. For more description, type: > ?meuse

Example 1 - predicting means and kriging error

In the example below, we are predicting the logarithms of the zinc concentrations of the meuse data set at the locations in meuse.grid. We are interested in the predicted mean and prediction variance of the log(zinc) values, in addition to the probability of the value being above 7 (the exceedance probability). We let the system choose the method itself, i.e. the methodName is set to automatic. There is still a time limitation though, the method has to be estimated to finish within the maximumTime of 30 seconds.

library(intamap)
data(meuse)
meuse$value = log(meuse$zinc)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid)=~x+y

spplot(meuse,"value",col.regions = bpy.colors())

output = interpolate(observations = meuse,
predictionLocations = meuse.grid,
outputWhat = list(mean=T, variance=T, excprob = 7),
maximumTime = 30,
methodName = "automatic",
optList = list() )

plot(output)
output$variogramModel
summary(output)
summary(t(output$outputTable))

The last three commands are there to check the result. First the result is plotted, the fitted variogram model can be examined, and there is a summary function that can give an overview of the final object. This object contains observations, predictions and some of the temporary results. The simplest predictions are mostly stored in the element called vardictions, whereas the more complicated results are stored in the element outputTable. This is particularly designed for the WPS, and needs to be transposed for interactive use.

Most of the arguments above are default values, and therefore not necessary to pass to the system. By calling: output = interpolate(meuse, meuse.grid) we only miss the probability of the value being above 7. The methodName can either be set to automatic, which gives the (presumebly) best interpolation method for the current problem, given the time available set with maximumTime, or it is possible to call one of the many methods directly. Possible options are:

optList is a parameter list where it is possible to give different arguments to the interpolation routine. Some of these are:

outputWhat is a list of the wanted prediction types. These include:

Example 2 - predicting means and kriging error using copula-method

 In the example below, we remove some of the observation locations to reduce computation time. Using the same maximumTime, the system will now choose the copula methods: observations = meuse[rnorm(155)>0.7,]
predictionLocations = coarsenGrid(meuse.grid,5)
output2 = interpolate(observations = observations,
predictionLocations = predictionLocations,
outputWhat = list(mean=T, variance=T, excprob = 7),
maximumTime = 30,
methodName = "automatic",
optList = list() )

plot(output2)

Example 3 - detailed calls

 The function interpolate is a wrapper around a series of S3 class functions in R. Among these are: preProcess
estimateParameters
spatialPredict
postProcess

It is also possible to call the single functions that are called from the interpolate function. First we have to create the interpolation object, and then do some sanity checks on this object (check if there are enough observations for interpolation, if some observation locations are duplicated, if the observations and prediction locations appear to be at different locations etc): krigingObject = createIntamapObject(
observations = meuse,
predictionLocations = meuse.grid,
formulaString = "value~1",
outputWhat = list(mean=TRUE, variance=TRUE, excprob = 7),
params = list(confProj = TRUE),
class = "automap")
checkSetup(krigingObject)
 This function takes many of the same arguments as the interpolate function. The optList above is here referred to as params. methodName is here referred to as class. This class is the one that helps R choose between the different methods for each of the S3 functions in the list above. The command methods() give an overview over methods for a certain function.

> methods(estimateParameters)
[1] estimateParameters.automap* estimateParameters.copula*
[3] estimateParameters.default* estimateParameters.idw*
[5] estimateParameters.linearVariogram* estimateParameters.transGaussian*
[7] estimateParameters.yamamoto*

Non-visible functions are asterisked
The functions that does the interpolation can then be called: > krigingObject = preProcess(krigingObject)
> krigingObject = estimateParameters(krigingObject)
> krigingObject = spatialPredict(krigingObject)
> krigingObject = postProcess(krigingObject)
> plot(krigingObject)

Example 4 - block kriging

 The intamap package is also able to perform block kriging. This can be done through the function interpolateBlock, either on a grid, assuming that the block size is equal to the grid resolution, or on a set of predefined polygons. For block kriging on a grid, it is particularly easy to see the difference of the prediction variance. predictionLocations = spsample(meuse.grid,50,"regular")
gridded(predictionLocations) = TRUE
meuse$value = log(meuse$zinc)

outputWhat = list(mean=TRUE,variance=TRUE,excprob = 7)
outputBlock = interpolateBlock(meuse,predictionLocations,outputWhat, methodName = "automap")

plot(outputBlock)

Example 5 - block kriging - other prediction types

 In addition to the prediction types possible for point predictions, it is also possible to request prediction types that are only sensible for block predictions, such as the fraction above threshold or the predicted maximum within a block. These estimates are either based on analytial predictions (if possible) or from simulations if necessary. outputWhat = list(mean = TRUE, variance = TRUE)
blockWhat = list(fat = 7)

outputBlock2 = interpolateBlock(meuse,predictionLocations,outputWhat = outputWhat,
blockWhat = blockWhat, methodName = "automap")
plot(outputBlock2)