Using Intamap
The INTAMAP System
Examples using the psgp package
The following example illustrates the use of psgp in R for the case of observations with noise measurements. We show how noise models can be defined (in term of noise variance) and attributed to each observation. For this particular case, we assume there is a different noise model for each observations, although in practice several observations could use the same noise model.
We start by loading the Meuse dataset and project it onto CRS coordinates. # Load our favourite data set data(meuse) coordinates(meuse) = ~x+y meuse$value = log(meuse$zinc) data(meuse.grid) gridded(meuse.grid) = ~x+y proj4string(meuse) = CRS("+init=epsg:28992") proj4string(meuse.grid) = CRS("+init=epsg:28992") We then indicate that each observation will use its own noise model. This is done by setting a meuse$oeid vector of observation identifiers. Each identifier must be the index of a given noise model (as defined below). # Specify a different observation error model for each observation nobs = length(meuse$value) # Number of observations meuse$oeid = seq(1:nobs) # One error model per observation
We now define the noise models, in this case one for each observation. Each noise model is determined by its noise variance (the error is assumed Gaussian with zero mean). The variances must be specified in the vector meuse$oevar, where the i-th element is the variance for the error model i. Here, the variance are generated at random. # Indicate the variance for each of these error models meuse$oevar <- abs( rnorm( max(meuse$oeid) ) ) One must ensure that the indexes in meuse$oeid do not exceed the number of error models in meuse$oevar.
To specify that observation j will use error model i, we would simply use: meuse$oeid[j] = i
We can then create the intamap object, passing "psgp" as the object class. # Set up intamap object obj = createIntamapObject( observations = meuse, predictionLocations = meuse.grid, targetCRS = "+init=epsg:3035", class = "psgp" ) The optimal psgp parameters are determined using the estimateParameters() method, and predictions at a set of unobserved locations are done using the spatialPredict() method. Both methods require the intamap object as an argument. # Estimate parameters and predict at new locations (interpolation) obj = conformProjections(obj) obj = estimateParameters(obj) obj = spatialPredict(obj)
Last, we plot the interpolation results. # Plot results plotIntamap(obj)
