Using Intamap

The INTAMAP System

Examples using the intamapInteractive package

The package is on CRAN and it can be most easily downloaded in an R-session with: > install.packages("intamapInteractive") The R-package installation gives more details. The package can be loaded into a session by: > library(intamapInteractive) This also loads all the packages that the intamapInteractive package is dependent on, including the intamap package. 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, which can be accessed by the first command below. > package?intamapInteractive
 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 The other data set used is the sic2004 data set, also available from the gstat-package.

Example 1 - Bias identification

Many data sets can consist of data coming from a large number of different measurement networks, using different measurement devices or applying different methods for post-processing the observations. Some of these networks can exist in the same area, e.g. when different authorities are measuring the same, but at different locations (one of them in cities, the other one close to lakes), some networks will only exist as neighbouring networks (networks operated by a municipality or a country). Local networks can also be grouped together as one national data-base, which can again be merged into an international data-base.

The example below is a hypothetical example. In the first part of the example, the bounding box of the Meuse data set is divided into 12 regions, 8 of them with observations, defined in the object pBoundaries. The estimated bias parameters in object$regionalBias$regionalBias reflects that this method is might find biases where these should not be present, the bias for the North-western and South-eastern part are both in the order of +/- 0.5. This is one of the reasons why such a method is not included as a part of the automatic interpolation service in the intamap package, the results should be checked by the user.

library(intamapInteractive)
data(meuse)
data(meuse.grid)
observations = data.frame(x = meuse$x,y = meuse$y,value = log(meuse$zinc))
coordinates(observations) = ~x+y
gridded(meuse.grid) = ~x+y
pBoundaries = spsample(observations, 10, "regular",bb = bbox(observations) +
matrix(c(-400,-400,400,400),ncol=2),offset=c(0,0))
gridded(pBoundaries) = TRUE
cs = pBoundaries@grid@cellsize[1]/2

Srl = list()
nb = dim(coordinates(pBoundaries))[1]
for (i in 1:nb) {
pt1 = coordinates(pBoundaries)[i,]
x1 = pt1[1]-cs
x2 = pt1[1]+cs
y1 = pt1[2]-cs
y2 = pt1[2]+cs

boun = data.frame(x=c(x1,x2,x2,x1,x1),y=c(y1,y1,y2,y2,y1))
coordinates(boun) = ~x+y
boun = Polygon(boun)
Srl[[i]] = Polygons(list(boun),ID = as.character(i))
}
pBoundaries = SpatialPolygonsDataFrame(SpatialPolygons(Srl),
data = data.frame(ID=c(1:nb)))
observations$ID = overlay(observations,pBoundaries)
blines = findBoundaryLines(pBoundaries,regCode = "ID")

object = createIntamapObject(observations,meuse.grid,boundaryLines = blines,
params = list(removeBias = "regionalBias"))
object = biasCorr(object,regCode= "ID")
object$regionalBias$regionalBias

Example 2 - Optimization

It is in many cases necessary to change an existing measurement network. Either it is possible to install more there are more probes available, or it is necessary to remove observations location. This is possible using the function optimizeNetwork. There are several different methods available, and possible to either add or delete observation locations in a measurement network.

library(intamapInteractive)
library(maptools)
# use SIC2004 dataset
data(sic2004)

# load the input for the observation and prediction locations
# including the covariates for UK, and changing to more intuitive
# variable names
coordinates(sic.val) = ~x+y
observations = sic.val["dayx"]

#Make a polygon for the candidate locations
coordinates(sic.grid)=~x+y
predGrid = sic.grid
bb = bbox(predGrid)
boun = SpatialPoints(data.frame(x=c(bb[1,1],bb[1,2],bb[1,2],bb[1,1],bb[1,1]),
y=c(bb[2,1],bb[2,1],bb[2,2],bb[2,2],bb[2,1])))
Srl = Polygons(list(Polygon(boun)),ID = as.character(1))
candidates = SpatialPolygonsDataFrame(SpatialPolygons(list(Srl)),
data = data.frame(ID=1))

# limit the number of prediction locations to have faster UK
# computations
nGrid = dim(coordinates(predGrid))[1]
predGrid = predGrid[sample(seq(1,nGrid),1000),]
# fit the variogram model (using function fit.variogram from package
# gstat)
model = fit.variogram(variogram(dayx~x+y, observations), vgm(50, "Sph", 250000, 250))
plot(variogram(dayx~x+y, observations), model=model, formulaString = dayx~x+y)
# compute the Mukv of the current network
initMukv <- calculateMukv(observations, predGrid, model)
print(initMukv)
# delete optimally 10 stations from current network with method "ssa"
# (spatial simulated annealing) and criterion "mukv"
optim1 = optimizeNetwork(observations, predGrid, candidates, method = "ssa",
action = "del", nDiff = 10, model = model, criterion = "MUKV", plot = TRUE )
MukvDel1 <- calculateMukv(optim1, predGrid, model, formulaString = dayx~x+y)
print(MukvDel1)
# add optimally 10 stations from current network with method "spcov"
# (spatial coverage method)
optim2 = optimizeNetwork(observations, predGrid, candidates, method = "ssa",
action = "add", nDiff = 10, model = model, criterion = "MUKV", plot = TRUE )
MukvAdd1 <- calculateMukv(optim2, predGrid, model, formulaString = dayx~x+y)
print(MukvAdd1)

# delete optimally 10 stations from current network with method "manual"
optim3 = optimizeNetwork(observations, predGrid, candidates, method = "manual",
action = "del", nDiff = 10, model = model, criterion = "MUKV", plot = TRUE )
MukvDel2 <- calculateMukv(optim3, predGrid, model, formulaString = dayx~x+y)
print(MukvDel2)
print(MukvDel1)
print(MukvAdd1)
print(initMukv)