| Title: | Procedures for Automated Interpolation |
|---|---|
| Description: | Geostatistical interpolation has traditionally been done by manually fitting a variogram and then interpolating. Here, we introduce classes and methods that can do this interpolation automatically. Pebesma et al (2010) gives an overview of the methods behind and possible usage <doi:10.1016/j.cageo.2010.03.019>. |
| Authors: | Edzer Pebesma [aut], Jon Olav Skoien [aut, cre], Olivier Baume [ctb], A Chorti [ctb], Dionisis Hristopulos [ctb], Hannes Kazianka [ctb], Stepahnie Melles [ctb], Giannis Spiliopoulos [ctb] |
| Maintainer: | Jon Olav Skoien <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.5-7 |
| Built: | 2024-10-24 06:53:03 UTC |
| Source: | CRAN |
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):
library(intamap)
library(sf)
# set up intamap object, either manually:
data(meuse)
coordinates(meuse) = ~x + y
data(meuse.grid)
gridded(meuse.grid) = ~x + y
obj = list(
observations = meuse,
predictionLocations = meuse.grid,
targetCRS = "+init=epsg:3035",
params = getIntamapParams()
)
class(obj) = c("idw")
# or using createIntamapObject
obj = createIntamapObject(
observations = meuse,
predictionLocations = meuse.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)
The command library(intamap) loads the R code of the intamap
package to the current R session, along with the packages required for
this (sp, gstat, akima, automap, mvtnorm, evd, MASS). After the retirement
of rgdal, it is recommended to use sf-objects from the sf package.
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)
# set up intamap object:
obj = createIntamapObject(
observations = meuse,
predictionLocations = meuse.grid,
targetCRS = "+init=epsg:3051",
class = "idw"
)
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 postgis
the tables eurdep.data and 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 preProcess,
parameterEstimate etc will be used: intamap provides methods
for each of the generic functions
preProcess,
estimateParameters,
spatialPredict,
postProcess.
Although it would be possible to apply two classes in this case (dataType in addition to idw),
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
be called preProcess.idw, estimateParameters.idw,
spatialPredict.idw, and postProcess.idw (and eventually
also preProcess.eurdep and preProcess.eurdep).
To see that, we get in an interactive session:
> library(intamap) Loading required package: sp Loading required package: gstat Geospatial Data Abstraction Library extensions to R successfully loaded > methods(estimateParameters) [1] estimateParameters.automap* estimateParameters.copula* [3] estimateParameters.default* estimateParameters.idw* [5] estimateParameters.linearVariogram* estimateParameters.transGaussian* [7] 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
preProcess.BayesianKriging,
estimateParameters.BayesianKriging,
spatialPredict.BayesianKriging, and
postProcess.BayesianKriging, which would be called automatically
when using their generic form (preProcess etc).
It is also possible to provide a method that calls another
method. Further, for each generic there is a default method. For
estimateParameter and 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, preProcess.default
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
priors from estimateParameters to spatialPredict.
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.
observations object of class
SpatialPointsDataFrame, containing
a field value that is the target variable.
predictionLocations object extending class Spatial, containing
prediction locations.
targetCRScharacter; target CRS or missing
formulaStringformula string for parameter estimation and prediction functions
paramslist 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
(localBias and regionalBias implemented
Defining which biases to be added in the postProcess function.
This has not yet been implemented.
"LM"
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 predict.gstat.
By default, predictions or simulations refer to the support of the data values.
"gaussian" 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 preProcess,
calling the function conformProjections.
The number of clusters to use, if applying to a method which can
run processes in parallel. Currently implemented for methods automap,
copula and psgp.
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 intamap
package
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.
Calculates predictive mean, predictive variance, predictive quantiles and exceedance probabilities for certain thresholds in the spatial copula model.
bayesCopula(obj,estimates,search=10,calc=list(mean=TRUE,variance=TRUE),testMean=FALSE)bayesCopula(obj,estimates,search=10,calc=list(mean=TRUE,variance=TRUE),testMean=FALSE)
obj |
Intamap object including observations and predictionLocations, see
|
estimates |
List of estimated parameters (typically obtained by calling |
search |
local prediction: number of observed locations considered for prediction at each unknown point |
calc |
list of what prediction type is required:
|
testMean |
Whether or not the predictive means (if calculated) should be tested for being reasonable. |
bayesCopula is used for plug-in prediction at unobserved spatial locations. The name of the function is somewhat
misleading since no Bayesian approach is implemented so far. It is possible to calculate numerically the predictive mean
and variance for both the Gaussian and the chi-square spatial copula model. Exceedance probabilities and predictive
quantiles are only supported for the Gaussian copula model. Note that it may occur that the predictive distribution has
no finite moments. In this case, a possible predictor is the median of the predictive distribution. If testMean=TRUE and
the predictive means have no reasonable values, the median is automatically calculated and a warning is produced.
The copula prediction method is computationally demanding.
There is a possibility of running it as a parallel process by setting the parameter
nclus > 1 for the interpolation process. This requires a previous installation
of the package doParallel.
List with the following elements:
Mean of the predictive distribution. NULL if not calculated.
Variance of the predtictive distribution. NULL if not calculated.
Quantiles of the predictive distribution NULL if not calculated.
Probabilities for the predictive distribution to exceed predefined thresholds. NULL if not calculated.
Hannes Kazianka
[1] Kazianka, H. and Pilz, J. (2009), Spatial Interpolation Using Copula-Based Geostatistical Models. GeoENV2008 - Geostatistics for Environmental Application (P. Atkinson, C. Lloyd, eds.), Springer, New York
[2] 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.
copulaEstimation, spatialPredict, estimateParameters
## Not run: data(intamapExampleObject) ## estimate parameters for the copula model copula <- list(method="norm") anisotropy <- list(lower = c(0,1), upper = c(pi, Inf), params = c(pi/3, 2)) correlation <- list(model = "Ste", lower=c(0.01, 0.01, 0.01), upper = c(0.99, Inf, 20), params = c(0.05, 4, 3)) margin <- list(name = "gev", lower = c(0.01, -Inf), upper = c(Inf, Inf), params = c(30, 0.5)) trend <- list(F = as.matrix(rep(1, 196)), lower = -Inf, upper = Inf, params = 40) estimates <- copulaEstimation(intamapExampleObject, margin, trend, correlation, anisotropy, copula) ## make predictions at unobserved locations predictions<-bayesCopula(intamapExampleObject, estimates, search = 25, calc = list(mean = TRUE, variance = TRUE, excprob = 40, quantile = 0.95)) ## End(Not run)## Not run: data(intamapExampleObject) ## estimate parameters for the copula model copula <- list(method="norm") anisotropy <- list(lower = c(0,1), upper = c(pi, Inf), params = c(pi/3, 2)) correlation <- list(model = "Ste", lower=c(0.01, 0.01, 0.01), upper = c(0.99, Inf, 20), params = c(0.05, 4, 3)) margin <- list(name = "gev", lower = c(0.01, -Inf), upper = c(Inf, Inf), params = c(30, 0.5)) trend <- list(F = as.matrix(rep(1, 196)), lower = -Inf, upper = Inf, params = 40) estimates <- copulaEstimation(intamapExampleObject, margin, trend, correlation, anisotropy, copula) ## make predictions at unobserved locations predictions<-bayesCopula(intamapExampleObject, estimates, search = 25, calc = list(mean = TRUE, variance = TRUE, excprob = 40, quantile = 0.95)) ## End(Not run)
blockPredict is a generic method for prediction of
spatially aggregated variables within the intamap-package package.
blockPredict(object, ...)blockPredict(object, ...)
object |
a list object of the type described in |
... |
other arguments that will be passed to the requested interpolation method.
See the individual interpolation methods for more information. The following arguments
from
|
The function blockPredict is a wrapper around the spatialPredict.block function
within the intamap-package package, to simplify the calls for block predictions.
Block predictions are spatial predictions assumed to be valid for a certain area.
The blocks can either be given by passing SpatialPolygons as the
predicitonLocations or by passing the block-argument through the parameters of the
object or through the ...-argument.
There are esentially two ways to solve the problems of block predictions.
block predictions can be found directly by block kriging
block predictions can be found through numerical simulations over a set of points within the block, the requested output is found by averaging over these simulations
The analytical solutions are used when applicable. This is typically for ordinary kriging based methods and prediction types that can be found by linear aggregation (e.g. block mean).
If the prediction type necessitates simulations, this is done by subsampling the blocks. This can either be done block-wise, with a certain number of points within each block, with a certain cellsize, or with a certain number of points
automap Uses function autoKrige in the
automap package.
If object already includes a variogram model,
krige in the gstat-package will be called directly.
a list object similar to object, but extended with predictions at
a the set of locations defined object.
Jon Olav Skoien
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.
# This example skips some steps that might be necessary for more complicated # tasks, such as estimateParameters and pre- and postProcessing of the data 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") # set up intamap object: obj = createIntamapObject( observations = meuse, predictionLocations = meuse.grid[sample(1:length(meuse.grid),10),], targetCRS = "+init=epsg:3035", class = "automap" ) # do interpolation step: obj = conformProjections(obj) obj = estimateParameters(obj) obj = blockPredict(obj,block=c(100,100)) # blockPredict # intamap object for which simulation is needed: meuse$value = meuse$zinc obj = createIntamapObject( observations = meuse, predictionLocations = meuse.grid[sample(1:length(meuse.grid),5),], params = list(ngrid = 16), class = "transGaussian" # trans-Gaussian kriging method ) obj = estimateParameters(obj, lambda = 0) # lambda is optional, lambda = 0 gives lognormal kriging obj = blockPredict(obj,block=c(100,100)) # blockPredict# This example skips some steps that might be necessary for more complicated # tasks, such as estimateParameters and pre- and postProcessing of the data 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") # set up intamap object: obj = createIntamapObject( observations = meuse, predictionLocations = meuse.grid[sample(1:length(meuse.grid),10),], targetCRS = "+init=epsg:3035", class = "automap" ) # do interpolation step: obj = conformProjections(obj) obj = estimateParameters(obj) obj = blockPredict(obj,block=c(100,100)) # blockPredict # intamap object for which simulation is needed: meuse$value = meuse$zinc obj = createIntamapObject( observations = meuse, predictionLocations = meuse.grid[sample(1:length(meuse.grid),5),], params = list(ngrid = 16), class = "transGaussian" # trans-Gaussian kriging method ) obj = estimateParameters(obj, lambda = 0) # lambda is optional, lambda = 0 gives lognormal kriging obj = blockPredict(obj,block=c(100,100)) # blockPredict
checkSetup will do some sanity checks on input data provided through object.
checkSetup(object, quiet = FALSE)checkSetup(object, quiet = FALSE)
object |
object, to be passed to |
quiet |
logical; TRUE to suppress OK statement |
checkSetup is a function that makes certain tests on the intamap object to
make sure that it is suited for interpolation. Particularly, it will issue a warning
or an error if one of the following conditions are met:
observations is not an element of object
observations contain less than 20 observations
Some of the observation locations are duplicated
formulaString is not an element of object
None of the columns of observations has a name that corresponds to the independent variable of formulaString
predictionLocations is not an element of object
predictionLocations is not a Spatial object
targetCRS is given, but observations and predictionLocations do not have CRS set
addBias includes biases that are not part of removeBias
The function will issue a warning if it appears that predictionLocations
and observations share a small region. This warning is given as it is
a likely cause of errors, although it can also happen if predictionLocations
are limited to one small cluster.
returns TRUE if check passes, will halt with error when some some error condition is met.
Edzer J. Pebesma
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.
coarsenGrid is a function that resamples a SpatialGridDataFrame.
coarsenGrid(object,coarse=2,offset = sample(c(0:(coarse-1)),2,replace=TRUE))coarsenGrid(object,coarse=2,offset = sample(c(0:(coarse-1)),2,replace=TRUE))
object |
a |
coarse |
an integer telling how much the grid should be coarsened |
offset |
integer giving the relative offset of the first point, see details below for a closer description |
The function coarsenGrid is a function that samples from a
SpatialGridDataFrame.
The argument coarse indicates that every coarse row and column will
be sampled, starting with the row and column represented in offset. offset = c(0,0)
implies that the smallest x- and y-coordinates will be a part of the resampled
data set, offset = c(1,1) implies that sampling will start on the second row and column.
Jon Olav Skoien
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 of an interoperable automated interpolation Web Service. Computers and Geosciences 37 (3), 2011.
data(meuse.grid) gridded(meuse.grid) = ~x+y newMeuseGrid = coarsenGrid(meuse.grid,coarse=2,offset=c(1,1))data(meuse.grid) gridded(meuse.grid) = ~x+y newMeuseGrid = coarsenGrid(meuse.grid,coarse=2,offset=c(1,1))
Getting a conformed projection for a set of Spatial* elements
necessary for interpolation in the intamap-package.
conformProjections(object)conformProjections(object)
object |
an object of the type described in |
conformProjections is a function that attempts to reproject all projected
elements in object to one common projection.
The function is usually called with an intamap object as argument
from createIntamapObject if the parameter confProj = TRUE.
Thus it is a function that is usually not necessary to call separately.
The need for this function is because
several of the functions in a typical spatial interpolation work flow inside the
intamap-package require that the elements have a common projection. In
addition, there are some functions which are not able to deal with
unprojected spatial objects, i.e. objects with coordinates given in lattitude and
longitude. conformProjections will hence also attempt to reproject all
elements that
have coordinates in lattitude and longitude, even in the cases where they
all have the same projections.
If only one of observations or predictionLocations has a projection (or is longlat), the other one is assumed to be equal. A warning is issued in this case.
The common projection depends on the object that is passed to conformProjections.
First of all, if intCRS (see below) is present as an element of the object, all elements
will be reprojected to this projection. If not, intCRS will be set equal to
the first projection possible in the list below.
Can be given as a component in object - and is the
user-defined common projection used for interpolation
Can be given as a component in object - and is
the user-defined target projections
The projection of the predictionLocations in object
The projection of the observations
A list of the parameters to be included in the object described in intamap-package
Jon Olav Skoien
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.
data(meuse) coordinates(meuse) = ~x+y proj4string(meuse) <- CRS("+proj=stere +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m") predictionLocations = spsample(meuse, 50, "regular") krigingObject = createIntamapObject( observations = meuse, predictionLocations = predictionLocations, formulaString = as.formula("log(zinc)~1"), intCRS = "+init=epsg:3035" ) krigingObject = conformProjections(krigingObject) proj4string(meuse) proj4string(krigingObject$observations)data(meuse) coordinates(meuse) = ~x+y proj4string(meuse) <- CRS("+proj=stere +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m") predictionLocations = spsample(meuse, 50, "regular") krigingObject = createIntamapObject( observations = meuse, predictionLocations = predictionLocations, formulaString = as.formula("log(zinc)~1"), intCRS = "+init=epsg:3035" ) krigingObject = conformProjections(krigingObject) proj4string(meuse) proj4string(krigingObject$observations)
Estimates parameters of the spatial copula model using maximum likelihood.
copulaEstimation(obj,margin,trend,correlation,anisotropy,copula,tol=0.001,...)copulaEstimation(obj,margin,trend,correlation,anisotropy,copula,tol=0.001,...)
obj |
Intamap object, see description in |
margin |
list with the following elements:
|
trend |
list with the following elements:
|
correlation |
list with the following elements:
|
anisotropy |
list with the following elements:
|
copula |
list with the following elements:
|
tol |
Tolerance level for the optimization process. |
... |
Arguments to be passed to |
copulaEstimation performs maximum likelihood estimation of all possible parameters included in the Gaussian and
chi-squared spatial copula model: parameters of the predefined family of marginal distributions (including spatial trend
or external drift), correlation function parameters, parameters for geometric anisotropy and parameters for the copula
(only used for the chi-squared copula model). Due to the large number of variables that need to be optimized, a
profile-likelihood approach is used. Although convergence to a global optimum is not assured, the profile-likelihood method
makes it less likely that the optimization routine, optim, gets stuck in a local optimum. The result of
copulaEstimation is a list containing all parameter point estimates that are needed for plug-in spatial
prediction. It is advisable to check the output of the algorithm by trying different starting values for the optimization.
A list with the following elements:
margin |
Same as the input except that the list element "params" now consists of the optimized parameters of the marginal distribution function. |
trend |
Same as the input except that the list element "params" now consists of the optimized parameters of the trend model. |
correlation |
Same as the input except that the list element "params" now consists of the optimized parameters of the correlation function model. |
anisotropy |
Same as the input except that the list element "params" now consists of the optimized parameters of geometric anisotropy. |
copula |
Same as the input except that the list element "params" now consists of the optimized copula parameters. |
Hannes Kazianka
[1] Kazianka, H. and Pilz, J. (2009), Spatial Interpolation Using Copula-Based Geostatistical Models. GeoENV2008 - Geostatistics for Environmental Application (P. Atkinson, C. Lloyd, eds.), Springer, New York
[2] 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.
bayesCopula, spatialPredict, estimateParameters
data(intamapExampleObject) ## estimate parameters for the copula model ## Not run: copula<-list(method="norm") anisotropy <- list(lower = c(0, 1), upper = c(pi, Inf), params = c(pi/3, 2)) correlation <- list(model = "Ste", lower = c(0.01, 0.01, 0.01), upper = c(0.99, Inf, 20), params = c(0.05, 4, 3)) margin <- list(name = "gev", lower = c(0.01, -Inf), upper = c(Inf, Inf), params = c(30, 0.5)) trend <- list(F = as.matrix(rep(1, 196)), lower = -Inf, upper = Inf, params = 40) estimates <- copulaEstimation(intamapExampleObject, margin, trend, correlation, anisotropy, copula) ## make predictions at unobserved locations predictions <- bayescopula(intamapExampleObject, estimates, search = 25, calc = list(mean = TRUE, variance = TRUE, excprob = 40, quantile = 0.95)) ## End(Not run)data(intamapExampleObject) ## estimate parameters for the copula model ## Not run: copula<-list(method="norm") anisotropy <- list(lower = c(0, 1), upper = c(pi, Inf), params = c(pi/3, 2)) correlation <- list(model = "Ste", lower = c(0.01, 0.01, 0.01), upper = c(0.99, Inf, 20), params = c(0.05, 4, 3)) margin <- list(name = "gev", lower = c(0.01, -Inf), upper = c(Inf, Inf), params = c(30, 0.5)) trend <- list(F = as.matrix(rep(1, 196)), lower = -Inf, upper = Inf, params = 40) estimates <- copulaEstimation(intamapExampleObject, margin, trend, correlation, anisotropy, copula) ## make predictions at unobserved locations predictions <- bayescopula(intamapExampleObject, estimates, search = 25, calc = list(mean = TRUE, variance = TRUE, excprob = 40, quantile = 0.95)) ## End(Not run)
This is a help function for creating an object (see intamap-package
to be used for interpolation within the intamap-package
createIntamapObject(observations, obsChar, formulaString, predictionLocations=100, targetCRS, boundaries, boundaryLines, intCRS, params=list(), boundFile, lineFile, class="idw", outputWhat, blockWhat = "none",...)createIntamapObject(observations, obsChar, formulaString, predictionLocations=100, targetCRS, boundaries, boundaryLines, intCRS, params=list(), boundFile, lineFile, class="idw", outputWhat, blockWhat = "none",...)
observations |
a
|
obsChar |
list with observation characteristics, used by some interpolation methods |
formulaString |
formula that defines the dependent variable as a linear model
of independent variables; suppose the dependent variable has name |
predictionLocations |
either a |
targetCRS |
the wanted projection for the interpolated map |
boundaries |
|
boundaryLines |
|
intCRS |
a particular projection requested for the interpolation |
params |
parameters for the interpolation, given as exceptions to the
default parameters set in the function |
boundFile |
Filename where boundaries can be found, e.g. a shapefile |
lineFile |
Filename where paired points on boundaries can be found |
class |
setting the class(es) of the object, see |
outputWhat |
List defining the requested type of output. Parameters:
The list defaults to list (mean = TRUE) for objects of class IDW and list(mean=TRUE, variance = TRUE) for all other objects. |
blockWhat |
List defining particular output for block predictions. These include:
|
... |
|
This function is a help function for creating an object (see
intamap-package) for interpolation within the
intamap-package.
The function uses some default values if certain elements are not included.
If createIntamapObject is called without predictionLocations, or if a number
is given, the function will sample a set of predictionLocations. These will
be sampled from a regular grid.
targetCRS and intCRS are not mandatory variables, but are recommended if the user wants predictions of a certain projection. intCRS is not necessary if the targetCRS is given and has a projection (is not lat-long). It is recommended to include the argument intCRS if all projected elements are lat-long, as many of the interpolation methods do not work optimal with lat-long data.
The ...-argument can be used for arguments necessary for new methods not being
a part of the intamap-package. It is also a method for reusing previously calculated
elements that can be assumed to be unchanged for the second interpolation.
An object with observations, prediction locations, parameters and possible
additional elements for automatic interpolation. The object will have class
equal to the value of argument class, and methods in the
intamap-package will dispatch on the
object according to this class.
Jon Olav Skoien
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.
intamap-package and getIntamapParams
# set up data: 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") # set up intamap object: idwObject = createIntamapObject( observations = meuse, predictionLocations = meuse.grid, targetCRS = "+init=epsg:3035", class = "idw" )# set up data: 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") # set up intamap object: idwObject = createIntamapObject( observations = meuse, predictionLocations = meuse.grid, targetCRS = "+init=epsg:3035", class = "idw" )
This function estimates geometric anisotropy parameters for 2-D scattered data using the CTI method.
estimateAnisotropy(object,depVar, formulaString)estimateAnisotropy(object,depVar, formulaString)
object |
(i) An Intamap type object (see |
depVar |
name of the dependent variable; this is used only in case (ii). |
formulaString |
formula that defines the dependent variable as a linear model
of independent variables, only used for case (ii); suppose the dependent variable has name |
Given the input object that defines N coordinate pairs (x,y) and observed values (z),
this method estimates of the geometric anisotropy parameters.
Geometric anisotropy is a statistical property, which implies that the
iso-level contours of the covariance function are elliptical. In this case the anisotropy is determined from
the anisotropic ratio (R) and the orientation angle () of the ellipse.
Assuming a Cartesian coordinate system of axes x and y, represents the angle
between the horizontal axis and PA1, where PA1 is one of the principal axes of the ellipse,
arbitrarily selected (PA2 will denote the other axis). R represents the ratio of the correlation along PA1 divided by
the correlation length PA2. Note that the returned value of R is always greater than one (see value below.)
The estimation is based on the Covariance Tensor Identity (CTI) method. In CTI, the
Hessian matrix of the covariance function is estimated from sample derivatives. The
anisotropy parameters are estimated by explicit solutions of nonlinear equations that link
(R,) with ratios of the covariance Hessian matrix elements.
To estimate the sample derivatives from scattered data, a background square lattice is used. The lattice extends in the horizontal direction from x.min to x.max where x.min (x.max) is equal to the minimum (maximum) x-coordinate of the data, and similarly in the vertical direction. The cell step in each direction is equal to the length of the lattice to the respective direction divided by the square root of N.
BiLinear interpolation, as implemented in akima package, is used to interpolate the
field's z values at the nodes of the lattice.
The CTI method is described in detail in (Chorti and Hristopulos, 2008).
Note that to be compatible with gstat the returned estimate of the anisotropy ratio is always
greater than 1.
For observations assumed to have a trend, the trend is first subtracted from the data using universal kriging. This is an approximation, as the trend subtraction does not take anisotropy into account.
(i) If the input is an Intamap object, the value is a modification of the input object,
containing a list element anisPar with the estimated anisotropy parameters.
(ii)if the input is a SpatialPointsDataFrame, then only the list anisPar is returned.
The list anisPar contains the following elements:
ratio |
The estimate of the anisotropy ratio parameter. Using the degeneracy of the anisotropy under simultaneous ratio inversion and axis rotation transformations, the returned value of the ratio is always greater than 1. |
direction |
The estimate of the anisotropy orientation angle. It returns the angle between the major anisotropy axis and the horizontal axis, and its value is in the interval (-90,90) degrees. |
Q |
A 3x1 array containing the sample estimates of the diagonal and off-diagonal elements (Q11,Q22,Q12) of the covariance Hessian matrix evaluated at zero lag. |
doRotation |
Boolean value indicating if the estimated anisotropy is statistically significant. This value is based on a statistical test of the isotropic (R= 1) hypothesis using a non-parametric approximation for the 95 percent confidence interval for R. This approximation leads to conservative (wider than the true) estimates of the confidence interval. If doRotation==TRUE then an isotropy restoring transformation (rotation and rescaling) is performed on the coordinates. If doRotation==FALSE no action is taken. |
This function uses akima package to perform "bilinear" interpolation. The source code also allows
other interpolation methods, but this option is not available when the function is called from within INTAMAP.
In the gstat package, the anisotropy ratio is defined in the interval (0,1) and the orientation
angle is the angle between the
vertical axis and the major anisotropy axis, measured in the clockwise direction.
If one wants to use ordinary kriging inside INTAMAP
the necessary transformations are performed in the function estimateParameters.automap.
If one wants to use ordinary kriging
in the gstat package (but outside INTAMAP) the required transformations can be found in the
source code of the estimateParameters.automap
function.
A.Chorti, D.T.Hristopulos,G. Spiliopoulos
[1] 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.
[2] A. Chorti and D. T. Hristopulos (2008). Non-parametric Identification of Anisotropic (Elliptic) Correlations in Spatially Distributed Data Sets, IEEE Transactions on Signal Processing, 56(10), 4738-4751 (2008).
[3] Em.Petrakis and D. T. Hristopulos (2009). A non-parametric test of statistical isotropy for Differentiable Spatial Random Fields in Two Dimensions. Work in progress. email: [email protected]
library(gstat) data(sic2004) coordinates(sic.val)=~x+y sic.val$value=sic.val$dayx params=NULL estimateAnisotropy(sic.val,depVar = "joker")library(gstat) data(sic2004) coordinates(sic.val)=~x+y sic.val$value=sic.val$dayx params=NULL estimateAnisotropy(sic.val,depVar = "joker")
Function to estimate correlation structure parameters. The actual parameters depend on the method used.
## S3 method for class 'automap' estimateParameters(object, ... ) ## S3 method for class 'copula' estimateParameters(object, ... ) ## Default S3 method: estimateParameters(object, ...) ## S3 method for class 'idw' estimateParameters(object, ... ) ## S3 method for class 'linearVariogram' estimateParameters(object, ...) ## S3 method for class 'transGaussian' estimateParameters(object, ... ) ## S3 method for class 'yamamoto' estimateParameters(object, ... )## S3 method for class 'automap' estimateParameters(object, ... ) ## S3 method for class 'copula' estimateParameters(object, ... ) ## Default S3 method: estimateParameters(object, ...) ## S3 method for class 'idw' estimateParameters(object, ... ) ## S3 method for class 'linearVariogram' estimateParameters(object, ...) ## S3 method for class 'transGaussian' estimateParameters(object, ... ) ## S3 method for class 'yamamoto' estimateParameters(object, ... )
object |
an intamap object of the type described in |
... |
other arguments that will be passed to the requested interpolation method. See the individual methods for more information. Some parameters that are particular for some methods:
|
The function estimateParameters is a wrapper around different
methods for estimating correlation parameters to be used for the spatial
prediction method spatialPredict.
Below are some details about and/or links to the different methods currently implemented
in the intamap-package.
automapIt is possible but not necessary to estimate variogram parameters for
this method. If estimateParameters is called with an object of class automap,
autofitVariogram will be called.
If object already includes a variogram model when
spatialPredict is called,
krige in the gstat-package will be called directly.
The user can submit an argument model with the model(s) to be fitted.
copulafinding the best copula parameters using copulaEstimation
defaulta default method is not really implemented, this function is only created to give a sensible error message if the function is called with an object for which no method exist
idwfits the best possible idw-power to the data set by brute force searching within
the idpRange
linearVariogramthis function just returns the original data, no parameter fitting is necessary for linear variogram kriging
transGaussianFinding the best model parameters for transGaussian kriging
(krigeTg). This means finding the best lambda for
the boxcox-transformation and the fitted variogram
parameters for the transformed variable. If significant = TRUE
will lambda only be estimated
if the data show some deviation from normality, i.e., that at least one
of the tests described under interpolate is TRUE. Note that
transGaussian kriging is only possible for data with strictly positive values.
yamamotoa wrapper around estimateParameters.automap, only to assure that there is a method
also for this class, difference to automap is more important in spatialPredict
It is also possible to add to the above methods with functionality from other packages, if wanted. You can also check which methods are available from other packages by calling
>methods(estimateParameters)
a list object similar to object, but extended with correlation parameters.
Jon Olav Skoien
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.
createIntamapObject, spatialPredict, intamap-package
set.seed(13131) # set up data: 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") # set up intamap object: idwObject = createIntamapObject( observations = meuse, formulaString=as.formula(zinc~1), predictionLocations = meuse.grid, class = "idw" ) # run test: checkSetup(idwObject) # do interpolation steps: idwObject = estimateParameters(idwObject, idpRange = seq(0.25,2.75,.25), nfold=3) # faster idwObject$inverseDistancePowerset.seed(13131) # set up data: 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") # set up intamap object: idwObject = createIntamapObject( observations = meuse, formulaString=as.formula(zinc~1), predictionLocations = meuse.grid, class = "idw" ) # run test: checkSetup(idwObject) # do interpolation steps: idwObject = estimateParameters(idwObject, idpRange = seq(0.25,2.75,.25), nfold=3) # faster idwObject$inverseDistancePower
Function that takes time samples function that can read intamap objects
estimateTimeModel(FUN, class, formulaString, debug.level, ...)estimateTimeModel(FUN, class, formulaString, debug.level, ...)
FUN |
A string with function's name |
class |
class of intamapObject, which interpolation method to be used. |
formulaString |
the formula of the request, mainly to see if the request has independent variables. |
debug.level |
if |
... |
other arguments as defined in the createIntamapObject function |
This function uses createIntamapObject function to create synthetic data, in order to take time samples for the function with string name "FUN". The calculated model is stored, as an element of a list, in a local file (workspace for now) and it's used in order to give quick time estimates.
The function does not return a variable but stores the result in an element list with the same name.
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.
function that generates time models and saves them in workspace.
generateTimeModels(genClasses = NULL, noGenClasses = NULL, nSam = 1, test = FALSE, debug.level = 0)generateTimeModels(genClasses = NULL, noGenClasses = NULL, nSam = 1, test = FALSE, debug.level = 0)
genClasses |
list of particular classes for which time models should be generated |
noGenClasses |
list of particular classes for which time models should not be generated |
nSam |
number of attempts to be tried for each combination of predictions and observations, defaults to 1, higher number should be used for better accuracy. nSam/2 is used for copulas, to reduce computation time. |
test |
logical; if true, the time models are generated based on fewer iterations, for speed |
debug.level |
if |
This function calculates a time model for different interpolation types
in the intamap-package and returns a list object
with the estimated models. It's users responsibility to store the model in the
workspace. The normal procedure would be to run the function without arguments. However,
it is both possible to define a list for which classes the user want to
generate models, or a list of classes that are not of interest.
The time model is based on creation of a set of synthetical data sets
of different size, both regarding number of observations and prediction locations.
The function will estimate parameters and make predictions with the different combinations,
and for each method, fit a local polynomial regression model (loess)
This model can then be used by predictTime to estimate the
prediction time for an interpolation request with a certain number
of observations and prediction locations.
The function generates a timeModels object, which can be used to estimate
prediction times for different requests to the interpolate function
in the intamap-package, via predictTime.
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.
## Not run: timeModels=generateTimeModels() q("yes") ## restart R in the same directory ## End(Not run)## Not run: timeModels=generateTimeModels() q("yes") ## restart R in the same directory ## End(Not run)
This function sets a range of the parameters for the intamap-package,
to be included in the object described in intamap-package
getIntamapParams(oldPar,newPar,...)getIntamapParams(oldPar,newPar,...)
oldPar |
An existing set of parameters for the interpolation process,
of class |
newPar |
A |
... |
Individual parameters for updating
|
A list of the parameters with class intamapParams to be included in the
object described in intamap-package
This function will mainly be called by createIntamapObject, but
can also be called by the user to create a parameter set or update an
existing parameter set. If none of the arguments is a list of class
IntamapParams), the function will assume that the argument(s) are
modifications to the default set of parameters.
If the function is called with two lists of parameters (but the first one is
not of class IntamapParams) they are both seen as modifications to the
default parameter set. If they share some parameters, the parameter values from
the second list will be applied.
Jon Olav Skoien
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.
# Create a new set of intamapParameters, with default parameters: params = getIntamapParams() # Make modifications to the default list of parameters params = getIntamapParams(newPar=list(removeBias = "local", secondParameter = "second")) # Make modifications to an existing list of parameters params = getIntamapParams(oldPar = params,newPar = list(predictType = list(exc=TRUE)))# Create a new set of intamapParameters, with default parameters: params = getIntamapParams() # Make modifications to the default list of parameters params = getIntamapParams(newPar=list(removeBias = "local", secondParameter = "second")) # Make modifications to an existing list of parameters params = getIntamapParams(oldPar = params,newPar = list(predictType = list(exc=TRUE)))
get interpolation method names
getInterpolationMethodNames()getInterpolationMethodNames()
none
character array with names for available interpolation methods
Edzer Pebesma
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.
getInterpolationMethodNames()getInterpolationMethodNames()
Salted (slightly modified) observations from the European gamma radiation network
data(intamap)data(intamap)
The data set is a salted data set from the European gamma radiation network https://remap.jrc.ec.europa.eu/GammaDoseRates.aspx. Salted does in this context mean that all locations and observation values have been randomized through addition of a random value.
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.
Intamap object of class "copula" containing a simulated data set with 196 spatial locations.
data(intamapExampleObject)data(intamapExampleObject)
The data set is a realization of a random field generated using a Gaussian copula and generalized extreme value distributed margins (location=40,shape=0.5, scale=30). The correlation function is Matern (Stein's representation) with range=4, kappa=3 and nugget=0.05. Furthermore, there is geometric anisotropy with direction=pi/3 and ratio=2.
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.
spatialPredict.copula, estimateParameters.copula
## Not run: data(intamapExampleObject) ## estimate parameters for the copula model intamapExampleObject<-estimateParameters(intamapExampleObject) ## make predictions at unobserved locations intamapExampleObject<-spatialCopula(intamapExampleObject) ## End(Not run)## Not run: data(intamapExampleObject) ## estimate parameters for the copula model intamapExampleObject<-estimateParameters(intamapExampleObject) ## make predictions at unobserved locations intamapExampleObject<-spatialCopula(intamapExampleObject) ## End(Not run)
interpolate is a function that interpolates spatial data
interpolate(observations, predictionLocations, outputWhat = list(mean = TRUE, variance = TRUE), obsChar = NA, methodName = "automatic", maximumTime = 30, optList = list(), cv = FALSE, ...) interpolateBlock(observations, predictionLocations, outputWhat, blockWhat = "none", obsChar = NA, methodName = "automatic", maximumTime = 30, optList = list(), ...)interpolate(observations, predictionLocations, outputWhat = list(mean = TRUE, variance = TRUE), obsChar = NA, methodName = "automatic", maximumTime = 30, optList = list(), cv = FALSE, ...) interpolateBlock(observations, predictionLocations, outputWhat, blockWhat = "none", obsChar = NA, methodName = "automatic", maximumTime = 30, optList = list(), ...)
observations |
observation data, object of class |
predictionLocations |
prediction locations, object of class |
outputWhat |
list with names what kind of output is expected, e.g.
|
blockWhat |
List defining particular output for block predictions. See |
obsChar |
list with observation characteristics, used by some interpolation methods |
methodName |
name of interpolation method to be used, see spatialPredict for more details, or |
maximumTime |
the maximum time available for interpolation, will be compared to the result of |
optList |
list; further options, mainly passed to
|
cv |
If cross-validation should be done |
... |
other arguments to be passed to other functions |
The functions interpolate and interpolateBlock are particularly implemented for being called by
a Web Processing Server (WPS), but they can also be used interactively. The only necessary
arguments are observations and predictionLocations. It is also recommended
to set outputWhat, and blockWhat if necessary. If outputWhat
contains nsim, the return table will also contain a number of realisations,
for methods able to return simulations.
interpolate can use different interpolation methods for the result. The function
will internally call the following functions which can be method specific.
An indication of available methods can be given by methods(estimateParameters)
or methods(spatialPredict).
The method can be set through the argument methodName, or through the
built-in automatic selection method. There are different criteria that helps
in selecting the right method for a particular data set. There are four
methods that are available for the automatic choice:
automap, psgp (from the separate package psgp)
copula and transgaussian are the possibilities.
First of all, if observation errors are present, the psgp method is preferred.
If not, it is checked whether the data appear to deviate significantly from normality.
This is assumed to be the case if any of the tests below are TRUE:
test[1] = length(boxplot.stats(dataObs)$out)/length(dataObs) > 0.1
test[2] = fivenum(dataObs)[3] - fivenum(dataObs)[2] < IQR(dataObs)/3
test[3] = fivenum(dataObs)[4] - fivenum(dataObs)[3] < IQR(dataObs)/3
g = boxcox(dataObs ~ 1,lambda=seq(-2.5,2.5,len=101),plotit=FALSE)$y
test[4] = g[71] < sort(g)[91]
where fivenum defines the Tukey five number statistic and
IQR finds the interquartile range of the data. If the minimum of
dataObs is <= 0, min(dataObs) + sdev(dataObs) is added to all values.
At last, the function calls predictTime for an estimate of the
prediction time. If any of the tests above were true and the estimated prediction time
for copula prediction is below maximumTime, the copula
method is chosen. If any of the
tests were TRUE and the estimated prediction time is too long, transGaussian
kriging is chosen, as long as all values are above zero. If any of
the tests are true for a set of observations with negative or zero-values,
automap is chosen, but a warning is issued.
The element methodParameters in the object being returned is a string that makes it possible
to regenerate the variogram model or the copula parameters in createIntamapObject.
This is particularly useful when the function is called through a WPS, when
the element with the estimated parameters cannot be preserved in a state
that makes it possible to use them for a later call to interpolate.
The possibility
for doing parallel processing 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 downloaded, additionally the parameter nclus must be set to
a value larger than 1. Parallel computation is not taken into account when
estimating the prediction times.
An intamap object, which is a list with elements, see intamap-package.
The exact number and names of these elements might vary due to different methods applied,
but the list below shows the most typical:
observations |
the observations, as a |
predictionLocations |
the prediction locations, as a |
formulaString |
the relationship between independent and dependent variables,
|
outputWhat |
a list of the prediction types to return |
anisPar |
the estimated anisotropy parameters |
variogramModel |
the estimated parameter for the method, can also be e.g. |
methodParameters |
a string, that when parsed, can be used to regenerate the variogram model or copula parameters. Useful for repeated calls to interpolate when it is not necessary to reestimate the parameters. |
predictions |
a |
outputTable |
a matrix, organized in a convenient way for the calling WPS;
first row: x-coordinates, second row: y-coordinates; further rows:
output elements as specified by |
processDescription |
some textual descriptions of the interpolation process, including warnings |
Edzer Pebesma
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.
createIntamapObject, estimateParameters,
spatialPredict, intamap-package
data(meuse) coordinates(meuse) = ~x+y meuse$value = meuse$zinc data(meuse.grid) gridded(meuse.grid) = ~x+y x = interpolate(meuse, meuse.grid, list(mean=TRUE, variance=TRUE)) summary(t(x$outputTable))data(meuse) coordinates(meuse) = ~x+y meuse$value = meuse$zinc data(meuse.grid) gridded(meuse.grid) = ~x+y x = interpolate(meuse, meuse.grid, list(mean=TRUE, variance=TRUE)) summary(t(x$outputTable))
function that generates a parsable string of identified method parameters for an intamap interpolation object
## Default S3 method: methodParameters(object) ## S3 method for class 'copula' methodParameters(object) ## S3 method for class 'idw' methodParameters(object)## Default S3 method: methodParameters(object) ## S3 method for class 'copula' methodParameters(object) ## S3 method for class 'idw' methodParameters(object)
object |
a list object. Most arguments necessary for interpolation
are passed through this object. See |
The function creates a text-string that makes it possible to add the
the method parameters (anisotropy and idw-parameter, variogram model or
copula parameters) to the object in a later call to
createIntamapObject or interpolate
without having to re-estimate the parameters.
This function is particularly useful when interpolate is
called from a Web Processing Service, and the user wants to reuse the
method parameters. The function is mainly assumed to be called from
within interpolate.
The default method assumes a variogram model of gstat type, e.g. a variogram
similar to what can be created with a call to vgm.
Also psgp uses this variogram model.
A string that, when parsed, will recreate the methodParameters
Jon Olav Skoien
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.
sessionInfo() data(meuse) coordinates(meuse) = ~x+y meuse$value = log(meuse$zinc) # set up intamap object: krigingObject = createIntamapObject( observations = meuse, formulaString = as.formula('value~1'),class = "automap") # do estimation steps: krigingObject = estimateParameters(krigingObject) krigingObject = methodParameters(krigingObject) # Create a new object krigingObject2 = createIntamapObject(observations = meuse, formulaString = as.formula('value~1'), params = list(methodParameters = krigingObject$methodParameters)) krigingObject$variogramModel krigingObject2$variogramModelsessionInfo() data(meuse) coordinates(meuse) = ~x+y meuse$value = log(meuse$zinc) # set up intamap object: krigingObject = createIntamapObject( observations = meuse, formulaString = as.formula('value~1'),class = "automap") # do estimation steps: krigingObject = estimateParameters(krigingObject) krigingObject = methodParameters(krigingObject) # Create a new object krigingObject2 = createIntamapObject(observations = meuse, formulaString = as.formula('value~1'), params = list(methodParameters = krigingObject$methodParameters)) krigingObject$variogramModel krigingObject2$variogramModel
Plotting function for intamap-objects of the type described in
intamap-package
plotIntamap(object, zcol = "all", sp.layout = NULL, plotMat = c(2,2), ...) ## S3 method for class 'copula' plot(x, ...) ## S3 method for class 'idw' plot(x, ...) ## S3 method for class 'automap' plot(x, ...) ## S3 method for class 'linearVariogram' plot(x, ...) ## S3 method for class 'transGaussian' plot(x, ...) ## S3 method for class 'yamamoto' plot(x, ...)plotIntamap(object, zcol = "all", sp.layout = NULL, plotMat = c(2,2), ...) ## S3 method for class 'copula' plot(x, ...) ## S3 method for class 'idw' plot(x, ...) ## S3 method for class 'automap' plot(x, ...) ## S3 method for class 'linearVariogram' plot(x, ...) ## S3 method for class 'transGaussian' plot(x, ...) ## S3 method for class 'yamamoto' plot(x, ...)
object |
a list object. Most arguments necessary for interpolation
are passed through this object. See |
x |
intamap object, when plot is called directly |
zcol |
a list of column names to be plotted; if equal to all, the column names
will correspond to all possible column names from |
sp.layout |
an object that can contain lines, points and polygons that function as extra layout;
see |
plotMat |
an array of lengt two with the number of rows and columns of plots per page |
... |
other parameters that can be passed to other plot functions (e.g.
|
All the plot methods above are simple wrapper functions around the plotIntamap function.
A plot of some of the elements of object. This will typically be the
sample variogram and the fitted variogram model (if a method based on variograms
has been used) and all the predictions.
Jon Olav Skoien
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.
data(meuse) meuse$value = log(meuse$zinc) data(meuse.grid) coordinates(meuse) = ~x+y coordinates(meuse.grid) = ~x+y object = interpolate(meuse, meuse.grid, outputWhat = list(mean = TRUE, variance = TRUE, excprob = 7, excprob = 8, quantile = 0.9, quantile = 0.95), methodName = "automap") plot(object)data(meuse) meuse$value = log(meuse$zinc) data(meuse.grid) coordinates(meuse) = ~x+y coordinates(meuse.grid) = ~x+y object = interpolate(meuse, meuse.grid, outputWhat = list(mean = TRUE, variance = TRUE, excprob = 7, excprob = 8, quantile = 0.9, quantile = 0.95), methodName = "automap") plot(object)
post-processing of data for the intamap-package. The function will typically call
other functions for adding back biases, aggregation etc.
## Default S3 method: postProcess(object, ...) ## S3 method for class 'idw' postProcess(object, ...)## Default S3 method: postProcess(object, ...) ## S3 method for class 'idw' postProcess(object, ...)
object |
a list object. Most arguments necessary for interpolation
are passed through this object. See |
... |
other parameters that can be passed to functions called from |
The function postProcess includes code for postprocessing an object after interpolation.
The function can easily be replaced by more specific methods relevant for a
certain data set, doing more data specific things in addition to what is done in the default method.
An object of same type as above, but with new elements. Most important from the default
method is the outputTable, a matrix, organized in a convenient way for the calling WPS;
first row: x-coordinates, second row: y-coordinates; further rows:
output elements as specified by outputWhat (see createIntamapObject
Edzer J. Pebesma
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.
Functions that gives a time estimate for an intamap function given the number of observations and predictionLocations
predictTime(nObs, nPred, class, formulaString, calibration=FALSE, outputWhat, FUN = "spatialPredict",...)predictTime(nObs, nPred, class, formulaString, calibration=FALSE, outputWhat, FUN = "spatialPredict",...)
nObs |
An integer or an array of integers containing the number of observations. |
nPred |
An integer or an array of integers containing the number of predictions. |
class |
class of intamapObject, which interpolation method to be used |
formulaString |
the formula of the request, mainly to see if the request has independent variables |
calibration |
enables or disables time calibration - not properly implemented yet |
outputWhat |
List defining the requested type of output, see
|
FUN |
A string with the intamap functions name, now obsolete |
... |
other arguments needed to define the intamap object. |
The function is based on timeModels being available in the workspace.
This is a loess-model that has been fitted to different calls to a range of
interpolation requests with synthetically generated data in
generateTimeModels.
An integer or an array of integers with the predicted times.
RUN FIRST generateTimeModels() or estimateTimeModel() in order to save time Models to workspace. It might take some time!
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.
pre-processing of the data for the intamap-package package.
## Default S3 method: preProcess(object, ...) ## S3 method for class 'idw' preProcess(object, ...)## Default S3 method: preProcess(object, ...) ## S3 method for class 'idw' preProcess(object, ...)
object |
see |
... |
Additional parameters |
The function preProcess includes code for preprocessing an object before interpolation.
The function can easily be replaced by more specific methods relevant for a
certain data set. Functions can be called from preProcess according to the settings in parameters
in the object, set by the function getIntamapParams.
The input object is returned, after its components have been pre-processed.
Jon Olav Skoien
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.
This function applies an isotropic transformation of
the coordinates specified in object.
rotateAnisotropicData(object,anisPar)rotateAnisotropicData(object,anisPar)
object |
(i) An Intamap type object (see |
anisPar |
An array containing the anisotropy parameters (anisotropy ratio and axes orientation)
(see |
This function performs a rotation and rescaling of the
coordinate axes in order to obtain a new coordinate system, in which
the observations become statistically isotropic. This assumes that
the estimates of the anisotropy ratio and the orientation angle
provided in anisPar are accurate.
(i) A modified object with transformed coordinates if
rotateAnisotropicData is called with an Intamap object as input (see
intamap-package) or (ii) the transformed coordinates if a
SpatialPointsDataFrame is used as input or (iii)
the transformed coordinates if a SpatialPoints
object is the input.
Hristopulos Dionisis, Spiliopoulos Giannis
[1] 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.
[2] A. Chorti and D. T. Hristopulos (2008). Non-parametric Identification of Anisotropic (Elliptic) Correlations in Spatially Distributed Data Sets, IEEE Transactions on Signal Processing, 56(10), 4738-4751 (2008).
estimateAnisotropy
library(gstat) data(sic2004) coordinates(sic.val)=~x+y sic.val$value=sic.val$dayx anisPar <- estimateAnisotropy(sic.val) print(anisPar) rotatedObs <- rotateAnisotropicData(sic.val,anisPar) newAnisPar <- estimateAnisotropy(rotatedObs) print(newAnisPar)library(gstat) data(sic2004) coordinates(sic.val)=~x+y sic.val$value=sic.val$dayx anisPar <- estimateAnisotropy(sic.val) print(anisPar) rotatedObs <- rotateAnisotropicData(sic.val,anisPar) newAnisPar <- estimateAnisotropy(rotatedObs) print(newAnisPar)
spatialPredict is a generic method for spatial predictions
within the intamap-package.
A series of methods have been implemented,
partly based on other R-packages (as krige),
other methods have been developed particularly for the INTAMAP project. The object
has to include a range of variables, further described in
intamap-package. The prediction method is
chosen based on the class of the object.
## S3 method for class 'automap' spatialPredict(object, nsim = 0, ...) ## S3 method for class 'copula' spatialPredict(object, ...) ## Default S3 method: spatialPredict(object, ...) ## S3 method for class 'idw' spatialPredict(object, ...) ## S3 method for class 'linearVariogram' spatialPredict(object, nsim = 0, ...) ## S3 method for class 'transGaussian' spatialPredict(object, nsim = 0, ...) ## S3 method for class 'yamamoto' spatialPredict(object, nsim = 0, ...)## S3 method for class 'automap' spatialPredict(object, nsim = 0, ...) ## S3 method for class 'copula' spatialPredict(object, ...) ## Default S3 method: spatialPredict(object, ...) ## S3 method for class 'idw' spatialPredict(object, ...) ## S3 method for class 'linearVariogram' spatialPredict(object, nsim = 0, ...) ## S3 method for class 'transGaussian' spatialPredict(object, nsim = 0, ...) ## S3 method for class 'yamamoto' spatialPredict(object, nsim = 0, ...)
object |
a list object. Most arguments necessary for interpolation
are passed through this object. See |
nsim |
number of simulations to return, for methods able to return simulations |
... |
other arguments that will be passed to the requested interpolation method. See the individual interpolation methods for more information. |
The function spatialPredict is a wrapper around different
spatial interpolation methods found within the intamap-package or within other packages
in R. It is for most of the
methods necessary to have parameters of the correlation structure
included in object to be able to carry out the spatial prediction.
Below are some details
about particular interpolation methods
defaulta default method is not really implemented, this function is only created to give a sensible error message if the function is called with an object for which no method exist
automapIf the object already has an element variogramModel with
variogram parameters,
krige is called. If the this is not a part of the object,
estimateParameters is called to create this element.
copulaspatial prediction using bayesCopula
idwapplies inverse distance modelling with the idp-power found by estimateParameters.idw
linearVariogramthis function estimates the process using an unfitted linear variogram; although variance is returned it can not be relied upon
transGaussianspatial prediction using krigeTg
yamamotospatial prediction using yamamotoKrige
It is also possible to add to the above methods with functionality from other packages, if wanted. You can also check which methods are available from other packages by calling
>methods(spatialPredict)
a list object similar to object, but extended with predictions at
a the set of locations defined object.
Jon Olav Skoien
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.
gstat,autoKrige,
createIntamapObject, estimateParameters,
intamap-package
# This example skips some steps that might be necessary for more complicated # tasks, such as estimateParameters and pre- and postProcessing of the data 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") # set up intamap object: obj = createIntamapObject( observations = meuse, predictionLocations = meuse.grid, targetCRS = "+init=epsg:3035", params = getIntamapParams(), class = "linearVariogram" ) # do interpolation step: obj = spatialPredict(obj) # spatialPredict.linearVariogram# This example skips some steps that might be necessary for more complicated # tasks, such as estimateParameters and pre- and postProcessing of the data 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") # set up intamap object: obj = createIntamapObject( observations = meuse, predictionLocations = meuse.grid, targetCRS = "+init=epsg:3035", params = getIntamapParams(), class = "linearVariogram" ) # do interpolation step: obj = spatialPredict(obj) # spatialPredict.linearVariogram
summary function for intamap-objects of the type described in
intamap-package
summaryIntamap(object, ...) ## S3 method for class 'copula' summary(object, ...) ## S3 method for class 'idw' summary(object, ...) ## S3 method for class 'automap' summary(object, ...) ## S3 method for class 'linearVariogram' summary(object, ...) ## S3 method for class 'transGaussian' summary(object, ...) ## S3 method for class 'yamamoto' summary(object, ...)summaryIntamap(object, ...) ## S3 method for class 'copula' summary(object, ...) ## S3 method for class 'idw' summary(object, ...) ## S3 method for class 'automap' summary(object, ...) ## S3 method for class 'linearVariogram' summary(object, ...) ## S3 method for class 'transGaussian' summary(object, ...) ## S3 method for class 'yamamoto' summary(object, ...)
object |
a list object. Most arguments necessary for interpolation
are passed through this object. See |
... |
parameters to be passed to the default summary function for each element |
A summary of some of the elements of object.
Jon Olav Skoien
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.
This is a standard model for estimating the prediction time within the
intamap-package. It was created by the function
generateTimeModels, on a 64 bits Linux server running
R version 2.9.0 and intamap-package version 1.1.15.
The prediction time will depend on the speed of the local machine, on the version
of R and intamap-package, and on the installed libraries.
It is therefore strongly recommended to run generateTimeModels on the local
machine and store the result in the workspace if the predicted interpolation
time is of real interest. timeModels in the workspace will be chosen if available.
It is not necessary to load the data set, this happens automatically in predictTime
if timeModels if the object is not already existing in the workspace.
data(timeModels)data(timeModels)
Jon Olav Skoien
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.
unbiasedKrige is a function for modifying a kriging prediction
to a prediction that can be assumed to be unbiased for a certain threshold.
unbiasedKrige(object, formulaString, observations, predictionLocations, model, outputWhat, yamamoto, iwqmaxit = 500, iwqCpAddLim = 0.0001, debug.level, ...)unbiasedKrige(object, formulaString, observations, predictionLocations, model, outputWhat, yamamoto, iwqmaxit = 500, iwqCpAddLim = 0.0001, debug.level, ...)
object |
either an object of the intamap type (see |
formulaString |
formula that defines the dependent variable as a linear model of independent variables; suppose the dependent variable has name z, for ordinary and simple kriging use the formula z~1; for universal kriging, suppose z is linearly dependent on x and y, use the formula z~x+y |
observations |
a |
predictionLocations |
the predictionLocations, only necessary if the method is "IWQSEL" and formulaString contains independent variables. Should preferentally be a grid if the method is "IWQSEL" |
model |
variogram model of dependent variable (or its residuals), defined
by a call to |
outputWhat |
Argument with type of unbiasedness method ("MOK" or "IWQSEL") and the thresholds. |
yamamoto |
logical describing if the yamamoto approach )is to be used in simulations.
Defaults to yamamoto = FALSE when object is a |
iwqmaxit |
maximum number of iterations in iwqsel |
iwqCpAddLim |
convergence criteria in iwqsel |
debug.level |
debug level, passed to subfunctions |
... |
other arguments that will be passed to subfunctions. These include
|
It is a fact that predictions from kriging tend to be biased towards the mean of
the process. The function unbiasedKrige is a function that adds one or more predictions
to the original output, which are assumed to be unbiased relative to a certain
threshold. The two methods supported are the IWQSEL-method (Craigmile, 2006) and
MOK (Skoien et al, 2008).
an object of type intamap, as described in intamap-package, or a
Spatial*DataFrame with one or more new prediction columns, representing different
methods and thresholds.
Jon Olav Skoien
Craigmile, P. F., N. Cressie, T. J. Santner, and Y. Rao. 2006. A loss function approach to identifying environmental exceedances. Extremes, 8, 143-159.
Skoien, J. O., G. B. M. Heuvelink, and E. J. Pebesma. 2008. Unbiased block predictions and exceedance probabilities for environmental thresholds. In: J. Ortiz C. and X. Emery (eds). Proceedings of the eight international geostatistics congress. Gecamin, Santiago, Chile, pp. 831-840.
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.
library(automap) library(gstat) data(meuse) data(meuse.grid) coordinates(meuse) = ~x+y gridded(meuse.grid) = ~x+y predictionLocations = meuse.grid[sample(1:length(meuse.grid),50),] vmod = autofitVariogram(log(zinc)~1,meuse)$var_model prediction = krige(log(zinc)~1,meuse,predictionLocations,vmod) summary(prediction) prediction <- unbiasedKrige(prediction,log(zinc)~1, meuse, model = vmod, outputWhat = list(MOK = 6.0, MOK = 7.0, IWQSEL=7.0), iwqmaxit = 100, iwqCpAddLim = 0.01) summary(prediction)library(automap) library(gstat) data(meuse) data(meuse.grid) coordinates(meuse) = ~x+y gridded(meuse.grid) = ~x+y predictionLocations = meuse.grid[sample(1:length(meuse.grid),50),] vmod = autofitVariogram(log(zinc)~1,meuse)$var_model prediction = krige(log(zinc)~1,meuse,predictionLocations,vmod) summary(prediction) prediction <- unbiasedKrige(prediction,log(zinc)~1, meuse, model = vmod, outputWhat = list(MOK = 6.0, MOK = 7.0, IWQSEL=7.0), iwqmaxit = 100, iwqCpAddLim = 0.01) summary(prediction)
ordinary kriging and simulation with an alternative kriging variance
yamamotoKrige(formula, Obs, newPoints, model, nsim = 0, nmax = 20, maxdist = Inf)yamamotoKrige(formula, Obs, newPoints, model, nsim = 0, nmax = 20, maxdist = Inf)
formula |
formula that defines the dependent variable as a linear model of
independent variables; suppose the dependent variable has name
|
Obs |
|
newPoints |
|
model |
variogram model - of the type that can be found by a call to
|
nsim |
integer; if set to a non-zero value, conditional simulation is used instead of kriging interpolation. For this, sequential Gaussian simulation is used, following a single random path through the data. |
nmax |
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, all observations are used. |
maxdist |
maximum number of neighbours to use in local kriging, defaults to Inf |
The term yamamotoKrige comes from the paper of Yamamoto (2000) where
he suggests using local variance around the kriging estimate (weighted with
the kriging weights) as an alternative kriging variance. This as a
solution to more reliable estimates of the kriging variance also when the
stationarity assumption has been violated. The method was applied by
Skoien et al. (2008), who showed that it can have advantages for cases where the
stationarity assumption behind kriging is violated.
If the number of observations is high, it is recommended have nmax lower.
This is partly because the method relies on positive kriging weights. The method
to do this adds the norm of the largest negative weight to all weights, and rescales.
This tends to smooth the weights, giving a prediction closer to the average if a
too large number of observation locations is used.
Either a Spatial*DataFrame with predictions and prediction variance,
in the columns var1.pred and var1.var, together with the
classical ordinary kriging variance in var1.ok, or simulations with
column names simx where x is the number of the simulation.
Jon Olav Skoien
Skoien, J. O., G. B. M. Heuvelink, and E. J. Pebesma. 2008. Unbiased block predictions and exceedance probabilities for environmental thresholds. In: J. Ortiz C. and X. Emery (eds). Proceedings of the eight international geostatistics congress. Santiago, Chile: Gecamin, pp. 831-840.
Yamamoto, J. K. 2000. An alternative measure of the reliability of ordinary kriging estimates. Mathematical Geology, 32 (4), 489-509.
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.
library(gstat) library(automap) data(sic2004) coordinates(sic.val) = ~x+y coordinates(sic.test) = ~x+y variogramModel = autofitVariogram(joker~1,sic.val)$var_model newData = yamamotoKrige(joker~1,sic.val,sic.test,variogramModel,nmax = 20) summary(newData) plot(sqrt(var1.ok)~var1.pred,newData) # Kriging variance the same in regions with extreme values plot(sqrt(var1.var)~var1.pred,newData) # Kriging standard deviation higher for high predictions (close to extreme values)library(gstat) library(automap) data(sic2004) coordinates(sic.val) = ~x+y coordinates(sic.test) = ~x+y variogramModel = autofitVariogram(joker~1,sic.val)$var_model newData = yamamotoKrige(joker~1,sic.val,sic.test,variogramModel,nmax = 20) summary(newData) plot(sqrt(var1.ok)~var1.pred,newData) # Kriging variance the same in regions with extreme values plot(sqrt(var1.var)~var1.pred,newData) # Kriging standard deviation higher for high predictions (close to extreme values)