| Title: | Interpolation of Data with Variable Spatial Support |
|---|---|
| Description: | Data with irregular spatial support, such as runoff related data or data from administrative units, can with 'rtop' be interpolated to locations without observations with the top-kriging method. A description of the package is given by Skøien et al (2014) <doi:10.1016/j.cageo.2014.02.009>. |
| Authors: | Jon Olav Skøien [aut, cre], Qingyun Duan [ctb] (For the original FORTRAN code of SCE-UA, translated and modified to R in this package.) |
| Maintainer: | Jon Olav Skøien <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.6-9 |
| Built: | 2024-10-27 03:04:11 UTC |
| Source: | https://github.com/cran/rtop |
This package provides geostatistical methods for analysis and interpolation of data that has an irregular support, such as as runoff characteristics or population health data. The methods in this package are based on the top-kriging approach suggested in Skoien et al (2006), with some extensions from Gottschalk (1993). This package can be used as an add-on package for the automatic interpolation package developed within the intamap project (www.intamap.org).
The work flow within the package suggests that the user is interested in a prediction of a process at a series of locations where observations have not been made. The example below shows a regionalization of mean annual runoff in Austria.
Although it is possible to perform each step with all necessary arguments,
the easiest interface to
the method is to store all variables (such as observations, prediction locations
and parameters) in an rtop-object, which is created by a call to
createRtopObject. The element params below consists of
changes to the default parameters. A further description can be found
in getRtopParams. The changes below means that
the functions will use geostatistical distance instead of full regularization,
and that the variogram model will be fitted to the variogram cloud.
Most other functions in the rtop-package can take this object as an argument,
and will add the results as one or more new element(s) to this object.
The data in the example below are stored as shape-files
in the extdata-directory of the rtop-pacakge, use the directory of your own data instead.
The observations consist of mean summer runoff
from 138 catchments in Upper Austria. The predictionLocations are 863 catchments
in the same region. observations and predictionLocations can either be stored as
SpatialPolygonsDataFrame-objects, or as sf-polygons.
rpath = system.file("extdata",package="rtop")
library(sf)
observations = st_read(rpath, "observations")
predictionLocations = st_read(rpath,"predictionLocations")
# Create a column with the specific runoff:
observations$obs = observations$QSUMMER_OB/observations$AREASQKM
params = list(gDist = TRUE, cloud = TRUE)
rtopObj = createRtopObject(observations, predictionLocations,
params = params)
There are help-methods available in cases when data are not available as
shape-files, or when the observations are not part of the shape-files.
See readAreaInfo and readAreas.
A call to rtopVariogram adds the sample variogram to the object,
whereas rtopFitVariogram fits a variogram model. The last function
will call rtopVariogram if rtopObj does not contain a sample variogram.
rtopObj = rtopVariogram(rtopObj) rtopObj = rtopFitVariogram(rtopObj, maxn = 2000)
The function checkVario is useful to produce some
diagnostic plots for the sample variogram and the fitted variogram model.
checkVario(rtopObj)
The interpolation function (rtopKrige) solves the kriging system based on the
computed regularized semivariances. The covariance matrices are created in a
separate regularization function (varMat), and are stored in
the rtop-object for easier access if it is necessary to redo parts of the
analysis, as this is the computationally expensive part of the interpolation.
Cross-validation can be called with the argument cv=TRUE, either in
params or in the call to rtopKrige.
rtopObj = rtopKrige(rtopObj)
if (is(rtopObj$predictions, "Spatial")) {
spplot(rtopObj$predictions, col.regions = bpy.colors(), c("var1.pred"))
} else {
# the plotting order to get small polygons on top is not automatic with sf,
# but here is a method that works without modifying the predictions
library(dplyr)
# Arrange according to areas attribute in descending order
rtopObj$predictions |> arrange(desc(AREASQKM)) |>
# Make ggplot and set fill color to var1.pred
ggplot(aes(fill = var1.pred)) + geom_sf()
}
rtopObj = rtopKrige(rtopObj, cv = TRUE)
if (is(rtopObj$predictions, "Spatial")) {
spplot(rtopObj$predictions, col.regions = bpy.colors(), c("var1.pred","var1.var"))
} else {
# Here is an alternative method for plotting small polygons on top of the larger ones,
# modifying the predictions
rtopObj$predictions = rtopObj$predictions[order(rtopObj$predictions$AREASQ, decreasing = TRUE), ]
# It is also possible to change the order of the polygons
ggplot(rtopObj$predictions) + aes(fill = var1.pred) + geom_sf() + scale_fill_distiller(palette = "YlOrRd")
}
L. Gottschalk. Interpolation of runoff applying objective methods. Stochastic Hydrology and Hydraulics, 7:269-281, 1993.
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
The function will create diagnostic plots for analysis of the variograms fitted to sample variograms of data with support
## S3 method for class 'rtop' checkVario(object, acor = 1, log = "xy", cloud = FALSE, gDist = TRUE, acomp = NULL, curveSmooth = FALSE, params = list(), ...) ## S3 method for class 'rtopVariogramModel' checkVario(object, sampleVariogram = NULL, observations = NULL, areas = NULL, dists = NULL, acomp = NULL, params = list(), compVars = list(), acor = 1, log = "xy", legx = NULL, legy = NULL, plotNugg = TRUE, curveSmooth = FALSE, ...)## S3 method for class 'rtop' checkVario(object, acor = 1, log = "xy", cloud = FALSE, gDist = TRUE, acomp = NULL, curveSmooth = FALSE, params = list(), ...) ## S3 method for class 'rtopVariogramModel' checkVario(object, sampleVariogram = NULL, observations = NULL, areas = NULL, dists = NULL, acomp = NULL, params = list(), compVars = list(), acor = 1, log = "xy", legx = NULL, legy = NULL, plotNugg = TRUE, curveSmooth = FALSE, ...)
object |
either: object of class |
acor |
unit correction factor in the key, e.g. to see numbers more easily interpretable for large areas. As an example, ucor = 0.000001 when area is given in square meters and should rather be shown as square kilometers. Note that this parameter also changes the value of the nugget to the new unit. |
log |
text variable for log-plots, default to log-log |
cloud |
logical; whether to look at the cloud variogram instead of the binned variogram |
gDist |
logical; whether to use ghosh-distance for semivariogram regularization instead of full integration of the semivariogram |
sampleVariogram |
a sample variogram of the data |
observations |
a set of observations |
areas |
either an array of areas that should be used as examples, or
the number of areas per order of magnitude (similar to the parameter |
dists |
either an array of distances that should be used as examples, or
the number of distances per order of magnitude(similar to the parameter |
acomp |
either a matrix with the area bins that should be visualized, or a number
giving the number of pairs to show. If a sample variogram is given, the |
curveSmooth |
logical or numerical; describing whether the curves in the last plot should be smoothed or not. If numeric,
it gives the degrees of freedom (df) for the splines used for smoothing. See also |
params |
list of parameters to modify the default parameters of rtopObj or
the default parameters found from |
compVars |
a list of variograms of |
legx |
x-coordinate of the legend for fine-tuning of position, see x-argument of |
legy |
y-coordinate of the legend for fine-tuning of position, see y-argument of |
plotNugg |
logical; whether the nugget effect should be added to the plot or not |
... |
arguments to lower level functions |
The function gives diagnostic plots for the fitted variograms, where the regularized variograms are shown together with the sample variograms and possibly also user defined variograms. In addition, if an rtopObject is submitted, the function will also give plots of the relationship between variance and area size and a scatter plot of the fit of the observed and regularized variogram values. The sizes of the dots are relative to the number of pairs in each group.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
library(gstat) rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM params = list(cloud = TRUE, gDist = TRUE) rtopObj = createRtopObject(observations, predictionLocations, params = params) # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj, maxn = 2000) checkVario(rtopObj, compVars = list(first = vgm(5e-6, "Sph", 30000,5e-8), second = vgm(2e-6, "Sph", 30000,5e-8))) rtopObj = checkVario(rtopObj, acor = 0.000001, acomp = data.frame(acl1 = c(2,2,2,2,3,3,3,4,4), acl2 = c(2,3,4,5,3,4,5,4,5))) rtopObj = checkVario(rtopObj, cloud = TRUE, identify = TRUE, acor = 0.000001)library(gstat) rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM params = list(cloud = TRUE, gDist = TRUE) rtopObj = createRtopObject(observations, predictionLocations, params = params) # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj, maxn = 2000) checkVario(rtopObj, compVars = list(first = vgm(5e-6, "Sph", 30000,5e-8), second = vgm(2e-6, "Sph", 30000,5e-8))) rtopObj = checkVario(rtopObj, acor = 0.000001, acomp = data.frame(acl1 = c(2,2,2,2,3,3,3,4,4), acl2 = c(2,3,4,5,3,4,5,4,5))) rtopObj = checkVario(rtopObj, cloud = TRUE, identify = TRUE, acor = 0.000001)
This is a help function for creating an object (see rtop-package
to be used for interpolation within the rtop package
createRtopObject(observations, predictionLocations, formulaString, params=list(), overlapObs, overlapPredObs, ...)createRtopObject(observations, predictionLocations, formulaString, params=list(), overlapObs, overlapPredObs, ...)
observations |
|
predictionLocations |
a |
formulaString |
formula that defines the dependent variable as a linear model
of independent variables; suppose the dependent variable has name |
params |
parameters to modify the default parameters of the rtop-package,
set internally in this function by a call to |
overlapObs |
matrix with observations that overlap each other |
overlapPredObs |
matrix with |
... |
Extra parameters to |
An object of class rtop with observations, prediction locations,
parameters and possible other elements useful for interpolation in the rtop-package.
Most other externally visible functions in the package will be able to
work with this object, and add the results as a new element.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Setting some parameters params = list(gDist = TRUE, cloud = FALSE) # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Build an object rtopObj = createRtopObject(observations, predictionLocations, params = params)rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Setting some parameters params = list(gDist = TRUE, cloud = FALSE) # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Build an object rtopObj = createRtopObject(observations, predictionLocations, params = params)
Download additional example data from Vienna University of Technology
downloadRtopExampleData(folder = system.file("extdata", package="rtop"))downloadRtopExampleData(folder = system.file("extdata", package="rtop"))
folder |
the folder to which the downloaded data set will be copied |
The function will have as a side effect that additional example data is
downloaded from Vienna University of Techology. This will for the default
case replace the existing example data-set in the rtop package. Alternatively
the user can specify a separate directory for the data set.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
## Not run: downloadRtopExampleData() rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath,"observations") ## End(Not run)## Not run: downloadRtopExampleData() rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath,"observations") ## End(Not run)
Calculate geostatistical distances (Ghosh-distances) between areas
## S3 method for class 'rtop' gDist(object, params = list(), ...) ## S3 method for class 'SpatialPolygonsDataFrame' gDist(object, object2 = NULL, ...) ## S3 method for class 'SpatialPolygons' gDist(object, object2 = NULL, ...) ## S3 method for class 'list' gDist(object, object2 = NULL, diag = FALSE, debug.level = 0, ...)## S3 method for class 'rtop' gDist(object, params = list(), ...) ## S3 method for class 'SpatialPolygonsDataFrame' gDist(object, object2 = NULL, ...) ## S3 method for class 'SpatialPolygons' gDist(object, object2 = NULL, ...) ## S3 method for class 'list' gDist(object, object2 = NULL, diag = FALSE, debug.level = 0, ...)
object |
object of class |
params |
a set of parameters, used to modify the default parameters for
the |
object2 |
an object of same type as |
diag |
logical; if TRUE only calculate the geostatistical distances between each element and itself, only when the objects are lists of discretized areas and object2 = object or object2 = NULL |
debug.level |
debug.level = 0 will suppress output from the call to varMat, done for calculation of the geostatistical distances |
... |
other parameters, for |
If called with one list of discretized elements, a matrix with the geostatistical distances
between the elements within the list. If called with two lists of discretized elements,
a matrix with the geostatistical distances between the elements in the two lists.
If called with diag = TRUE, the function returns an array of the geostatistical
distance within each of the elements in the list.
If called with one SpatialPolygons or SpatialPolygonsDataFrame
or the function returns a list with one matrix with geostatistical distances between
the elements of the object. If called with two objects, the list will also containt
a matrix of the geostatistical distances between the elements of the two objects, and an array
of the geostatistical distances within the elements of the second object.
If called with an rtop-object, the function will return the object, amended with the list above.
The geostatistical distance can be seen as the average distance between points in two elements, or the average distance within points in a single element. The distance measure is also sometimes referred to as Ghosh-distance, from Ghosh (1951) who found analytical expressions for these distances between blocks with regular geometry.
The use of geostatistical distances within rtop is based on an idea
from Gottschalk (1993), who suggested
to replace the traditional regularization of variograms within block-kriging
(as done in the original top-kriging application of Skoien et al (2006))
with covariances of the geostatistical distance. The covariance between two
areas can then be found as C(a1,a2) = cov(gd) where gd is the geostatistical
distance between the two areas a1 and a2, instead of an integration
of the covariance function between the two areas.
rtop is based on semivariograms
instead of covariances, and the semivariogram value between the two areas
can be found as gamma(a1,a2) = g(gd) - 0.5 (g(gd1) + g(gd2)) where
g is a semivariogram valid for point support, gd1) and gd2
are the geostatistical distances within each of the two areas.
Jon Olav Skoien
Ghosh, B. 1951. Random distances within a rectangle and between two rectangles. Bull. Calcutta Math. Soc., 43, 17-24.
Gottschalk, L. 1993. Correlation and covariance of runoff. Stochastic Hydrology and Hydraulics, 7, 85-101.
Skoien, J. O., R. Merz, and G. Bloschl. 2006. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10, 277-287.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") gDist = gDist(observations)rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") gDist = gDist(observations)
This function sets a range of the parameters for the intamap package,
to be included in the object described in rtop-package
getRtopParams(params,newPar, observations, formulaString, ...)getRtopParams(params,newPar, observations, formulaString, ...)
params |
An existing set of parameters for the interpolation process,
of class |
newPar |
A |
observations |
|
formulaString |
formula that defines the dependent variable as a linear model
of independent variables, see e.g. |
... |
Individual parameters for updating
|
A list of the parameters with class rtopParams to be included in the
object described in rtop-package
This function will mainly be called by createRtopObject, 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
rtopParams, the function will assume that the argument(s) are
modifications to the default set of parameters. The function can also be called
by other functions in the rtop-package if the users chooses not to work with
an object of class rtop.
If the function is called with two lists of parameters (but the first one is
not of class rtopParams) 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.
Parallel processing has been included for some of the functions. The default is no parallel procesing, and the package also attempts to decide whether it is sensible to start a set of clusters and distribute jobs to them based on the size of the job. The default limit might not be the best for every system.
Jon Olav Skoien
Cressie, N. 1985. Fitting variogram models by weighted least squares. Mathematical Geology, 17 (5), 563-586
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
createRtopObject and rtop-package
# Create a new set of intamapParameters, with default parameters: params = getRtopParams() # Make modifications to the default list of parameters params = getRtopParams(newPar = list(gDist = TRUE, nugget = FALSE)) # Make modifications to an existing list of parameters params = getRtopParams(params = params, newPar = list(gDist = TRUE, nugget = FALSE))# Create a new set of intamapParameters, with default parameters: params = getRtopParams() # Make modifications to the default list of parameters params = getRtopParams(newPar = list(gDist = TRUE, nugget = FALSE)) # Make modifications to an existing list of parameters params = getRtopParams(params = params, newPar = list(gDist = TRUE, nugget = FALSE))
Plot a sample variogram cloud, possibly with identification of individual point pairs
## S3 method for class 'rtopVariogramCloud' plot(x, ...)## S3 method for class 'rtopVariogramCloud' plot(x, ...)
x |
object of class |
... |
parameters that are passed through to
|
This function is mainly a wrapper around plot.variogramCloud,
necessary because of different column names and different class names. The
description of arguments and value can therefore be found in
the help page of plot.variogramCloud.
Jon Olav Skoien
rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Create the sample variogram rtopVario = rtopVariogram(observations, params = list(cloud = TRUE)) plot(rtopVario)rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Create the sample variogram rtopVario = rtopVariogram(observations, params = list(cloud = TRUE)) plot(rtopVario)
readAreaInfo will read a text file with observations and descriptions of data with a spatial support.
readAreaInfo(fname = "ainfo.txt", id = "id", iobs = "iobs", obs = "obs", unc = "unc", filenames= "filenames", sep = "\t", debug.level = 1, moreCols = list(NULL))readAreaInfo(fname = "ainfo.txt", id = "id", iobs = "iobs", obs = "obs", unc = "unc", filenames= "filenames", sep = "\t", debug.level = 1, moreCols = list(NULL))
fname |
name of file with areal information |
id |
name of column with observation id |
iobs |
name of column with number of observations |
obs |
name of column with observations |
unc |
name of column with possible uncertainty of observation |
filenames |
name of column with filenames of areas if different names than id should be used. |
sep |
separator in csv-file |
debug.level |
used for giving additional output |
moreCols |
name of other column names the user wants included in ainfo |
The function is of particular use when data are not available as
shape-files, or when the observations are not part of the shape-files.
This function is mainly for compatibility with the former FORTRAN-version.
The simplest way to read the data in that case is through st_read. See
also rtop-package.
SpatialPointDataFrame with information about observations and/or
predictionLocations.
Jon Olav Skoien
readAreas will read area-files, add observations and convert the result to SpatialPolygonsDataFrame
readAreas(object, adir=".",ftype = "xy",projection = NA, ...)readAreas(object, adir=".",ftype = "xy",projection = NA, ...)
object |
either name of file with areal information or |
adir |
directory where the files with areal information are to be found |
ftype |
type of file, the only type supported currently is "xy", referring to x- and y-coordinates of boundaries |
projection |
add projection to the object if input is boundary-files |
... |
further parameters to be passed to |
If object is a file name, readAreaInfo will be called.
If it is a SpatialPointsDataFrame with observations and/or
predictionLocations, the function will read areal data from files according
to the ID associated with each observation/predictionLocation.
The function is of particular use when data are not available as
shape-files, or when the observations are not part of the shape-files.
This function is mainly for compatibility with the former FORTRAN-version.
The simplest way to read the data in that case is through st_read. See
also rtop-package.
The function creates a SpatialPolygonsDataFrame of observations
and/or predictionLocations, depending on the information given in object.
Jon Olav Skoien
Convenience function for using parallel computation with rtop. The function is usually not called by the user.
rtopCluster(nclus, ..., action = "start", type, outfile = NULL)rtopCluster(nclus, ..., action = "start", type, outfile = NULL)
nclus |
The number of workers in the cluster |
... |
Arguments for |
action |
Defines the action of the function. There are three options:
|
type |
The type of cluster; see |
outfile |
File to direct the output, |
It is usually not necessary for the user to call this function for starting or accessing a cluster. This is done automatically
by the different rtop-functions when needed if the parameter nclus is larger than one (see getRtopParams). If the user actually starts the cluster by a call to this function, it will
also be necessary to set the nclus parameter to a value larger than one for the cluster to be used by different functions.
Restarting the cluster might be necessary if the cluster has a problem (e.g. does not return memory) or if the user wants to change to a different cluster type.
Stopping the cluster is useful when the user does not want to continue with parallel computation and wants to close down the workers.
If the function is called with action equal to "start" or "restart", the result is a cluster with nclus workers.
The cluster is also added to the global options with the name rtopCluster
(getOption("rtopCluster")).
If the function is called with action equal to "stop", the function stops the cluster, sets the rtopCluster of options to NULL and returns NULL to the user.
Jon Olav Skoien
rtopDisc will discretize an area for
regularization or calculation of Ghosh-distance
## S3 method for class 'rtop' rtopDisc(object, params = list(),...) ## S3 method for class 'SpatialPolygonsDataFrame' rtopDisc(object, params = list(), bb = bbox(object), ...) ## S3 method for class 'SpatialPolygons' rtopDisc(object, params = list(), bb = bbox(object), ...) ## S3 method for class 'rtopVariogram' rtopDisc(object, params = list(), ...)## S3 method for class 'rtop' rtopDisc(object, params = list(),...) ## S3 method for class 'SpatialPolygonsDataFrame' rtopDisc(object, params = list(), bb = bbox(object), ...) ## S3 method for class 'SpatialPolygons' rtopDisc(object, params = list(), bb = bbox(object), ...) ## S3 method for class 'rtopVariogram' rtopDisc(object, params = list(), ...)
object |
object of class |
bb |
boundary box, usually modified to be the common boundary box for two spatial object |
params |
possibility to pass parameters to modify the default parameters for
the
|
... |
Possibility to pass individual parameters |
There are different options for discretizing the objects. When the areas
from the bins are discretized, the options are random or regular sampling,
regular sampling is the default.
For the real areas, regular sampling appears to have computational advantages compared
with random sampling. In addition to the traditional regular sampling, rtop
also offers a third type of sampling which assures that the same discretization
points are used for overlapping areas.
Starting with a coarse grid covering the region of interest, this will for a certain support be refined till a requested minimum number of points from the grid is within the support. In this way, for areal supports, the number of points in the area with the largest number of points will be approximately four times the requested minimum number of points. This methods also assure that points used to discretize a large support will be reused when discretizing smaller supports within the large one, e.g. subcatchments within larger catchments.
The function returns a list of discretized areas, or if called with an rtop-object as argument, the object with lists of discretizations of the observations and prediction locations (if part of the object). If the function is called with an rtopVariogram (usually this is an internal call), the list contains discretized pairs of hypothetical objects from each bin of the semivariogram with a centre-to-centre distance equal to the average distance between the objects in a certain bin.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
rtopFitVariogram will fit a variogram model to the estimated binned variogram or cloud variogram of data with an areal support.
## S3 method for class 'rtop' rtopFitVariogram(object, params = list(), iprint = 0, ...) ## S3 method for class 'SpatialPolygonsDataFrame' rtopFitVariogram(object, params=list(), ...) ## S3 method for class 'SpatialPointsDataFrame' rtopFitVariogram(object, params=list(), ...) ## S3 method for class 'rtopVariogram' rtopFitVariogram(object, observations, dists = NULL, params=list(), mr = FALSE, aOver = NULL, iprint = 0, ...) ## S3 method for class 'rtopVariogramCloud' rtopFitVariogram(object, observations, dists = NULL, aOver = NULL, params=list(), mr = FALSE, iprint = 0, ...)## S3 method for class 'rtop' rtopFitVariogram(object, params = list(), iprint = 0, ...) ## S3 method for class 'SpatialPolygonsDataFrame' rtopFitVariogram(object, params=list(), ...) ## S3 method for class 'SpatialPointsDataFrame' rtopFitVariogram(object, params=list(), ...) ## S3 method for class 'rtopVariogram' rtopFitVariogram(object, observations, dists = NULL, params=list(), mr = FALSE, aOver = NULL, iprint = 0, ...) ## S3 method for class 'rtopVariogramCloud' rtopFitVariogram(object, observations, dists = NULL, aOver = NULL, params=list(), mr = FALSE, iprint = 0, ...)
object |
object of class The object can also be of class |
observations |
the observations, passed as a Spatial*DataFrame object, if
object is an |
params |
a set of parameters, used to modify the default parameters for
the |
dists |
either a matrix with geostatistical distances (created by a call
to the function |
mr |
logical; defining whether the function should return a list with discretized elements and geostatistical distances, even if it was not called with an rtop-object as argument. |
aOver |
a matrix with the overlapping areas of the observations, used for computation of the nugget effect. It will normally be recomputed by the function if it is NULL and necessary |
iprint |
print flag that is passed to |
... |
Other parameters to functions called from |
The function creates an object with the fitted
variogram Model (variogramModel)
and a data.frame (varFit) with the
differences between the sample semivariances and the regularized semivariances.
If mr = TRUE, the function also returns other objects (discretized elements
and geostatistical distances, if created) as a part of the returned object.
If the function is called with an rtop-object as argument, it will return an
rtop-object with variogramModel and varFit added to the object,
in addition to other objects created.
There are several options for fitting of the variogramModel, where the parameters
can be set in params, which is a list of parameters for modification
of the default parameters of the rtop-package given in a call to
getRtopParams. The first choice is between individual fitting and binned
fitting. This is based on the type of variogram submitted, individual fitting is done
if a cloud variogram (of class rtopVariogramCloud) is passed as argument,
and binned fitting if the submitted variogram is of class rtopVariogram.
If the function is called with an object of class rtop, having both
variogram and variogramCloud among its arguments, the variogram
model is fitted to the variogram which is consistent with the parameter cloud.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O. and G. Bloschl. Spatio-Temporal Top-Kriging of Runoff Time Series. Water Resources Research 43:W09419, 2007.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Setting some parameters params = list(gDist = TRUE, cloud = FALSE) # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Build an object rtopObj = createRtopObject(observations,predictionLocations, params = params) # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj) rtopObj$variogramModelrpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Setting some parameters params = list(gDist = TRUE, cloud = FALSE) # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Build an object rtopObj = createRtopObject(observations,predictionLocations, params = params) # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj) rtopObj$variogramModel
rtopKrige perform spatial interpolation or cross validation of data with areal support.
## S3 method for class 'rtop' rtopKrige(object, varMatUpdate = FALSE, params = list(), ...) ## S3 method for class 'SpatialPolygonsDataFrame' rtopKrige(object, predictionLocations = NULL, varMatObs, varMatPredObs, varMat, params = list(), formulaString, sel, ...) ## S3 method for class 'STSDF' rtopKrige(object, predictionLocations = NULL, varMatObs, varMatPredObs, varMat, params = list(), formulaString, sel, olags = NULL, plags = NULL, lagExact = TRUE, ...) ## Default S3 method: rtopKrige(object, predictionLocations = NULL, varMatObs, varMatPredObs, varMat, params = list(), formulaString, sel, wret = FALSE, ...)## S3 method for class 'rtop' rtopKrige(object, varMatUpdate = FALSE, params = list(), ...) ## S3 method for class 'SpatialPolygonsDataFrame' rtopKrige(object, predictionLocations = NULL, varMatObs, varMatPredObs, varMat, params = list(), formulaString, sel, ...) ## S3 method for class 'STSDF' rtopKrige(object, predictionLocations = NULL, varMatObs, varMatPredObs, varMat, params = list(), formulaString, sel, olags = NULL, plags = NULL, lagExact = TRUE, ...) ## Default S3 method: rtopKrige(object, predictionLocations = NULL, varMatObs, varMatPredObs, varMat, params = list(), formulaString, sel, wret = FALSE, ...)
object |
object of class |
varMatUpdate |
logical; if TRUE, also existing variance matrices will
be recomputed, if FALSE, only missing variance matrices will be computed,
see also |
predictionLocations |
|
varMatObs |
covariance matrix of observations, where diagonal must consist
of internal variance, typically generated from call
to |
varMatPredObs |
covariance matrix between observation locations and
prediction locations, typically generated from call
to |
varMat |
list covariance matrices including the two above |
params |
a set of parameters, used to modify the default parameters for
the |
formulaString |
formula that defines the dependent variable as a linear model
of independent variables, see e.g. |
sel |
array of prediction location numbers, if only a limited number of locations are to be interpolated/crossvalidated |
wret |
logical; if TRUE, return a matrix of weights instead of the predictions, useful for batch processing of time series, see also details |
olags |
A vector describing the relative lag which should be applied for the observation locations. See also details |
plags |
A vector describing the relative lag which should be applied for the predicitonLocations. See also details |
lagExact |
logical; whether differences in lagtime should be computed exactly or approximate |
... |
from |
This function is the interpolation routine of the rtop-package.
The simplest way of calling the function is with an rtop-object that
contains the fitted variogram model and all the other necessary data (see
createRtopObject or rtop-package).
The function will, if called with covariance matrices between observations
and between observations and prediction locations, use these for the interpolation.
If the function is called without these matrices, varMat will be
called to create them. These matrices can therefore be reused if necessary,
an advantage as it is computationally expensive to create them.
The interpolation that takes part within rtopKrige.default is based on
the semivariance matrices between observations and between observations and prediction
locations. It is therefore possible to use this function also to interpolate
data where the matrices have been created in other ways, e.g. based on distances
in physiographical space or distances along a stream.
The function returns the weights rather than the predictions if wret = TRUE.
This is useful for batch processing of time series, e.g. once the weights are
created, they can be used to compute the interpolated values for each time step.
rtop is able to take some advantage of multiple CPUs, which can be invoked with the
parameter nclus. When it gets a number larger than one, rtopKrige will start a cluster with nclus workers,
if the parallel-package has been installed.
The parameter singularSolve can be used when some areas are almost completely overlapping. In this case, the discretization of them might be equal, and the covariances to other areas will also be equal. The kriging matrix will in this case be singular. When singularSolve = TRUE, rtopKrige will remove one of the neighbours, and instead work with the mean of the two observations. An overview of removed neighbours can be seen in the resulting object, under the name removed.
Kriging of time series is possible when observations and predictionLocations
are spatiotemporal objects of type STSDF. The interpolation is
still spatial, in the sense that the regularisation of the variograms are just done
using the spatial extent of the observations, not a possible temporal extent, such as
done by Skoien and Bloschl (2007). However, it is possible to make predictions based on observations
from different time steps, through the use of the lag-vectors. These vectors describe a typical "delay"
for each observation and prediction location. This delay could for runoff related variables be similar
to travel time to each gauging location. For a certain prediction location, earlier time steps would be picked for neighbours with shorter travel time and later time steps for neighbours with slower travel times.
The lagExact parameter indicates whether to use a weighted average of two time steps, or just the time step which is closest to the difference in lag times.
The use of lag times should in theory increase the computation time, but might, due to different computation methods, even speed up the computation when the number of neighbours to be used (parameter nmax) is small compared to the number of observations. If computation is slow, it can be useful to test olags = rep(0, dim(observations)[1]) and similar for predictionLocations.
If called with SpatialPolygonsDataFrame, the function returns a SpatialPolygonsDataFrame with predictions, either at the
locations defined in predictionLocations, or as leave-one-out
cross-validation predicitons at the same locations as in object if
cv = TRUE
If called with an rtop-object, the function returns the same object with the predictions added to the object.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O. and G. Bloschl. Spatio-Temporal Top-Kriging of Runoff Time Series. Water Resources Research 43:W09419, 2007.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
# The following command will download the complete example data set # downloadRtopExampleData() # observations$obs = observations$QSUMMER_OB/observations$AREASQKM rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Setting some parameters; nclus > 1 will start a cluster with nclus # workers for parallel processing params = list(gDist = TRUE, cloud = FALSE, nclus = 1, rresol = 25) # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Build an object rtopObj = createRtopObject(observations, predictionLocations, params = params) # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj) # Predicting at prediction locations rtopObj = rtopKrige(rtopObj) # Cross-validation rtopObj = rtopKrige(rtopObj,cv=TRUE) cor(rtopObj$predictions$observed,rtopObj$predictions$var1.pred)# The following command will download the complete example data set # downloadRtopExampleData() # observations$obs = observations$QSUMMER_OB/observations$AREASQKM rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Setting some parameters; nclus > 1 will start a cluster with nclus # workers for parallel processing params = list(gDist = TRUE, cloud = FALSE, nclus = 1, rresol = 25) # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Build an object rtopObj = createRtopObject(observations, predictionLocations, params = params) # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj) # Predicting at prediction locations rtopObj = rtopKrige(rtopObj) # Cross-validation rtopObj = rtopKrige(rtopObj,cv=TRUE) cor(rtopObj$predictions$observed,rtopObj$predictions$var1.pred)
rtopSim will conditionally or unconditionally simulate data with areal support. This function should be seen as experimental, some issues are described below.
## S3 method for class 'rtop' rtopSim(object, varMatUpdate = FALSE, beta = NA, largeFirst = TRUE, replace = FALSE, params = list(), dump = NULL, debug.level, ...) ## Default S3 method: rtopSim(object = NULL, predictionLocations, varMatObs, varMatPredObs, varMatPred, variogramModel, ...)## S3 method for class 'rtop' rtopSim(object, varMatUpdate = FALSE, beta = NA, largeFirst = TRUE, replace = FALSE, params = list(), dump = NULL, debug.level, ...) ## Default S3 method: rtopSim(object = NULL, predictionLocations, varMatObs, varMatPredObs, varMatPred, variogramModel, ...)
object |
object of class |
varMatUpdate |
logical; if TRUE, also existing variance matrices will
be recomputed, if FALSE, only missing variance matrices will be computed,
see also |
beta |
The expected mean of the data, for unconditional simulations |
largeFirst |
Although the simulation method follows a random path around the predictionLocations, simulating the largest area first will assure that the true mean of the simulated values will be closer to beta |
replace |
logical; if observation locations are also present as predictions, should they be replaced? This is particularly when doing conditional simulations for a set of observations with uncertainty. |
params |
a set of parameters, used to modify the standard parameters for
the |
dump |
file name for saving the updated object, after adding variance matrices. Useful if there are problems with the simulation, particularly if it for some reason crashes. |
debug.level |
logical that controls some output, will override the object parameters |
predictionLocations |
|
varMatObs |
covariance matrix of possible observations, where diagonal must consist
of internal variance, typically generated from call
to |
varMatPredObs |
covariance matrix between possible observation locations and
simulation locations, typically generated from call
to |
varMatPred |
covariance matrix between simulation locations, typically generated
from a call to |
variogramModel |
a variogram model of type |
... |
possible modification of the object parameters or default parameters. |
This function can do constrained or unconstrained simulation for areas.
The simplest way of calling the function is with an rtop-object that
contains the fitted variogram model and all the other necessary data (see
createRtopObject or rtop-package). rtopSim
is the only function in rtop which does not need observations.
However, a variogram model is still necessary to perform simulations.
The arguments beta and largeFirst are only used for unconditional simulations.
The function is still in an experimental stage, and might change in the future. There are some issues with the current implementation:
Numerical issues can in some cases give negative estimation variances, which will result in an invalid distribution for the simulation. This will result in simulated NA values for these locations.
The variability of simulated values for small areas (such as small headwater catchments) will be relatively high based on the statistical uncertainty. This could be overestimated compared to the uncertainty which is possible based on rainfall.
If called with SpatialPolygons or sf as predictionLocations
and either SpatialPolygonsDataFrame, sf or NULL for
observations, the function returns aSpatialPolygonsDataFrame or sf with simulations at the
locations defined in predictionLocations
If called with an rtop-object, the function returns the same object with the simulations added to the object.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
# The following command will download the complete example data set # downloadRtopExampleData() rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Setting some parameters; nclus > 1 will start a cluster with nclus # workers for parallel processing params = list(gDist = TRUE, cloud = FALSE, nclus = 1, rresol = 25) # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Build an object rtopObj = createRtopObject(observations, predictionLocations, params = params, formulaString = "obs~1") # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj) # Conditional simulations for two new locations rtopObj10 = rtopSim(rtopObj, nsim = 5) rtopObj11 = rtopObj # Unconditional simulation at the observation locations # (These are moved to the predictionLocations) rtopObj11$predictionLocations = rtopObj11$observations rtopObj11$observations = NULL # Setting varMatUpdate to TRUE, to make sure that covariance # matrices are recomputed rtopObj12 = rtopSim(rtopObj11, nsim = 10, beta = 0.01, varMatUpdate = TRUE) summary(data.frame(rtopObj10$simulations)) summary(data.frame(rtopObj12$simulations))# The following command will download the complete example data set # downloadRtopExampleData() rpath = system.file("extdata",package="rtop") library(sf) observations = st_read(rpath, "observations") predictionLocations = st_read(rpath,"predictionLocations") # Setting some parameters; nclus > 1 will start a cluster with nclus # workers for parallel processing params = list(gDist = TRUE, cloud = FALSE, nclus = 1, rresol = 25) # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM # Build an object rtopObj = createRtopObject(observations, predictionLocations, params = params, formulaString = "obs~1") # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj) # Conditional simulations for two new locations rtopObj10 = rtopSim(rtopObj, nsim = 5) rtopObj11 = rtopObj # Unconditional simulation at the observation locations # (These are moved to the predictionLocations) rtopObj11$predictionLocations = rtopObj11$observations rtopObj11$observations = NULL # Setting varMatUpdate to TRUE, to make sure that covariance # matrices are recomputed rtopObj12 = rtopSim(rtopObj11, nsim = 10, beta = 0.01, varMatUpdate = TRUE) summary(data.frame(rtopObj10$simulations)) summary(data.frame(rtopObj12$simulations))
rtopVariogram will create binned variogram or cloud variogram of data with an areal support.
## S3 method for class 'rtop' rtopVariogram(object, params = list(), ...) ## S3 method for class 'SpatialPolygonsDataFrame' rtopVariogram(object, ...) ## S3 method for class 'SpatialPointsDataFrame' rtopVariogram(object, formulaString, params=list(), cloud, abins, dbins, ...) ## S3 method for class 'STSDF' rtopVariogram(object, formulaString, params=list(), cloud, abins, dbins, data.table = FALSE, ...)## S3 method for class 'rtop' rtopVariogram(object, params = list(), ...) ## S3 method for class 'SpatialPolygonsDataFrame' rtopVariogram(object, ...) ## S3 method for class 'SpatialPointsDataFrame' rtopVariogram(object, formulaString, params=list(), cloud, abins, dbins, ...) ## S3 method for class 'STSDF' rtopVariogram(object, formulaString, params=list(), cloud, abins, dbins, data.table = FALSE, ...)
object |
object of class |
formulaString |
formula that defines the dependent variable as a linear model
of independent variables; suppose the dependent variable has name |
params |
a set of parameters, used to modify the default parameters for
the |
cloud |
logical; if TRUE, calculate the semivariogram cloud, can be used to overrule the cloud parameter in params. |
abins |
possibility to set areal bins (not yet implemented) |
dbins |
possibility to set distance bins (not yet implemented) |
data.table |
an option to use data.table internally for the variogram computation for
|
... |
parameters to other functions called, e.g. gstat's
|
The function creates a variogram, either of type rtopVariogram or
rtopVariogramCloud. This variogram is based on the variogram
function from gstat, but has additional information about the spatial size or
length of the observations. An rtop-object with the variogram added is
returned if the function is called with an rtop-object as argument.
For spatio-temporal objects (STSDF), the variogram is the spatially variogram, averaged for all time steps. There is a possibility to use data.table internally in this function, which can improve computation time for some cases.
The variogram cloud is similar to the variogram cloud from gstat,
with the area/length added to the resulting data.frame. The binned variogram is
also based on the area or length, in addition to the distance between observations.
The bins equally distanced in the log10-space of the distances and areas (lengths).
The size of the bins is decided from the parameters amul and dmul,
defining the number of bins per order of magnitude (1:10, 10:100, and so on).
The distances between areas are in this function based on the centre of gravity.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
## Not run: library(sf) rpath = system.file("extdata",package="rtop") observations = st_read(rpath,"observations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM vario = rtopVariogram(observations, cloud = TRUE) ## End(Not run)## Not run: library(sf) rpath = system.file("extdata",package="rtop") observations = st_read(rpath,"observations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM vario = rtopVariogram(observations, cloud = TRUE) ## End(Not run)
Calibration function which searches a parameter set which is minimizing the value of an objective function
sceua(OFUN, pars, lower, upper, maxn = 10000, kstop = 5, pcento = 0.01, ngs = 5, npg = 2 * length(pars) + 1, nps = length(pars) + 1, nspl = 2 * length(pars) + 1, mings = ngs, iniflg = 1, iprint = 0, iround = 3, peps = 0.0001, plog = rep(FALSE,length(pars)), implicit = NULL, timeout = NULL, ...)sceua(OFUN, pars, lower, upper, maxn = 10000, kstop = 5, pcento = 0.01, ngs = 5, npg = 2 * length(pars) + 1, nps = length(pars) + 1, nspl = 2 * length(pars) + 1, mings = ngs, iniflg = 1, iprint = 0, iround = 3, peps = 0.0001, plog = rep(FALSE,length(pars)), implicit = NULL, timeout = NULL, ...)
OFUN |
A function to be minimized, with first argument the vector of parameters over which minimization is to take place. It should return a scalar result as an indicator of the error for a certain parameter set |
pars |
a vector with the initial guess the parameters |
lower |
the lower boundary for the parameters |
upper |
the upper boundary for the parameters |
maxn |
the maximum number of function evaluations |
kstop |
number of shuffling loops in which the criterion value must change by the given percentage before optimization is terminated |
pcento |
percentage by which the criterion value must change in given number (kstop) of shuffling loops to continue optimization |
ngs |
number of complexes in the initial population |
npg |
number of points in each complex |
nps |
number of points in a sub-complex |
nspl |
number of evolution steps allowed for each complex before complex shuffling |
mings |
minimum number of complexes required, if the number of complexes is allowed to reduce as the optimization proceeds |
iniflg |
flag on whether to include the initial point in population. iniflg = 0, not included. iniflg= 1, included |
iprint |
flag for controlling print-out after each shuffling loop. iprint < 0: no output. iprint = 1: print information on the best point of the population. iprint > 0: print information on every point of the population |
iround |
number of significant digits in print-out |
peps |
convergence level for parameter set (lower number means smaller difference between parameters of the population required for stop) |
plog |
whether optimization should be done in log10-domain. Either a single TRUE value for all parameters, or a vector with TRUE/FALSE for the different parameters |
implicit |
function for implicit boundaries for the parameters (e.g. sum(pars[4]+pars[5]) < 1). See below for details |
timeout |
if different from NULL: maximum time in seconds for execution before the optimization returns with the parameters so far. |
... |
arguments for the objective function, must be named |
sceua is an R-implementation of the Shuffle Complex Evolution - University of Arizona (Duan et al., 1992), a global optimization method which "combines the strengths of the simplex procedure of Nelder and Mead (1965) with the concepts of controlled random search (Price, 1987), competetive evolusion (Holland, 1975)" with the concept of complex shuffling, developed by Duan et al. (1992).
This implementation follows the Fortran implementation relatively close, but adds the possibility of searching in log-space for one or more of the parameters, and it uses the capability of R to pass functions as arguments, making it possible to pass implicit conditions to the parameter selection.
The objective function OFUN is a function which should give an error value for each parameter set.
It should never return non-numeric values such as NA, NULL, or Inf. If some parameter combinations can
give such values, the return value should rather be a large number.
The function works with fixed upper and lower boundaries for the parameters. If the possible range of
a parameter might span several orders of magnitude, it might be better to search in log-space for the optimal parameter,
to reduce the risk of being trapped in local optima. This can be set with the argument plog, which is either
a single value (FALSE/TRUE) or a vector for all parameters.
plog = c(TRUE, FALSE, FALSE, TRUE, TRUE) means that the search for parameters 1,4 and 5 should be in log10-space,
whereas the search for parameters 2 and 3 are in normal space.
Implicit boundaries can be evoked by passing a function implicit to sceua.
This function should give 0 when parameters are acceptable
and 1 if not. If, for example, the condition is that the following sum of parameters four and five should be limited:
sum(pars[4]+pars[5]) <= 1
then the function will be implicit = function(pars) (2*pars[4] + pars[5]) > 1
The function returns a list with the following elements
- a vector of the best parameters combination
- the value of the objective function for this parameter set
- a list of two values
- the function convergence relative to pcento
- the parameter convergence relative to peps
- the number of function evaluations
- the number of shuffling loops
- logical; TRUE if the optimization was aborted because the timeout time was reached, FALSE otherwise
There are also two elements returned as attributes:
- the entire set of parameters from the last evolution step
- the values of the objective function from the last evolution step
The last two can be accessed as attr(sceuares, "parset") and attr(sceuares, "xf"), if the
result is stored as sceuares.
Jon Olav Skoien
Duan, Q., Sorooshian, S., and Gupta, V.K., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28 (4), 1015?1031.
Holland, H.H., 1975. Adaptation in natural and artificial systems, University of Michigan Press, Ann Arbor.
Nelder, J.A. and Mead, R., 1965. A simplex method for function minimization, Comput. J., 7(4), 308-313.
Price, W.L., 1987. Global optimization algorithms for a CAD workstation, J. Optim. Theory Appl., 55(1), 133-146.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
set.seed(1) # generate example data from a function with three parameters # with some random noise fun = function(x, pars) pars[2]*sin(x*pars[1])+pars[3] x = rnorm(50, sd = 3) y = fun(x, pars = c(5, 2, 3)) + rnorm(length(x), sd = 0.3) plot(x,y) # Objective function, summing up squared differences OFUN = function(pars, x, yobs) { yvals = fun(x, pars) sum((yvals-yobs)^2) } sceuares = sceua(OFUN, pars = c(0.1,0.1,0.1), lower = c(-10,0,-10), upper = c(10,10,10), x = x, yobs = y) sceuares xx = seq(min(x), max(x), 0.1) lines(xx, fun(xx, pars = sceuares$par))set.seed(1) # generate example data from a function with three parameters # with some random noise fun = function(x, pars) pars[2]*sin(x*pars[1])+pars[3] x = rnorm(50, sd = 3) y = fun(x, pars = c(5, 2, 3)) + rnorm(length(x), sd = 0.3) plot(x,y) # Objective function, summing up squared differences OFUN = function(pars, x, yobs) { yvals = fun(x, pars) sum((yvals-yobs)^2) } sceuares = sceua(OFUN, pars = c(0.1,0.1,0.1), lower = c(-10,0,-10), upper = c(10,10,10), x = x, yobs = y) sceuares xx = seq(min(x), max(x), 0.1) lines(xx, fun(xx, pars = sceuares$par))
intamap package
This function makes it possible to use rtop-objects in the functions of the package.
It is necessary to load the intamap-package before calling this function.
useRtopWithIntamap()useRtopWithIntamap()
The function will have as side effect that the intamap package is loaded, and that rtop-methods are registered for the intamap-functions estimateParameters, spatialPredict and 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.
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
library(intamap) useRtopWithIntamap()library(intamap) useRtopWithIntamap()
This gives an easier interface to the parameters of the variogram model
rtopVariogramModel(model = "Ex1", sill = NULL, range = NULL, exp = NULL, nugget = NULL, exp0 = NULL, observations = NULL, formulaString = obs~1) ## S3 method for class 'rtop' updateRtopVariogram(object, ...) ## S3 method for class 'rtopVariogramModel' updateRtopVariogram(object, action = "mult", ..., checkVario = FALSE, sampleVariogram = NULL, observations = NULL)rtopVariogramModel(model = "Ex1", sill = NULL, range = NULL, exp = NULL, nugget = NULL, exp0 = NULL, observations = NULL, formulaString = obs~1) ## S3 method for class 'rtop' updateRtopVariogram(object, ...) ## S3 method for class 'rtopVariogramModel' updateRtopVariogram(object, action = "mult", ..., checkVario = FALSE, sampleVariogram = NULL, observations = NULL)
model |
variogram model, currently "Ex1" is the only implemented, see Skoien et al (2006) |
sill |
sill of variogram |
range |
range of variogram |
exp |
the exponent of the fractal part of the variogram, see Skoien et al (2006) |
exp0 |
gives the angle of the first part of the variogram in a log-log plot (weibull type), should be between 0 and 2. See Skoien et al (2006) |
nugget |
nugget of point variogram |
formulaString |
formula that defines the dependent variable as a linear model
of independent variables, see e.g. |
object |
either: object of class |
action |
character variable defining whether the new parameters should
be |
checkVario |
logical, will issue a call to |
sampleVariogram |
a sample variogram of the data |
observations |
a set of observations |
... |
parameters to lower level functions |
The function helps creating and updating the parameters of the variogram, by using common names and simple update methods. This is mainly for manual fitting of the variogram. The automatic call to checkVario makes it easier to visualize the effect of the changes to the variogram
Jon Olav Skoien
## Not run: library(sf) rpath = system.file("extdata",package="rtop") observations = st_read(rpath,"observations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM predictionLocations = st_read(rpath,"predictionLocations") rtopObj = createRtopObject(observations,predictionLocations) # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj) rtopObj = updateRtopVariogram(rtopObj, exp = 1.5, action = "mult", checkVario = TRUE) ## End(Not run)## Not run: library(sf) rpath = system.file("extdata",package="rtop") observations = st_read(rpath,"observations") # Create a column with the specific runoff: observations$obs = observations$QSUMMER_OB/observations$AREASQKM predictionLocations = st_read(rpath,"predictionLocations") rtopObj = createRtopObject(observations,predictionLocations) # Fit a variogram (function also creates it) rtopObj = rtopFitVariogram(rtopObj) rtopObj = updateRtopVariogram(rtopObj, exp = 1.5, action = "mult", checkVario = TRUE) ## End(Not run)
varMat will create a semivariogram matrix between all the supports in a set of locations (observations or prediction locations) or semivariogram matrices between all the supports in one or two sets of locations, and also between them.
## S3 method for class 'rtop' varMat(object, varMatUpdate = FALSE, fullPred = FALSE, params = list(), ...) ## S3 method for class 'SpatialPolygonsDataFrame' varMat(object, object2 = NULL,...) ## S3 method for class 'SpatialPolygons' varMat(object, object2 = NULL, variogramModel, overlapObs, overlapPredObs, ...) ## S3 method for class 'list' varMat(object, object2 = NULL, coor1, coor2, maxdist = Inf, variogramModel, diag = FALSE, sub1, sub2, debug.level = ifelse(interactive(), 1, 0), ...)## S3 method for class 'rtop' varMat(object, varMatUpdate = FALSE, fullPred = FALSE, params = list(), ...) ## S3 method for class 'SpatialPolygonsDataFrame' varMat(object, object2 = NULL,...) ## S3 method for class 'SpatialPolygons' varMat(object, object2 = NULL, variogramModel, overlapObs, overlapPredObs, ...) ## S3 method for class 'list' varMat(object, object2 = NULL, coor1, coor2, maxdist = Inf, variogramModel, diag = FALSE, sub1, sub2, debug.level = ifelse(interactive(), 1, 0), ...)
object |
either: 1) an object of class |
varMatUpdate |
logical; if TRUE, also existing variance matrices will be recomputed, if FALSE, only missing variance matrices will be computed |
fullPred |
logical; whether to create the full covariance matrix also for the predictions, mainly used for simulations |
params |
a set of parameters, used to modify the default parameters for
the |
object2 |
if |
variogramModel |
variogramModel to be used in calculation of the semivariogram matrix (matrices) |
... |
typical parameters to modify from the default parameters of the
rtop-package (or modifications of the previously set parameters for the
|
overlapObs |
matrix with observations that overlap each other |
overlapPredObs |
matrix with |
coor1 |
coordinates of centroids of |
coor2 |
coordinates of centre-of-gravity of |
maxdist |
maximum distance between areas for inclusion in semivariogrma matrix |
diag |
logical; if TRUE only the semivariogram values along the diagonal will be calculated, typical for semivariogram matrix of prediction locations |
sub1 |
semivariogram array for subtraction of inner variances of areas |
sub2 |
semivariogram array for subtraction of inner variances of areas |
debug.level |
debug.level >= 1 will give output for every element |
The lower level versions of the function calculates a semivariogram matrix
between locations in object or between the locations in object
and the locations in object2. The method for object of type rtop
calculates semivariogram matrices between observation locations, between prediction locations,
and between observation locations and prediction locations, and adds these
to object.
The argument varMatUpdate is typically used to avoid repeated computations
of the same variance matrices. Default is FALSE, which will avoid recomputation
of the variance matrix for the observations if the procedure is cross-validation
before interpolation. Should be set to TRUE if the variogram Model has been
changed, or if observation and/or prediction locations have been changed.
If an rtop-object contains observations and/or predictionLocations of type
STSDF, the covariance matrix is computed based on the spatial
properties of the object.
Jon Olav Skoien
Skoien J. O., R. Merz, and G. Bloschl. Top-kriging - geostatistics on stream networks. Hydrology and Earth System Sciences, 10:277-287, 2006.
Skoien, J. O., Bloschl, G., Laaha, G., Pebesma, E., Parajka, J., Viglione, A., 2014. Rtop: An R package for interpolation of data with a variable spatial support, with an example from river networks. Computers & Geosciences, 67.
## Not run: library(sf) rpath = system.file("extdata",package="rtop") observations = st_read(rpath,"observations") vmod = list(model = "Ex1", params = c(0.00001,0.007,350000,0.9,1000)) vm = varMat(observations, variogramModel = vmod) ## End(Not run)## Not run: library(sf) rpath = system.file("extdata",package="rtop") observations = st_read(rpath,"observations") vmod = list(model = "Ex1", params = c(0.00001,0.007,350000,0.9,1000)) vm = varMat(observations, variogramModel = vmod) ## End(Not run)