Version:	0.11.3
Date:	2023-04-17
Title:	Geostatistics for Compositional Analysis
Depends:	R (≥ 2.10)
Suggests:	FNN, JADE, jointDiag, DescTools, knitr, rmarkdown (≥ 2.3), magrittr, readxl
Imports:	methods, gstat, compositions (≥ 2.0), sp, boot, foreach, utils, RColorBrewer
Description:	Support for geostatistical analysis of multivariate data, in particular data with restrictions, e.g. positive amounts, compositions, distributional data, microstructural data, etc. It includes descriptive analysis and modelling for such data, both from a two-point Gaussian perspective and multipoint perspective. The methods mainly follow Tolosana-Delgado, Mueller and van den Boogaart (2018) <doi:10.1007/s11004-018-9769-3>.
License:	CC BY-SA 4.0 | GPL-2 | GPL-3 [expanded from: CC BY-SA 4.0 | GPL (≥ 2)]
URL:	https://codebase.helmholtz.cloud/geomet/gmGeostats
Copyright:	(C) 2020 by Helmholtz-Zentrum Dresden-Rossendorf and Edith Cowan University; gsi.DS code by Hassan Talebi
RoxygenNote:	7.1.1
Encoding:	UTF-8
VignetteBuilder:	knitr
LazyData:	true
Collate:	'Anamorphosis.R' 'preparations.R' 'gstatCompatibility.R' 'gmAnisotropy.R' 'compositionsCompatibility.R' 'variograms.R' 'gmSpatialMethodParameters.R' 'abstractClasses.R' 'accuracy.R' 'closeup.R' 'data.R' 'exploratools.R' 'genDiag.R' 'geostats.R' 'gmDataFrameStack.R' 'gmSimulation.R' 'gmSpatialModel.R' 'gmValidationStrategy.R' 'grids.R' 'mask.R' 'spSupport.R' 'uncorrelationTest.R' 'zzz.R'
NeedsCompilation:	yes
Packaged:	2023-04-17 14:25:39 UTC; tolosa53
Author:	Raimon Tolosana-Delgado [aut], Ute Mueller [aut], K. Gerald van den Boogaart [ctb, cre], Hassan Talebi [ctb, cph], Helmholtz-Zentrum Dresden-Rossendorf [cph], Edith Cowan University [cph]
Maintainer:	K. Gerald van den Boogaart <support@boogaart.de>
Repository:	CRAN
Date/Publication:	2023-04-18 11:50:12 UTC

Combination of gmCgram variogram structures

Description

combination of nested structures of a gmCgram object

Usage

## S3 method for class 'gmCgram'
x + y

Arguments

Value

The combined nested structures

Examples

utils::data("variogramModels")
v1 = setCgram(type=vg.Gau, sill=diag(2), anisRanges = 3*diag(c(3,1)))
v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
vm = v1+v2

Force a matrix to be anisotropy range matrix,

Description

Force a matrix M to be considered an anisotropy range matrix, i.e with ranges and orientations, such that u = sqrt(h' * M^{-1} * h) allows to use an isotropic variogram.

Usage

AnisotropyRangeMatrix(x, checkValidity = TRUE)

## S3 method for class 'AnisotropyScaling'
as.AnisotropyRangeMatrix(x)

Arguments

Value

the same matrix with a class attribute

Methods (by generic)

• as.AnisotropyRangeMatrix: Convert from AnisotropyScaling

See Also

Other anisotropy: AnisotropyScaling(), anis_GSLIBpar2A(), as.AnisotropyRangeMatrix(), as.AnisotropyScaling(), is.anisotropySpecification()

Convert to anisotropy scaling matrix

Description

Convert an anisotropy specification to a scaling matrix

Usage

AnisotropyScaling(x)

Arguments

Value

An anisotropy scaling matrix A is such that for any lag vector h, the variogram model turns isotropic in terms of u'=h'\cdot A. This function does not check any special property for this matrix! You should probably be using anis_GSLIBpar2A() isntead, and leave AnisotropyScaling() for internal uses.

See Also

Other anisotropy: AnisotropyRangeMatrix(), anis_GSLIBpar2A(), as.AnisotropyRangeMatrix(), as.AnisotropyScaling(), is.anisotropySpecification()

Examples

( l = anis_GSLIBpar2A(angles=30, ratios=0.5))
( ll = unclass(l) )
AnisotropyScaling(l)

Create a parameter set specifying a LU decomposition simulation algorithm

Description

Create a parameter set describing a Cholesky (or LU) decomposition algorithm to two-point simulation, mostly for covariance or variogram-based gaussian random fields.

Usage

CholeskyDecomposition(nsim = 1, seed = NULL, ...)

Arguments

Value

an S3-list of class "gmCholeskyDecomposition" containing the few elements given as arguments to the function. This is just a compact way to provide further functions such as predict_gmSpatialModel with appropriate triggers for choosing a prediction method or another, in this case for triggering LU or Cholesky decomposition simulation.

Examples

(chols_local = CholeskyDecomposition(nsim=100, nBands=300))
## then run predict(..., pars=chols_local)

Create a parameter set specifying a direct sampling algorithm

Description

Create a parameter set describing a direct sampling algorithm to multipoint simulation. All parameters except nsim are optional, as they have default values reasonable according to experience.

Usage

DSpars(
  nsim = 1,
  scanFraction = 0.25,
  patternSize = 10,
  gof = 0.05,
  seed = NULL,
  ...
)

Arguments

Value

an S3-list of class "gmDirectSamplingParameters" containing the six elements given as arguments to the function. This is just a compact way to provide further functions such as predict_gmSpatialModel with appropriate triggers for choosing a prediction method or another, in this case for triggering direct sampling.

Examples

(dsp = DSpars(nsim=100, scanFraction=75, patternSize=6, gof=0.05))
## then run predict(..., pars=dsp)

Create a data frame stack

Description

Make a stacked data frame, that is, a stack of data.frames representing e.g. repeated measurements, parallel time series, or a stack of multivariate realisation of some random process. If a data frame is analogous to a matrix, a DataFrameStack is analogous to an array. It is highly recommendable to work with named dimensions in stacked data frames.

Usage

DataFrameStack(x, ...)

as.DataFrameStack(x, ...)

## S3 method for class 'data.frame'
DataFrameStack(x, stackDim = 2, dim = NULL, Dimnames = NULL, ...)

## S3 method for class 'data.frame'
as.DataFrameStack(x, stackDim = 2, dim = NULL, Dimnames = NULL, ...)

## S3 method for class 'array'
DataFrameStack(x, stackDim = 2, ...)

## S3 method for class 'array'
as.DataFrameStack(x, stackDim = 2, ...)

## S3 method for class 'list'
DataFrameStack(x, stackDimName = NULL, Dimnames = NULL, ...)

## S3 method for class 'list'
as.DataFrameStack(x, stackDimName = NULL, Dimnames = NULL, ...)

Arguments

Value

The data provided reorganised as a DataFrameStack, with several additional attributes.

Functions

• DataFrameStack: create a DataFrameStack from an array

• as.DataFrameStack: create a DataFrameStack from an array

• as.DataFrameStack.data.frame: create a DataFrameStack from an array

• DataFrameStack.array: create a DataFrameStack from an array

• as.DataFrameStack.array: create a DataFrameStack from an array

• DataFrameStack.list: create a DataFrameStack from a list

• as.DataFrameStack.list: create a DataFrameStack from an array

See Also

stackDim() to get or set the stacking dimension; noStackDim() to get (not set) the non-stacking dimension; as.array.DataFrameStack() and as.list.DataFrameStack() to convert the stack to an array or a list; dimnames.DataFrameStack() to get the dimension names; [⁠[.DataFrameStack⁠] to extract rows of a stack; getStackElement() and setStackElement (same page as getStackElement) to extract/modify an element of the stack; gmApply() to apply any function to the stack, typically element-wise; and swarmPlot() to combine plot elements for each stack element into a single plot.

Examples

ll = lapply(1:3, function(i) as.data.frame(matrix(1:10+(i-1)*10, ncol=2)))
dfs = DataFrameStack(ll, stackDimName="rep")
dimnames(dfs)
df = as.data.frame(matrix(1:30, ncol=6))
dfs = DataFrameStack(df, dimnames = list(obs=1:5, vars=c("A","B"), rep=1:3), stackDim = "rep")
dimnames(dfs)
ar = array(1:30, dim = c(5,2,3), dimnames=list(obs=1:5, vars=c("A","B"), rep=1:3))
dfs = DataFrameStack(ar, stackDim="rep")
dimnames(dfs)

Empirical structural function specification

Description

Abstract class, containing any specification of an empirical variogram (or covariance function, or variations). Members must implement a coercion method to class "gmEVario" (see gsi.EVario2D() for an example), and (possibly) coercion to class "gstatVariogram" (see gstat::variogram())

Superclass for grid or nothing

Description

Superclass for slots containing a grid topology or being empty

Create a parameter set of local for neighbourhood specification.

Description

Create a parameter set describing a kriging neighbourhood (local or global) for cokriging and cokriging based simulation. This heavily relies on the definitions of gstat::gstat(). All parameters are optional, as their default amounts to a global neihghbourhood.

Usage

KrigingNeighbourhood(
  nmax = Inf,
  nmin = 0,
  omax = 0,
  maxdist = Inf,
  force = FALSE,
  anisotropy = NULL,
  ...
)

Arguments

Value

an S3-list of class "gmKrigingNeighbourhood" containing the six elements given as arguments to the function. This is just a compact way to provide further functions such as predict_gmSpatialModel with appropriate triggers for choosing a prediction method or another, in this case for triggering cokriging (if alone) or eventually sequential simulation (see SequentialSimulation()).

Examples

data("jura", package="gstat")
X = jura.pred[,1:2]
summary(X)
Zc = jura.pred[,7:10]
ng_global = KrigingNeighbourhood()
ng_local = KrigingNeighbourhood(maxdist=1, nmin=4, 
                                omax=5, force=TRUE)
ng_local
ng_global
make.gmCompositionalGaussianSpatialModel(data = Zc, coords = X,
                                         V = "alr", ng = ng_local)

Create a anisotropic model for regionalized compositions

Description

Creates a (potentially anisotropic) variogram model for variation-variograms

Usage

LMCAnisCompo(
  Z,
  models = c("nugget", "sph", "sph"),
  azimuths = rep(0, length(models)),
  ranges = matrix(1, nrow = length(models), ncol = G),
  sillarray = NULL,
  V = NULL,
  tol = 1e-12,
  G = 2
)

Arguments

Value

an object of class "LMCAnisCompo" with all provided information appropriately structured, ready to be used for fitting, plotting or prediction.

Examples

data("jura", package="gstat")
Zc = compositions::acomp(jura.pred[,7:9])
LMCAnisCompo(Zc, models=c("nugget", "sph"), azimuths=c(0,45))

Specify the leave-one-out strategy for validation of a spatial model

Description

Function to specify the leave-one-out strategy for validation of a spatial model

Usage

LeaveOneOut()

Value

an object of class c("LeaveOneOut", "NfoldCrossValidation") to be used in a call to validate()

See Also

Other validation functions: NfoldCrossValidation, validate()

Examples

LeaveOneOut()

Generalised diagonalisations Calculate several generalized diagonalisations out of a data set and its empirical variogram

Description

Generalised diagonalisations Calculate several generalized diagonalisations out of a data set and its empirical variogram

Usage

Maf(x, ...)

## S3 method for class 'data.frame'
Maf(x, vg, i = 2, ...)

## S3 method for class 'rmult'
Maf(x, vg, i = 2, ...)

## S3 method for class 'aplus'
Maf(x, vg, i = 2, ...)

## S3 method for class 'rplus'
Maf(x, vg, i = 2, ...)

## S3 method for class 'ccomp'
Maf(x, vg, i = 2, ...)

## S3 method for class 'rcomp'
Maf(x, vg, i = 2, ...)

## S3 method for class 'acomp'
Maf(x, vg, i = 2, ...)

UWEDGE(x, ...)

## Default S3 method:
UWEDGE(x, ...)

## S3 method for class 'acomp'
UWEDGE(x, vg, i = NULL, ...)

## S3 method for class 'rcomp'
UWEDGE(x, vg, i = NULL, ...)

RJD(x, ...)

## Default S3 method:
RJD(x, ...)

## S3 method for class 'acomp'
RJD(x, vg, i = NULL, ...)

## S3 method for class 'rcomp'
RJD(x, vg, i = NULL, ...)

Arguments

Value

An object extending c("princomp.CLASSOF(x)","princomp") with classes "genDiag", and an extra class argument depending on the diagonalisation method chosen.

Function Maf results carry the extra class "maf", and they correspond to minimum/maximum autocorrelation factors (MAF) as proposed by Switzer and Green (1984). In this case, the slicer is typically just the index of one lag distance (defaults to i=2). MAF provides the analytical solution to the joint diagonalisation of two matrices, the covariance of increments provided by the slicing and the conventional covariance matrix (of the idt transformed values, if appropriate). Resulting factors are ordered in decreasing order of spatial continuity, which might produce surprising scree-plots for those who are used to see a screeplot of a principal component analysis.

Function UWEDGE (Uniformly Weighted Exhaustive Diagonalization with Gauss iterations; Tichavsky and Yeredor, 2009) results carry the extra class "uwedge". The function is a wrapper on jointDiag::uwedge from package jointDiag (Gouy-Pailler, 2017). In this case the slicer is typically just a subset of indices of lag distances to consider (defaults to the nearest indexes to mininum, maximum and mean lag distances of the variogram). UWEDGE iteratively seeks for a pair of matrices (a mixing and a demixing matrices) diagonalises the set of matrices M_1, M_2, \ldots, M_K given, by minimizing the target quantity

Q_{uwedge}(A, W) = \sum_{k=1}^K Tr[E_k^t\cdot E_k],

where E_k = (W^t\cdot M_k \cdot W- A\cdot \Lambda_k\cdot A^t) and \Lambda_k = diag(W^t\cdot M_k \cdot W) is the resulting diagonalized version of each matrix. Obtained factors are ordered in decreasing order of spatial continuity, which might produce surprising scree-plots for those who are used to see a screeplot of a principal component analysis.

Function RJD results carry the extra class "rjd". The function is a wrapper on JADE::rjd (Miettinen, Nordhausen and Taskinen, 2017), implementing the Rotational joint diagonalisation method (Cardoso and Souloumiac, 1996). In this case the slicer is typically just a subset of indices of lag distances to consider (defaults to the nearest indexes to mininum, maximum and mean lag distances). This algorithm also served for quasi-diagonalising a set of matrices as in UWEDGE, just that in this case the quantity to minimise is the sum of sequares of all off-diagonal elements of A^t\cdot M_k\cdot A for all k=1, 2, \ldots K.

All these functions produce output mimicking princomp, i.e. with elements

sdev

contrary to the output in PCA, this contains the square root of the metric variance of the predictions obtained for each individual factor; this is the quantity needed for screeplot to create plots of explained variance by factor

loadings

matrix of contributions of each (cdt-transformed) original variable to the new factors

center

center of the data set (eventually, represented through cdt), in compositional methods

scale

the scalings applied to each original variable

n.obs

number of observations

scores

the scores of the supplied data on the new factors

call

the call to the function (attention: it actually may come much later)

and additionally some of the following arguments, in different order

invLoadings

matrix of contributions of each factor onto each original variable

Center

compositional methods return here the cdt-backtransformed center

InvLoadings

compositional methods return here the clr-backtransformed inverse loadings, so that each column of this matrix can be understood as a composition on itself

DownInvLoadings

compositional methods return here the clr-backtransformed "minus inverse loadings", so that each column of this matrix can be understood as a composition on itself; details in princomp.acomp

C1, C2

Maf returns the two matrices that were diagonalised

eigenvalues

Maf returns the generalized eigenvalues of the diagonalisation of C1 and C2

gof

UWEDGE returns the values of the goodness of fit criterion across sweeps

diagonalized

RJD returns the diagonalized matrices, in an array of (K,D,D)-dimensions, being D the number of variables in x

type

a string describing which package and which function was used as a workhorse for the calculation

NOTE: if the arguments you provide to RJD and UWEDGE are not of the appropriate type (i.e. data.frames or equivalent) the default method of these functions will just attempt to call the basis functions JADE:rjd and JointDiag:uwedge respectively. This will be the case if you provide x as a "matrix", or as an "array". In those cases, the output will NOT be structured as an extension to princomp results; instead they will be native output from those functions.

Functions

• Maf: Generalised diagonalisations

• Maf.rmult: Generalised diagonalisations

• Maf.aplus: Generalised diagonalisations

• Maf.rplus: Generalised diagonalisations

• Maf.ccomp: Generalised diagonalisations

• Maf.rcomp: Generalised diagonalisations

• Maf.acomp: Generalised diagonalisations

• UWEDGE: Generalised diagonalisations

• UWEDGE.default: Generalised diagonalisations

• UWEDGE.acomp: Generalised diagonalisations

• UWEDGE.rcomp: Generalised diagonalisations

• RJD: Generalised diagonalisations

• RJD.default: Generalised diagonalisations

• RJD.acomp: Generalised diagonalisations

• RJD.rcomp: Generalised diagonalisations

References

Cardoso, J. K. and Souloumiac A. 1996. Jacobi angles for simultaneous diagonalization. SIAM Journal of Matrix Analysis and Applications 17(1), 161-164.

Gouy-Pailler C., 2017. jointDiag: Joint approximate diagonalization of a set of square matrices. R package version 0.3. https://CRAN.R-project.org/package=jointDiag

Miettinen J., Nordhausen K., and Taskinen, S., 2017. Blind source separation based on Joint diagonalization in R: The packages JADE and BSSasymp. Journal of Statistical Software 76(2), 1-31.

Switzer P. and Green A.A., 1984. Min/Max autocorrelation factors for multivariate spatial imaging, Stanford University, Palo Alto, USA, 14pp.

Tichavsky, P. and Yeredor, A., 2009. Fast approximate joint diagonalization incorporating weight matrices. IEEE Transactions on Signal Processing 57, 878 ??? 891.

See Also

Other generalised Diagonalisations: coloredBiplot.genDiag(), predict.genDiag()

Examples

require("magrittr")
require("gstat")
require("compositions")
data("jura", package="gstat")
gs = gstat(id="Cd", formula=log(Cd)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
gstat(id="Co", formula=log(Cd)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
gstat(id="Cr", formula=log(Cr)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
gstat(id="Cu", formula=log(Cu)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
gstat(id="Ni", formula=log(Ni)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
gstat(id="Pb", formula=log(Pb)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
gstat(id="Zn", formula=log(Zn)~1, locations=~Xloc+Yloc, data=jura.pred)
vg = variogram(gs)
mf = Maf(aplus(jura.pred[, -(1:6)]), vg=vg)
mf
mf$loadings
biplot(mf)

Structural function model specification

Description

Abstract class, containing any specification of a variogram (or covariance) model. Members must implement a coercion method to class "gmCgram" (see setCgram() for an example), and (possibly) coercion to class "variogramModel" or "variogramModelList" (see gstat::vgm())

National Geochemical Survey of Australia: soil data

Description

A dataset containing the geochemical composition of soil samples covering approx 2/3 of Australia, at a spatial resolution of 1 sample/5000 $km^2$. The original data is published by Geosciences Australia under the Creative Commons CC4 Licence + Attribution (CC-BY-SA-4.0). The provided data set has additional information and some pre-treatment with respect to the one available in the references below. These changes:

Usage

NGSAustralia

Format

An object of class tbl_df (inherits from tbl, data.frame) with 5259 rows and 76 columns.

Details

• Concentrations are all reported as mg/kg, also known as parts per million (ppm). Values below detection limit are reported as minus the detection limit, as standardised in package "compositions".

• The samples were complemented with information about their belonging to one of the major crustal blocks (MCB) of Australia.

• Easting and Northing coordinates were computed using the Lambert conformal conic projection of Australia (earth ellipsoid GRS80; standard parallels: 18S and 36S latitude; central meridian: 134E longitude).

#' @format A tibble, a data set of compound class c("tbl_df", "tbl", "data.frame") with 5259 observations and 76 variables:

ORDER

Entry order

SITEID

ID of the sampling site

DATE SAMPLED

Timestamp of the sampling

EAST

Easting coordinate of the sampling site

NORTH

Northing coordinate of the sampling site

LATITUDE

Latitude of the sampling site

LONGITUDE

Longitude of the sampling site

STATE

State of the sampling site, one of: "NWS" (New South Wales), "NT" (Northern Territory), "QLD" (Qeensland), "SA" (Southern Australia), "TAS" (Tasmania), "VIC" (Victoria) or "WA" (Western Australia)

REGION

One of the three geo-regions of the data set: "EAST", "WEST" or "EUCLA"

DUPLICATE CODE

Marker for duplicate samples, for quality control

SAMPLEID

ID of the sample

GRAIN SIZE

Particle size of the soil sample, one of "<2 mm" (coarse) or "<75 um" (fine)

DEPTH

Sampling depth, one of: "TOS" (top soil) or "BOS" (bottom soil)

CODE

A combination of ⁠Grain Size⁠ and DEPTH, one of "Tc" "Tf" "Bc" "Bf", standing for TOS coarse, TOS fine, BOS coarse and BOS fine, respectively

⁠Ag ICP-MS mg/kg 0.03⁠

concemtration of Silver in that sample, analysed with inductive coupled plasma mass spectrometry (ICP-MS), with a detection limit of 0.03ppm

⁠Al XRF mg/kg 26⁠

concentration of Aluminium, analysed with X-Ray Fluorescence (XRF), detection limit 26 ppm

⁠As ICP-MS mg/kg 0.4⁠

concentration of Arsenic measured with ICP-MS , detection limit: 0.4 mg/kg

⁠Au FA mg/kg 0.001⁠

concentration of Gold measured with fire assay (FA) , detection limit: 0.001 mg/kg

⁠Ba ICP-MS mg/kg 0.5⁠

concentration of Barium measured with ICP-MS , detection limit: 0.5 mg/kg

⁠Be ICP-MS mg/kg 1.1⁠

concentration of Berillium measured with ICP-MS , detection limit: 1.1 mg/kg

⁠Bi ICP-MS mg/kg 0.02⁠

concentration of Bismuth measured with ICP-MS , detection limit: 0.02 mg/kg

⁠Ca XRF mg/kg 14⁠

concentration of Calcium measured with XRF , detection limit: 14 mg/kg

⁠Cd ICP-MS mg/kg 0.1⁠

concentration of Cadmium measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Ce ICP-MS mg/kg 0.03⁠

concentration of Cerium measured with ICP-MS , detection limit: 0.03 mg/kg

⁠Cl XRF mg/kg 10⁠

concentration of Chlorine measured with XRF , detection limit: 10 mg/kg

⁠Co ICP-MS mg/kg 0.1⁠

concentration of Cobalt measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Cr ICP-MS mg/kg 0.5⁠

concentration of Cromium measured with ICP-MS , detection limit: 0.5 mg/kg

⁠Cs ICP-MS mg/kg 0.1⁠

concentration of Caesium measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Cu ICP-MS mg/kg 0.2⁠

concentration of Copper measured with ICP-MS , detection limit: 0.2 mg/kg

⁠Dy ICP-MS mg/kg 0.1⁠

concentration of Dysprosium measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Er ICP-MS mg/kg 0.03⁠

concentration of Erbium measured with ICP-MS , detection limit: 0.03 mg/kg

⁠Eu ICP-MS mg/kg 0.03⁠

concentration of Europium measured with ICP-MS , detection limit: 0.03 mg/kg

⁠F ISE mg/kg 20⁠

concentration of Fluoride measured with Ion Selective Electrode (ISE), detection limit: 20 mg/kg

⁠FeT XRF mg/kg 35⁠

concentration of total Iron measured with XRF , detection limit: 35 mg/kg

⁠Ga ICP-MS mg/kg 0.1⁠

concentration of Gallium measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Gd ICP-MS mg/kg 0.03⁠

concentration of Gadolinium measured with ICP-MS , detection limit: 0.03 mg/kg

⁠Ge ICP-MS mg/kg 0.04⁠

concentration of Germanium measured with ICP-MS , detection limit: 0.04 mg/kg

⁠Hf ICP-MS mg/kg 0.04⁠

concentration of Hafnium measured with ICP-MS , detection limit: 0.04 mg/kg

⁠Ho ICP-MS mg/kg 0.02⁠

concentration of Holmium measured with ICP-MS , detection limit: 0.02 mg/kg

⁠K XRF mg/kg 42⁠

concentration of Potassium measured with XRF , detection limit: 42 mg/kg

⁠La ICP-MS mg/kg 0.1⁠

concentration of Lantanum measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Lu ICP-MS mg/kg 0.02⁠

concentration of Lutetium measured with ICP-MS , detection limit: 0.02 mg/kg

⁠Mg XRF mg/kg 60⁠

concentration of Magnesium measured with XRF , detection limit: 60 mg/kg

⁠Mn XRF mg/kg 39⁠

concentration of Manganese measured with XRF , detection limit: 39 mg/kg

⁠Mo ICP-MS mg/kg 0.3⁠

concentration of Molybdenum measured with ICP-MS , detection limit: 0.3 mg/kg

⁠Na XRF mg/kg 74⁠

concentration of Sodium measured with XRF , detection limit: 74 mg/kg

⁠Nb ICP-MS mg/kg 0.03⁠

concentration of Niobium measured with ICP-MS , detection limit: 0.03 mg/kg

⁠Nd ICP-MS mg/kg 0.1⁠

concentration of Neodymium measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Ni ICP-MS mg/kg 0.5⁠

concentration of Nickel measured with ICP-MS , detection limit: 0.5 mg/kg

⁠P XRF mg/kg 22⁠

concentration of Phosphorus measured with XRF , detection limit: 22 mg/kg

⁠Pb ICP-MS mg/kg 0.1⁠

concentration of Lead measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Pd FA mg/kg 0.001⁠

concentration of Palladium measured with FA , detection limit: 0.001 mg/kg

⁠Pr ICP-MS mg/kg 0.02⁠

concentration of Praseodymium measured with ICP-MS , detection limit: 0.02 mg/kg

⁠Pt FA mg/kg 0.0005⁠

concentration of Platinum measured with FA , detection limit: 0.0005 mg/kg

⁠Rb ICP-MS mg/kg 0.2⁠

concentration of Rubidium measured with ICP-MS , detection limit: 0.2 mg/kg

⁠S XRF mg/kg 10⁠

concentration of Sulphur measured with XRF , detection limit: 10 mg/kg

⁠Sb ICP-MS mg/kg 0.4⁠

concentration of Antimony measured with ICP-MS , detection limit: 0.4 mg/kg

⁠Sc ICP-MS mg/kg 0.3⁠

concentration of Scandium measured with ICP-MS , detection limit: 0.3 mg/kg

⁠Si XRF mg/kg 47⁠

concentration of Silicon measured with XRF , detection limit: 47 mg/kg

⁠Sm ICP-MS mg/kg 0.04⁠

concentration of Samarium measured with ICP-MS , detection limit: 0.04 mg/kg

⁠Sn ICP-MS mg/kg 0.2⁠

concentration of Tin measured with ICP-MS , detection limit: 0.2 mg/kg

⁠Sr ICP-MS mg/kg 0.2⁠

concentration of Strontium measured with ICP-MS , detection limit: 0.2 mg/kg

⁠Ta ICP-MS mg/kg 0.02⁠

concentration of Tantalum measured with ICP-MS , detection limit: 0.02 mg/kg

⁠Tb ICP-MS mg/kg 0.02⁠

concentration of Terbium measured with ICP-MS , detection limit: 0.02 mg/kg

⁠Th ICP-MS mg/kg 0.02⁠

concentration of Thorium measured with ICP-MS , detection limit: 0.02 mg/kg

⁠Ti XRF mg/kg 30⁠

concentration of Titanium measured with XRF , detection limit: 30 mg/kg

⁠U ICP-MS mg/kg 0.02⁠

concentration of Uranium measured with ICP-MS , detection limit: 0.02 mg/kg

⁠V ICP-MS mg/kg 0.1⁠

concentration of Vanadium measured with ICP-MS , detection limit: 0.1 mg/kg

⁠W ICP-MS mg/kg 0.1⁠

concentration of Wolframium measured with ICP-MS , detection limit: 0.1 mg/kg

⁠Y ICP-MS mg/kg 0.05⁠

concentration of Yttrium measured with ICP-MS , detection limit: 0.05 mg/kg

⁠Yb ICP-MS mg/kg 0.04⁠

concentration of Ytterbium measured with ICP-MS , detection limit: 0.04 mg/kg

⁠Zn ICP-MS mg/kg 0.9⁠

concentration of Zink measured with ICP-MS , detection limit: 0.9 mg/kg

⁠Zr ICP-MS mg/kg 0.2⁠

concentration of Zirconium measured with ICP-MS , detection limit: 0.2 mg/kg

#'

LOI CALC mg/kg

concentration of Loss-on-Ignition, measured by calcination

⁠LOI CALC mg/kg⁠

concentration of Loss-On-Ignition measured with Calcination , detection limit: 0.2 mg/kg

#'

LOI CALC mg/kg

concentration of Loss-on-Ignition, measured by calcination

MRB

obtained from simplifying the major crustal boundaries of Korsch2015a, Korsch2015b

License

CC BY-SA 4.0

Source

https://www.ga.gov.au/about/projects/resources/national-geochemical-survey

Specify a strategy for validation of a spatial model

Description

Specify a strategy to validate a spatial model. Currently only leave-one-out and n-fold cross-validation are available, each specified by its own function. Leave-one-out takes no parameter.

Usage

NfoldCrossValidation(nfolds = 2, doAll = TRUE, ...)

Arguments

Value

An object, a list with an appropriate class, controlling the strategy specified. This can be of class "NfoldCrossValidation" or of class c("LeaveOneOut", "NfoldCrossValidation").

See Also

Other validation functions: LeaveOneOut, validate()

Examples

NfoldCrossValidation(nfolds=5, doAll=FALSE)

Predict method for objects of class 'gmSpatialModel'

Description

This is a one-entry function for several spatial prediction and simulation methods, for model objects of class gmSpatialModel. The several methods are chosen by means of pars objects of the appropriate class.

Usage

Predict(object, newdata, pars, ...)

predict(object, ...)

## S3 method for class 'gmSpatialModel'
predict(object, newdata, pars = NULL, ...)

## S4 method for signature 'gmSpatialModel'
predict(object, newdata, pars = NULL, ...)

## S4 method for signature 'gmSpatialModel,ANY,ANY'
Predict(object, newdata, pars, ...)

## S4 method for signature 'gmSpatialModel,ANY,gmNeighbourhoodSpecification'
Predict(object, newdata, pars, ...)

## S4 method for signature 'gmSpatialModel,ANY,gmTurningBands'
Predict(object, newdata, pars, ...)

## S4 method for signature 'gmSpatialModel,ANY,gmCholeskyDecomposition'
Predict(object, newdata, pars, ...)

## S4 method for signature 'gmSpatialModel,ANY,gmSequentialSimulation'
Predict(object, newdata, pars, ...)

## S4 method for signature 'gmSpatialModel,ANY,gmDirectSamplingParameters'
Predict(object, newdata, pars, ...)

Arguments

Details

Package "gmGeostats" aims at providing a broad series of algorithms for geostatistical prediction and simulation. All can be accesses through this interface, provided that arguments object and pars are of the appropriate kind. In object, the most important criterion is the nature of its slot model. In pars its class counts: for the creation of informative parameters in the appropriate format and class, a series of accessory functions are provided as well.

Classical (gaussian-based two-point) geostatistics are obtained if object@model contains a covariance function, or a variogram model. Argument pars can be created with functions such as KrigingNeighbourhood(), SequentialSimulation(), TurningBands() or CholeskyDecomposition() to respectively trigger a cokriging, as sequential Gaussian simulation, a turning bands simulation, or a simulation via Cholesky decomposition. The kriging neighbourhood can as well be incorporated in the "gmSpatialModel" object directly, or even be nested in a "SequentialSimulation" parameter object.

Conversely, to run a multipoint geostatistics algorithm, the first condition is that object@model contains a training image. Additionally, pars must describe the characteristics of the algorithm to use. Currently, only direct sampling is available: it can be obtained by providing some parameter object created with a call to DirectSamplingParameters(). This method requires newdata to be on a gridded set of locations (acceptable object classes are sp::gridTopology, sp::SpatialGrid, sp::SpatialPixels, sp::SpatialPoints or data.frame, for the last two a forced conversion to a grid will be attempted).

Value

Depending on the nature of newdata, the result will be a data container of the same kind, extended with the predictions or simulations. For instance, if we want to obtain predictions on the locations of a "SpatialPoints", the result will be a sp::SpatialPointsDataFrame(); if we want to obtain simulations on the coordinates provided by a "data.frame", the result will be a DataFrameStack() with the spatial coordinates stored as an extra attribute; or if the input for a simulation is a masked grid of class sp::SpatialPixels(), the result will be of class sp::SpatialPixelsDataFrame() which data slot will be a DataFrameStack.

See Also

Other gmSpatialModel: as.gmSpatialModel(), gmSpatialModel-class, make.gmCompositionalGaussianSpatialModel(), make.gmCompositionalMPSSpatialModel(), make.gmMultivariateGaussianSpatialModel()

Create a parameter set specifying a gaussian sequential simulation algorithm

Description

Create a parameter set describing a sequential simulation algorithm to two-point simulation, mostly for covariance or variogram-based gaussian random fields.

Usage

SequentialSimulation(
  nsim = 1,
  ng = NULL,
  rank = Inf,
  debug.level = 1,
  seed = NULL,
  ...
)

Arguments

Value

an S3-list of class "gmSequentialSimulation" containing the four elements given as arguments to the function. This is just a compact way to provide further functions such as predict_gmSpatialModel with appropriate triggers for choosing a prediction method or another, in this case for triggering sequential Gaussian simulation.

Examples

data("jura", package="gstat")
X = jura.pred[,1:2]
summary(X)
Zc = jura.pred[,7:10]
ng_local = KrigingNeighbourhood(maxdist=1, nmin=4, omax=5, force=TRUE)
(sgs_local = SequentialSimulation(nsim=100, ng=ng_local, debug.level=-1))
## then run predict(..., pars=sgs_local)

Create a parameter set specifying a turning bands simulation algorithm

Description

Create a parameter set describing a turning bands algorithm to two-point simulation, mostly for covariance or variogram-based gaussian random fields.

Usage

TurningBands(nsim = 1, nBands = 1000, seed = NULL, ...)

Arguments

Value

an S3-list of class "gmTurningBands" containing the few elements given as arguments to the function. This is just a compact way to provide further functions such as predict_gmSpatialModel with appropriate triggers for choosing a prediction method or another, in this case for triggering turning bands simulation.

Examples

(tbs_local = TurningBands(nsim=100, nBands=300))
## then run predict(..., pars=tbs_local)

Ore composition of a bench at a mine in Windarling, West Australia.

Description

A dataset containing the geochemical composition (and some extra variables) of a grade control study of a mine bench at Windarling, West Australia.

Usage

Windarling

Format

A data.frame with 1600 rows and 16 columns:

Hole_id

ID of the observation

Sample.West

indicator for belonging to a subsample on the West sector

Sample.East

indicator for belonging to a subsample on the East sector

West

indicator for belonging to the Western wing of the bench

East

indicator for belonging to the Eastern wing of the bench

Easting

Easting (X) coordinate of the sample

Northing

Northing (Y) coordinate of the sample

Lithotype

factor, indicating material type of the sample, one of: basalt, goethite, magnetite, schist

Fe

percent of Iron in the sample

P

percent of Phosphorus in the sample

SiO2

percent of Silica in the sample

Al2O3

percent of Alumina in the sample

S

percent of Sulphur in the sample

Mn

percent of Magnanese in the sample

CL

percent of Chlorine in the sample

LOI

percent Loss-on-Ignition of the sample

License

CC BY-SA 4.0

Source

Ward (2015) Compositions, logratios and geostatistics: An application to iron ore. MSc Thesis Edith Cowan University; https://ro.ecu.edu.au/theses/1581/

Extract rows of a DataFrameStack

Description

Extract rows (and columns if you know what you are doing) from a stacked data frame

Usage

## S3 method for class 'DataFrameStack'
x[i = NULL, j = NULL, ..., drop = FALSE]

Arguments

Value

the DataFrameStack or data.frame of the selected rows and columns.

Examples

ar = array(1:30, dim = c(5,2,3), dimnames=list(obs=1:5, vars=c("A","B"), rep=1:3))
dfs = DataFrameStack(ar, stackDim="rep")
dfs[1:2,]
stackDim(dfs[1:2,])

Subsetting of gmCgram variogram structures

Description

Extraction of some variables of a gmCgram object

Usage

## S3 method for class 'gmCgram'
x[i, j = i, ...]

Arguments

Details

This function can be used to extract the model for a a subset of variables. If only i is specified, j will be taken as equal to i. If you want to select all i's for certain j's or vice versa, give i=1:dim(x$nugget)[1] and j= your desired indices, respectively j=1:dim(x$nugget)[2] and i= your desired indices; replace x by the object you are giving. If i!=j, the output will be a c("gmXCgram","gmCgram") object, otherwise it will be a regular class "gmCgram" object. If you want to extract "slots" or "elements" of the variogram, use the $-notation. If you want to extract variables of the variogram matrices, use the [-notation.

Value

a gmCgram variogram object with the desired variables only.

See Also

Other gmCgram functions: [[.gmCgram(), as.function.gmCgram(), as.gmCgram.variogramModelList(), length.gmCgram(), ndirections(), plot.gmCgram(), variogramModelPlot()

Examples

utils::data("variogramModels")
v1 = setCgram(type=vg.Gau, sill=diag(2), anisRanges = 3*diag(c(3,1)))
v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
vm = v1+v2
vm[1,1]

Subsetting of logratioVariogram objects

Description

Subsetting of logratioVariogram objects

Usage

## S3 method for class 'logratioVariogramAnisotropy'
x[i = NULL, j = NULL, ...]

Arguments

Value

the selected variograms or lags, potentially of class "logratioVariogram" if only one direction is chosen

Examples

data("jura", package="gstat")
X = jura.pred[,1:2]
Zc = compositions::acomp(jura.pred[,7:9])
vg = logratioVariogram(data=Zc, loc=X)
vg[1]
vg = logratioVariogram(data=Zc, loc=X, azimuth=c(0,90))
vg[,1]
vg[1,1]
vg[1,]

Subsetting of gmCgram variogram structures

Description

Extraction or combination of nested structures of a gmCgram object

Usage

## S3 method for class 'gmCgram'
x[[i, ...]]

Arguments

Details

This function can be used to: extract the nugget (i=0), extract some structures (i=indices of the structures, possibly including 0 for the nugget), or filter some structures out (i=negative indices of the structures to remove; nugget will always removed in this case). If you want to extract "slots" or "elements" of the variogram, use the $-notation. If you want to extract variables of the variogram matrices, use the [-notation. The contrary operation (adding structures together) is obtained by summing (+) two gmCgram objects.

Value

a gmCgram variogram object with the desired structures only.

See Also

Other gmCgram functions: [.gmCgram(), as.function.gmCgram(), as.gmCgram.variogramModelList(), length.gmCgram(), ndirections(), plot.gmCgram(), variogramModelPlot()

Examples

utils::data("variogramModels")
v1 = setCgram(type=vg.Gau, sill=diag(2), anisRanges = 3*diag(c(3,1)))
v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
vm = v1+v2
vm[[1]]

Compute accuracy and precision

Description

Computes goodness-of-fit measures (accuracy, precision and joint goodness) adapted or extended from the definition of Deutsch (1997).

Usage

accuracy(object, ...)

## S3 method for class 'data.frame'
accuracy(
  object,
  observed = object$observed,
  prob = seq(from = 0, to = 1, by = 0.05),
  method = "kriging",
  outMahalanobis = FALSE,
  ivar,
  ...
)

## S3 method for class 'DataFrameStack'
accuracy(
  object,
  observed,
  ivars = intersect(colnames(observed), dimnames(object)[[noStackDim(object)]]),
  prob = seq(from = 0, to = 1, by = 0.05),
  method = ifelse(length(ivars) == 1, "simulation", "Mahalanobis"),
  outMahalanobis = FALSE,
  ...
)

Arguments

Details

For method "kriging", object must contain columns with names including the string "pred" for predictions and "var" for the kriging variance; the observed values can also be included as an extra column with name "observed", or else additionally provided in argument observed. For method "cokriging", the columns of object must contain predictions, cokriging variances and cokriging covariances in columns including the strings "pred", "var" resp. "cov", and observed values can only be provided via observed argument. Note that these are the natural formats when using gstat::predict.gstat() and other (co)kriging functions of that package.

For univariate and multivariate cokriging results (methods "kriging" and "cokriging"), the coverage values are computed based on the Mahalanobis square error, the (square) distance between prediction and true value, using as the positive definite bilinear form of the distance the variance-covariance cokriging matrix. The rationale is that, under the assumption that the random field is Gaussian, the distribution of this Mahalanobis square error should follow a \chi^2(\nu) with degrees of freedom \nu equal to the number of variables. Having this reference distribution allows us to compute confidence intervals for that Mahalanobis square error, and then count how many of the actually observed errors are included on each one of the intervals (the coverage). For a perfect adjustment to the distribution, the plot of coverage vs. nominal confidence (see plot.accuracy) should fall on the y=x line. NOTE: the original definition of Deutsch (1997) for univariate case did not make use of the \chi^2(1) distribution, but instead derived the desired intervals (symmetric!) from the standard normal distribution appearing by normalizing the residual with the kriging variance; the result is the same.

For method "simulation" and object object is a data.frame, the variable names containing the realisations must contain the string "sim", and observed must be a vector with as many elements as rows has object. If object is a DataFrameStack(), then it is assumed that the stacking dimension is running through the realisations; the true values must still be given in observed. In both cases, the method is based on ranks: with them we can calculate which is the frequency of simulations being more extreme than the observed value. This calculation is done considering bilateral intervals around the median of (realisations, observed value) for each location separately.

Method "mahalanobis" ("Mahalanobis" also works) is the analogous for multivariate simulations. It only works for object of class DataFrameStack(), and requires the stacking dimension to run through the realisations and the other two dimensions to coincide with the dimensions of observed, i.e. giving locations by rows and variables by columns. In this case, a covariance matrix will be computed and this will be used to compute the Mahalanobis square error defined above in method "cokriging": this Mahalanobis square error will be computed for each simulation and for the true value. The simulated Mahalanobis square errors will then be used to generate the reference distribution with which to derive confidence intervals.

Finally, highly experimental "flow" method requires the input to be in the same shape as method "mahalanobis". The method is mostly the same, just that before the Mahalanobis square errors are computed a location-wise flow anamorphosis (ana()) is applied to transform the realisations (including the true value as one of them) to joint normality. The rest of the calculations are done as if with method "mahalanobis".

Value

If outMahalanobis=TRUE (the primary use), this function returns a two-column dataset of class c("accuracy", "data.frame"), which first column gives the nominal probability cutoffs used, and the second column the actual coverage of the intervals of each of these probabilities. If outMahalanobis=FALSE, the output is a vector (for prediction) or matrix (for simulation) of Mahalanobis error norms.

Methods (by class)

• data.frame: Compute accuracy and precision

• DataFrameStack: Compute accuracy and precision

References

Mueller, Selia and Tolosana-Delgado (2023) Multivariate cross-validation and measures of accuracy and precision. Mathematical Geosciences (under review).

See Also

Other accuracy functions: mean.accuracy(), plot.accuracy(), precision(), validate(), xvErrorMeasures.default(), xvErrorMeasures()

Flow anamorphosis transform Compute a transformation that gaussianizes a certain data set

Description

Flow anamorphosis transform Compute a transformation that gaussianizes a certain data set

Usage

ana(Y, sigma0 = 0.1, sigma1 = 1, steps = 30, sphere = TRUE, weights = NULL)

Arguments

Value

a function with arguments (x, inv=FALSE), where x will be the data to apply the transformation to, and inv=FALSE will indicate if the direct or the inverse transformation is desired

Author(s)

K. Gerald van den Boogaart

See Also

anaForward, anaBackward, sphTrans

Examples

library(compositions)
data("jura", package="gstat")
Y = acomp(jura.pred[,c(10,12,13)])
plot(Y)
anafun = ana(Y)
class(anafun)
z = anafun(Y)
plot(z)
y = anafun(z, inv=TRUE)
plot(data.frame(orig=Y,recalc=y))

Backward gaussian anamorphosis backward transformation to multivariate gaussian scores

Description

Backward gaussian anamorphosis backward transformation to multivariate gaussian scores

Usage

anaBackward(
  x,
  Y,
  sigma0,
  sigma1 = 1 + sigma0,
  steps = 30,
  plt = FALSE,
  sphere = TRUE,
  weights = NULL
)

Arguments

Value

a matrix with the scores back-transformed to the same scale as Y; same dimensions of x

Author(s)

K. Gerald van den Boogaart, Raimon Tolosana-Delgado

See Also

ana() for defining a function that carries over the transformation (by means of a closure), anaBackward() for the explicit back-transformation, sphTrans() for defining a function that carries over the spherification of the data

Examples

data("jura", package="gstat")
Y = jura.pred[,c(10,12,13)]
plot(compositions::acomp(Y))
Ylr = compositions::alr(Y)
Xns = matrix(rnorm(500), ncol=2)
plot(Ylr)
points(Xns, col=2, pch=4)
Xlr = anaBackward(x=Xns, Y=Ylr, sigma0=0.1)
qqplot(Xlr[,1], Ylr[,1])
qqplot(Xlr[,2], Ylr[,2])
qqplot(Xlr[,1]+Xlr[,2], Ylr[,1]+Ylr[,2])

Forward gaussian anamorphosis forward transformation to multivariate gaussian scores

Description

Forward gaussian anamorphosis forward transformation to multivariate gaussian scores

Usage

anaForward(
  x,
  Y,
  sigma0,
  sigma1 = 1 + sigma0,
  steps = 30,
  plt = FALSE,
  sphere = TRUE,
  weights = NULL
)

Arguments

Value

a matrix with the gaussian scores; same dimensions of x

Author(s)

K. Gerald van den Boogaart, Raimon Tolosana-Delgado

See Also

ana() for defining a function that carries over the transformation (by means of a closure), anaBackward() for the explicit back-transformation, sphTrans() for defining a function that carries over the spherification of the data

Examples

data("jura", package="gstat")
Y = jura.pred[,c(10,12,13)]
plot(compositions::acomp(Y))
Ylr = compositions::alr(Y)
plot(Ylr)
z = anaForward(x=Ylr, Y=Ylr, sigma0=0.1)
plot(z, asp=1)
shapiro.test(z[,1])
shapiro.test(z[,2])

Produce anisotropy scaling matrix from angle and anisotropy ratios

Description

Produce anisotropy matrix (as the transposed of the Cholesky decomposition) from angle and anisotropy ratios

Usage

anis_GSLIBpar2A(ratios, angles, inv = FALSE)

anis2D_par2A(ratio, angle, inv = FALSE)

anis3D_par2A(ratios, angles, inv = FALSE)

Arguments

Value

a 3x3 matrix of anisotropy.

If inv=TRUE (or 1) the output is a matrix A such that norm(h %*% A) allows to use isotropic variograms, being h = c(hx, hy, hz) the lag vectors.

If inv=FALSE (or 0) the output is a matrix A such that norm(h %*% solve(A)) allows to use isotropic variograms.

Other values are meaningless.

Functions

• anis2D_par2A: 2D case

• anis3D_par2A: 3D case

See Also

Other anisotropy: AnisotropyRangeMatrix(), AnisotropyScaling(), as.AnisotropyRangeMatrix(), as.AnisotropyScaling(), is.anisotropySpecification()

Examples

## ratio=0.5, azimuth 30?? (i.e. direction 60??)
A = anis2D_par2A(1, 30)
A
AAt = A %*% t(A)
 #  project the bisector 1:1 (i.e. 45??)
(k = c(1,1,0) %*% A)
atan2(k[2], k[1]) * 180/pi  # should be 15
sqrt(sum(k^2))
sqrt( c(1,1,0) %*% AAt %*% c(1,1,0) )
A = anis2D_par2A(0.5, 60)
rd = 60 * pi/180
A
A %*% t(A)
c(cos(rd), sin(rd),0) %*% A #  should be 1
c(-sin(rd), cos(rd),0) %*% A #  should be +/- sqrt(2)
c60 = cos(60*pi/180)
s60 = sin(60*pi/180)
c30 = cos(30*pi/180)
s30 = sin(30*pi/180)
#  in the new coordinates, 60cwN is (0,1,0)
R60p = anis3D_par2A(ratios=c(1,1), angles=c(60,0,0))
c(s60, c60, 0) %*% R60p
R6030 = anis3D_par2A(ratios=c(1,1), angles=c(60,30,0))
# the original X axis is positive on newX and newY, but negative on newZ
c(1,0,0) %*% R6030
# rotate first direction 60 degrees azimuth, then dip 30degrees upwards
c( c60*c30, -s60*c30, s30) %*% R6030
(Ranis = anis3D_par2A(ratios=c(0.5,0.25), angles=c(60,30,0)) )

Force a matrix to be anisotropy range matrix,

Description

Force a matrix M to be considered an anisotropy range matrix, i.e with ranges and orientations, such that u = sqrt(h' * M^{-1} * h) allows to use an isotropic variogram.

Usage

as.AnisotropyRangeMatrix(x)

## Default S3 method:
as.AnisotropyRangeMatrix(x)

## S3 method for class 'AnisotropyRangeMatrix'
as.AnisotropyRangeMatrix(x)

Arguments

Value

the same anisotropy, specified as M

Methods (by class)

• default: Default conversion to anisotropy range matrix

• AnisotropyRangeMatrix: identity conversion

See Also

Other anisotropy: AnisotropyRangeMatrix(), AnisotropyScaling(), anis_GSLIBpar2A(), as.AnisotropyScaling(), is.anisotropySpecification()

Convert to anisotropy scaling matrix

Description

Convert an anisotropy specification to a scaling matrix

Usage

as.AnisotropyScaling(x)

## S3 method for class 'AnisotropyScaling'
as.AnisotropyScaling(x)

## S3 method for class 'numeric'
as.AnisotropyScaling(x)

## S3 method for class 'AnisotropyRangeMatrix'
as.AnisotropyScaling(x)

Arguments

Details

Method as.AnisotropyScaling.numeric() expects a vector of two numbers in 2D, or a vector of 5 numbers in 3D. These are in 2D, the azimuth of maximum continuity (in degrees, clockwise from North) and the anisotropy ratio of short/long range. In 3D these are: 1,2) the azimuth and the dip of the direction of maximal continuity; 3) the angle of rotation around the axis of the first direction; 4,5) the anisotropy ratios of the ranges of the second/first and third/first directions of maximal continuity. All angles are given in degrees, all ratios must be smaller or equal to 1. This follows gstat (and hence GSlib) conventions; see gstat::vgm() for details.

Value

A matrix A such that for any lag vector h, the variogram model turns isotropic in terms of u'=h'\cdot A.

Methods (by class)

• AnisotropyScaling: identity method

• numeric: from a vector of numbers

• AnisotropyRangeMatrix: from an AnisotropicRangeMatrix

See Also

Other anisotropy: AnisotropyRangeMatrix(), AnisotropyScaling(), anis_GSLIBpar2A(), as.AnisotropyRangeMatrix(), is.anisotropySpecification()

Examples

( l = anis_GSLIBpar2A(angles=30, ratios=0.5))
( ll = unclass(l) )
as.AnisotropyScaling(ll)

Recast a model to the variogram model of package "compositions"

Description

Recast a variogram model specified in any of the models of "gstat" or "gmGeostats" in the format of compositions::CompLinModCoReg()

Usage

as.CompLinModCoReg(v, ...)

Arguments

Value

The variogram model recast to "CompLinModCoReg"

Recast compositional variogram model to format LMCAnisCompo

Description

Recast a compositional variogram model of any sort to a variation-variogram model of class "LMCAnisCompo".

Usage

## S3 method for class 'gstat'
as.LMCAnisCompo(m, ...)

## S3 method for class 'variogramModelList'
as.LMCAnisCompo(m, V = NULL, orignames = NULL, ...)

as.LMCAnisCompo(m, ...)

## S3 method for class 'LMCAnisCompo'
as.LMCAnisCompo(m, ...)

## S3 method for class 'gmCgram'
as.LMCAnisCompo(m, V = NULL, orignames = rownames(V), ...)

## S3 method for class 'CompLinModCoReg'
as.LMCAnisCompo(m, varnames, ...)

Arguments

Value

the variogram model recasted to class "LMCAnisCompo"

Methods (by class)

• gstat: Recast compositional variogram model to format LMCAnisCompo

• variogramModelList: Recast compositional variogram model to format LMCAnisCompo

• LMCAnisCompo: Recast compositional variogram model to format LMCAnisCompo

• gmCgram: Recast compositional variogram model to format LMCAnisCompo

• CompLinModCoReg: Recast a variogram model from package "compositions" to format LMCAnisCompo

Convert a stacked data frame into an array

Description

Recast a DataFrameStack() into a named array.

Usage

## S3 method for class 'DataFrameStack'
as.array(x, ...)

Arguments

Value

the data recasted as an array with appropriate names.

Examples

ll = lapply(1:3, function(i) as.data.frame(matrix(1:10+(i-1)*10, ncol=2)))
dfs = DataFrameStack(ll, stackDimName="rep")
as.array(dfs)

Express a direction as a director vector

Description

Internal methods to express a direction (in 2D or 3D) as director vector(s). These functions are not intended for direct use.

Usage

as.directorVector(x, ...)

## Default S3 method:
as.directorVector(x, ...)

## S3 method for class 'azimuth'
as.directorVector(x, D = 2, ...)

## S3 method for class 'azimuthInterval'
as.directorVector(x, D = 2, ...)

Arguments

Value

A 2- or 3- column matrix which rows represents the unit director vector of each direction specified.

Methods (by class)

• default: default method

• azimuth: method for azimuths

• azimuthInterval: method for azimuthIntervals

Convert a gmCgram object to an (evaluable) function

Description

Evaluate a gmCgram on some h values, or convert the gmCgram object into an evaluable function

Usage

## S3 method for class 'gmCgram'
as.function(x, ...)

## S3 method for class 'gmCgram'
predict(object, newdata, ...)

Arguments

Value

a function that can be evaluated normally, with an argument X and possibly another argument Y; both must have the same number of columns than the geographic dimension of the variogram (i.e. dim(x$M)[3]).

Functions

• predict.gmCgram: predict a gmCgram object on some lag vector coordinates

See Also

Other gmCgram functions: [.gmCgram(), [[.gmCgram(), as.gmCgram.variogramModelList(), length.gmCgram(), ndirections(), plot.gmCgram(), variogramModelPlot()

Examples

utils::data("variogramModels")
v1 = setCgram(type=vg.Gau, sill=diag(2)+0.5, anisRanges = 2*diag(c(3,0.5)))
v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
vm = v1+v2
vgf = as.function(vm)
(h = rbind(c(0,1), c(0,0), c(1,1)))
vgf(h)
predict(vm, h)

Convert theoretical structural functions to gmCgram format

Description

Convert covariance function or variogram models to the format gmCgram of package gmGeostats

Usage

## S3 method for class 'variogramModelList'
as.gmCgram(m, ...)

## S3 method for class 'variogramModel'
as.gmCgram(m, ...)

## S3 method for class 'LMCAnisCompo'
as.gmCgram(
  m,
  V = attr(m, "contrasts"),
  orignames = rownames(m["sill", 1]$sill),
  ...
)

as.gmCgram(m, ...)

## Default S3 method:
as.gmCgram(m, ...)

Arguments

Value

the covariance/variogram model, recasted to class gmCgram. This is a generic function. Methods exist for objects of class LMCAnisCompo (for compositional data) and variogramModelList (as provided by package gstat).

Methods (by class)

• variogramModelList: Convert theoretical structural functions to gmCgram format

• variogramModel: Convert theoretical structural functions to gmCgram format

• LMCAnisCompo: method for "LMCAnisCompo" variogram model objects

• default: Convert theoretical structural functions to gmCgram format

See Also

Other gmCgram functions: [.gmCgram(), [[.gmCgram(), as.function.gmCgram(), length.gmCgram(), ndirections(), plot.gmCgram(), variogramModelPlot()

Convert empirical structural function to gmEVario format

Description

Convert empirical covariance functions or variograms to the format gmEVario of package gmGeostats

Usage

## S3 method for class 'gstatVariogram'
as.gmEVario(vgemp, ...)

## S3 method for class 'logratioVariogram'
as.gmEVario(vgemp, ...)

## S3 method for class 'logratioVariogramAnisotropy'
as.gmEVario(vgemp, ...)

as.gmEVario(vgemp, ...)

## Default S3 method:
as.gmEVario(vgemp, ...)

Arguments

Value

the empirical covariance function or variogram, recasted to class gmEVario. This is a generic function. Methods exist for objects of class logratioVariogramlogratioVariogramAnisotropy (for compositional data) and gstatVariogram (from package gstat).

Methods (by class)

• gstatVariogram: gstatVariogram method not yet available

• logratioVariogram: logratioVariogram method not yet available

• logratioVariogramAnisotropy: logratioVariogramAnisotropy method not yet available

• default: default method

See Also

Other gmEVario functions: gsi.EVario2D(), gsi.EVario3D(), ndirections(), plot.gmEVario(), variogramModelPlot()

Recast spatial object to gmSpatialModel format

Description

Recast a spatial data object model to format gmSpatialModel

Usage

as.gmSpatialModel(object, ...)

## Default S3 method:
as.gmSpatialModel(object, ...)

## S3 method for class 'gstat'
as.gmSpatialModel(object, V = NULL, ...)

Arguments

Value

The same spatial object re-structured as a "gmSpatialModel", see gmSpatialModel

Methods (by class)

• default: Recast spatial object to gmSpatialModel format

• gstat: Recast spatial object to gmSpatialModel format

See Also

Other gmSpatialModel: Predict(), gmSpatialModel-class, make.gmCompositionalGaussianSpatialModel(), make.gmCompositionalMPSSpatialModel(), make.gmMultivariateGaussianSpatialModel()

Convert a regionalized data container to gstat

Description

Convert a regionalized data container to a "gstat" model object

Usage

as.gstat(object, ...)

## Default S3 method:
as.gstat(object, ...)

Arguments

Value

A regionalized data container of class "gstat", eventually with variogram model included. See gstat::gstat() for more info.

Functions

• as.gstat.default: default does nothing

Examples

data("jura", package = "gstat")
X = jura.pred[,1:2]
Zc = jura.pred[,7:13]
gg = make.gmCompositionalGaussianSpatialModel(Zc, X, V="alr", formula = ~1)
as.gstat(gg)

Represent an empirical variogram in "gstatVariogram" format

Description

Represent an empirical variogram in "gstatVariogram" format, from package "gstat"; see gstat::variogram() for details.

Usage

as.gstatVariogram(vgemp, ...)

## Default S3 method:
as.gstatVariogram(vgemp, ...)

## S3 method for class 'gmEVario'
as.gstatVariogram(vgemp, ...)

## S3 method for class 'logratioVariogram'
as.gstatVariogram(
  vgemp,
  V = NULL,
  dir.hor = 0,
  dir.ver = 0,
  prefix = NULL,
  ...
)

## S3 method for class 'logratioVariogramAnisotropy'
as.gstatVariogram(vgemp, V = NULL, ...)

Arguments

Value

The function returns an object of class "gstatVariogram" containing the empirical variogram provided. See gstat::variogram() for details.

Methods (by class)

• default: Represent an empirical variogram in "gstatVariogram" format

• gmEVario: Represent an empirical variogram in "gstatVariogram" format (not yet available)

• logratioVariogram: Represent an empirical variogram in "gstatVariogram" format

• logratioVariogramAnisotropy: Represent an empirical variogram in "gstatVariogram" format

Examples

data("jura", package = "gstat")
X = jura.pred[,1:2]
Zc = compositions::acomp(jura.pred[,7:13])
lrvg = gmGeostats::logratioVariogram(data=Zc, loc=X)
as.gstatVariogram(lrvg, V="alr")

Convert a stacked data frame into a list of data.frames

Description

Recast a DataFrameStack() into a list of data.frames

Usage

## S3 method for class 'DataFrameStack'
as.list(x, ...)

Arguments

Value

the data recasted as list of data.frames

Examples

ar = array(1:30, dim = c(5,2,3), dimnames=list(obs=1:5, vars=c("A","B"), rep=1:3))
dfs = DataFrameStack(ar, stackDim="rep")
as.list(dfs)

Recast empirical variogram to format logratioVariogram

Description

Recast an empirical compositional variogram of any sort to a variation-variogram of class "logratioVariogram".

Usage

## S3 method for class 'gstatVariogram'
as.logratioVariogram(
  vgemp,
  V = NULL,
  tol = 1e-12,
  orignames = NULL,
  symmetrize = FALSE,
  ...
)

as.logratioVariogram(vgemp, ...)

## S3 method for class 'logratioVariogram'
as.logratioVariogram(vgemp, ...)

## S3 method for class 'gmEVario'
as.logratioVariogram(vgemp, ...)

Arguments

Value

the same model in the new format.

Methods (by class)

• gstatVariogram: method for gstatVariogram objects

• logratioVariogram: identity method

• gmEVario: method for gmEVario objects (not yet available)

Convert empirical variogram to "logratioVariogramAnisotropy"

Description

Convert an empirical variogram from any format to class "logratioVariogramAnisotropy"

Usage

as.logratioVariogramAnisotropy(vgemp, ...)

## Default S3 method:
as.logratioVariogramAnisotropy(vgemp, ...)

## S3 method for class 'logratioVariogram'
as.logratioVariogramAnisotropy(vgemp, ...)

## S3 method for class 'logratioVariogramAnisotropy'
as.logratioVariogramAnisotropy(vgemp, ...)

Arguments

Value

The empirical variogram as a "logratioVariogramAnisotropy" object

Methods (by class)

• default: default method, making use of as.logratioVariogram()

• logratioVariogram: method for "logratioVariogram" class

• logratioVariogramAnisotropy: identity transformation

Convert an LMC variogram model to gstat format

Description

Convert a linear model of coregionalisation to the format of package gstat. See gstat::vgm() for details.

Usage

as.variogramModel(m, ...)

## Default S3 method:
as.variogramModel(m, ...)

## S3 method for class 'gmCgram'
as.variogramModel(m, ...)

## S3 method for class 'LMCAnisCompo'
as.variogramModel(m, V = NULL, prefix = NULL, ensurePSD = TRUE, ...)

## S3 method for class 'CompLinModCoReg'
as.variogramModel(m, V = "alr", prefix = NULL, ensurePSD = TRUE, ...)

Arguments

Value

The LMC model specified in the format of package gstat, i.e. as the result of using gstat::vgm()

Methods (by class)

• default: Convert an LMC variogram model to gstat format

• gmCgram: Convert an LMC variogram model to gstat format

• LMCAnisCompo: Convert an LMC variogram model to gstat format

• CompLinModCoReg: Convert an LMC variogram model to gstat format

Examples

data("jura", package = "gstat")
X = jura.pred[,1:2]
Zc = compositions::acomp(jura.pred[,7:13])
lrmd = compositions::CompLinModCoReg(formula=~nugget()+sph(1.5), comp=Zc)
as.variogramModel(lrmd, V="alr")

Colored biplot for gemeralised diagonalisations Colored biplot method for objects of class genDiag

Description

Colored biplot for gemeralised diagonalisations Colored biplot method for objects of class genDiag

Usage

## S3 method for class 'genDiag'
coloredBiplot(x, choices = 1:2, scale = 0, pc.biplot = FALSE, ...)

Arguments

Value

nothing. Function is called exclusively to produce the plot

References

Mueller, Tolosana-Delgado, Grunsky and McKinley (2020) Biplots for compositional data derived from generalised joint diagonalization methods. Applied Computational Geosciences 8:100044 doi: 10.1016/j.acags.2020.100044

See Also

Other generalised Diagonalisations: Maf(), predict.genDiag()

Examples

data("jura", package="gstat")
juracomp = compositions::acomp(jura.pred[, -(1:6)]) 
lrvg = logratioVariogram(data=juracomp, loc=jura.pred[,1:2])
mf = Maf(juracomp, vg=lrvg)
mf
compositions::coloredBiplot(mf, xlabs.col=as.integer(jura.pred$Rock)+2)

Constructs a mask for a grid

Description

Constructs a mask for a grid

Usage

constructMask(grid, method = "maxdist", maxval = NULL, x = NULL)

Arguments

Details

Method 'maxdist' defines the mask as all points within a maximum distance (must be given in maxval) from the reference data (given in x: this is expected to be the original complete data, with coordinates and variables). For method 'sillprop' the mask is defined by those points which total kriging variance is below a fixed proportion (given in maxval, default=0.99) of the total variogram model sill (variogram model given in x, of class "variogramModelList"). In this method, the argument grid is expected to be the output of a cokriging analysis. Finally, method 'point2poly' created the mask by taking the points internal to a "SpatialPolygon" object (given in x).

Value

a logical vector with as many elements as points in the grid, with TRUE for those points within the mask, and FALSE for those outside the mask.

See Also

Other masking functions: getMask(), print.mask(), setMask(), unmask()

Examples

## with data.frame
x = 1:23
y = 1:29
xy = expand.grid(x=x, y=y)
xyz.df = data.frame(xy, z = rnorm(29*23)*ifelse(abs(xy$x-xy$y)<3, 1, NA)+(xy$x+xy$y)/2)
mask.df = constructMask(grid = xy, method = "maxdist", maxval = 3, x=xyz.df)
image(mask.df)
oldpar = par(mfrow = c(1,1))
mask.df
xyz.df.masked = setMask(xyz.df, mask.df)
dim(xyz.df.masked)
summary(xyz.df.masked)
xyz.df.unmasked = unmask(xyz.df.masked)
dim(xyz.df.unmasked)
length(x)*length(y)
summary(xyz.df.unmasked)
## with SpatialGrid
library(sp)
library(magrittr)
xy.sp = sp::SpatialPoints(coords = xy)
meandiff = function(x) mean(diff(x))
xy.gt = GridTopology(c(min(x),min(y)), c(meandiff(x), meandiff(y)), c(length(x),length(y)))
aux = sp::SpatialPixelsDataFrame(grid = xy.gt, data=xyz.df, points = xy.sp)
xyz.sgdf = as(aux, "SpatialGridDataFrame")
image_cokriged(xyz.sgdf, ivar="z")
## reorder the data in the grid and plot again
par(mfrow=c(1,1))
ms = function(x) sortDataInGrid(x, grid=xy.gt)
mask.gt = constructMask(grid = xy.gt, method = "maxdist", maxval = 3, x=xyz.sgdf)
image(x,y,matrix(ms(xyz.sgdf@data$z), nrow=23, ncol=29))  
image(x,y,matrix(ms(mask.gt), nrow=23, ncol=29))  
image(mask.gt)
## work with the mask and plot again
par(mfrow=c(1,1))
xyz.sgdf.masked = setMask(x = xyz.sgdf, mask = mask.gt)
getMask(xyz.sgdf.masked)
image(x,y,matrix(ms(xyz.sgdf@data$z), nrow=23, ncol=29))  
points(xyz.sgdf.masked@coords)
par(oldpar)

Return the dimnames of a DataFrameStack

Description

Return the dimnames of a DataFrameStack(), i.e. the three dimensions

Usage

## S3 method for class 'DataFrameStack'
dimnames(x)

## S4 method for signature 'Spatial'
dimnames(x)

Arguments

Value

a list (possibly named) with the element names of each of the three dimensions

Functions

• dimnames,Spatial-method: dimnames method for all Spatial*DataFrame objects of package sp which data slot contains a DataFrameStack()

Fit an LMC to an empirical variogram

Description

Fit a linear model of coregionalisation to an empirical variogram

Usage

fit_lmc(v, ...)

## S3 method for class 'gstatVariogram'
fit_lmc(
  v,
  g,
  model,
  fit.ranges = FALSE,
  fit.lmc = !fit.ranges,
  correct.diagonal = 1,
  ...
)

## Default S3 method:
fit_lmc(v, g, model, ...)

## S3 method for class 'logratioVariogram'
fit_lmc(v, g, model, ...)

## S3 method for class 'logratioVariogramAnisotropy'
fit_lmc(v, g, model, ...)

Arguments

Value

Method fit_lmc.gstatVariogram is a wrapper around gstat::fit.lmc(), that calls this function and gives the resulting model its appropriate class (c("variogramModelList", "list")). Method fit_lmc.default returns the fitted lmc (this function currently uses gstat as a calculation machine, but this behavior can change in the future)

Methods (by class)

• gstatVariogram: wrapper around gstat::fit.lmc method

• default: flexible wrapper method for any class for which methods for as.gstatVariogram(), as.gstat() and as.variogramModel() exist. In the future there may be direct specialised implementations not depending on package gstat.

• logratioVariogram: method for logratioVariogram wrapping compositions::fit.lmc. In the future there may be direct specialised implementations, including anisotropy (not yet possible).

• logratioVariogramAnisotropy: method for logratioVariogram with anisotropry

Examples

data("jura", package = "gstat")
X = jura.pred[,1:2]
Zc = jura.pred[,7:13]
gg = make.gmCompositionalGaussianSpatialModel(Zc, X, V="alr", formula = ~1)
vg = variogram(gg)
md = gstat::vgm(model="Sph", psill=1, nugget=1, range=1.5)
gg = fit_lmc(v=vg, g=gg, model=md)
variogramModelPlot(vg, model=gg)

Get the mask info out of a spatial data object

Description

Retrieve the mask information from an object (if present). See constructMask() for examples.

Usage

getMask(x)

## Default S3 method:
getMask(x)

## S3 method for class 'SpatialPixelsDataFrame'
getMask(x)

## S3 method for class 'SpatialPixels'
getMask(x)

## S3 method for class 'SpatialPointsDataFrame'
getMask(x)

Arguments

Value

The retrieved mask information from x, an object of class "mask"

Methods (by class)

• default: Get the mask info out of a spatial data object

• SpatialPixelsDataFrame: Get the mask info out of a SpatialPixelsDataFrame data object

• SpatialPixels: Get the mask info out of a SpatialPixels object

• SpatialPointsDataFrame: Get the mask info out of a SpatialPointsDataFrame data object

See Also

Other masking functions: constructMask(), print.mask(), setMask(), unmask()

Set or get the i-th data frame of a data.frame stack

Description

Set or get one element of the DataFrameStack()

Usage

getStackElement(x, i, ...)

## Default S3 method:
getStackElement(x, i, ...)

## S3 method for class 'list'
getStackElement(x, i, ...)

## S3 method for class 'DataFrameStack'
getStackElement(x, i, MARGIN = stackDim(x), ...)

## Default S3 method:
setStackElement(x, i, value, ...)

## S3 method for class 'data.frame'
setStackElement(x, i, value, ...)

## S3 method for class 'list'
setStackElement(x, i, value, ...)

## S3 method for class 'DataFrameStack'
setStackElement(x, i, value, MARGIN = stackDim(x), ...)

Arguments

Value

For the getters, the result is the data.frame of the stack asked for. For the setters the result is the original DataFrameStack with the corresponding element replaced. Spatial methods return the corresponding spatial object, ie. the spatial information of the stack is transferred to the extracted element.

Methods (by class)

• default: Set or get one element of the DataFrameStack()

• list: Set or get one element of the DataFrameStack()

• DataFrameStack: Set or get one element of the DataFrameStack()

• default: Set or get one element of the DataFrameStack()

• data.frame: Set one element of the DataFrameStack() in data.frame form

• list: Set get one element of a DataFrameStack() in list form

• DataFrameStack: Set one element of the DataFrameStack()

Examples

ar = array(1:30, dim = c(5,2,3), dimnames=list(obs=1:5, vars=c("A","B"), rep=1:3))
dfs = DataFrameStack(ar, stackDim="rep")
dfs
stackDim(dfs)
getStackElement(dfs, 1)

Download the Tellus survey data set (NI)

Description

Download the A-soil geochemistry data of the Tellus survey (Northern Ireland), and if desired produce a training image of the geochemical subcomposition Mg-Al-Ca-Fe-Rest.

Usage

getTellus(wd = ".", destfile = "TellusASoil.RData", TI = FALSE, cleanup = TRUE)

Arguments

Details

The function is provided due to conflicting licenses. You download the data from the server https://www.bgs.ac.uk/gsni/tellus/index.html at your own accord. Actually a visit to the project data webpage is highly recommended, in particular for learning about the QA/QC process of the data acquisition, other data sources available and how to use this wealth of data.

Value

TRUE if everything went OK. The function is actually called for the side effect of downloading data from the Tellus survey project website, so be careful to call it specifying the right directory at wd.

Geochemical soil A survey

The function will download an excel file to your working directory, do some manipulations, save the resulting data.frame in a "TellusASoil.RData" file (or any other name you provide on the destfile argument) and eventually, remove the excel file (if you left cleanup=TRUE). This data set has 6862 obervations and 58 variables:

Sample

Sample ID

EASTING

X coordinate of each point

NORTHING

Y coordinate of each point

Flag

a flag marking some data as a study subset, fixed for reproducibility

Ag

Silver concentration in ppm

Cd

Cadmium concentration in ppm

In

Indium concentration in ppm

Sn

Tin concentration in ppm

Sb

Antimony concentration in ppm

Te

Tellurium concentration in ppm

I

Indim concentration in ppm

Cs

Caesium concentration in ppm

Ba

Barium concentration in ppm

La

Lanthanum concentration in ppm

Ce

Cerium concentration in ppm

Na2O

Sodium oxide in percent

MgO

Magnesium oxide in percent

Al2O3

Aluminium oxide in percent

SiO2

Silicium oxide in percent

P2O5

Phosphorous(V) oxide in percent

SO3

Sulphur(III) oxide in percent

K2O

Potasium oxide in percent

CaO

Calcium oxide in percent

TiO2

Titanium oxide in percent

MnO

Manganese oxide in percent

Fe2O3

Total Iron(III) oxide in percent

Cl

Chlorine concentration in ppm

Sc

Scandium concentration in ppm

V

Vanadium concentration in ppm

Cr

Chromium concentration in ppm

Co

Cobalt concentration in ppm

Ni

Nickel concentration in ppm

Cu

Copper concentration in ppm

Zn

Zinc concentration in ppm

Ga

Callium concentration in ppm

Ge

Germanium concentration in ppm

As

Arsenic concentration in ppm

Se

Selenium concentration in ppm

Br

Bromine concentration in ppm

Rb

Rubidium concentration in ppm

Sr

Strontium concentration in ppm

Y

Yttrium concentration in ppm

Zr

Zirconium concentration in ppm

Nb

Niobium concentration in ppm

Mo

Molybdenum concentration in ppm

Nd

Neodymium concentration in ppm

Sm

Samarium concentration in ppm

Yb

Ytterbium concentration in ppm

Hf

Hafnium concentration in ppm

Ta

Tantalum concentration in ppm

W

Tungsten concentration in ppm

Tl

Thallium concentration in ppm

Pb

Lead concentration in ppm

Bi

Bismuth concentration in ppm

Th

Thorium concentration in ppm

U

Uranium concentration in ppm

pH

soil acidity

LOI

loss on ignition in percent

Training image

Additionally, if you give TI!=FALSE, the function will produce an additional file "Tellus_TI.RData" (if TI=TRUE, or any other filename that you specify on the argument TI) with a data.frame with 13287 rows and 8 columns:

EASTING

X coordinate of each cell point

NORTHING

Y coordinate of each cell point

MgO

Magnesia proportion

Al2O3

Alumina proportion

CaO

Calcium oxide proportion

Fe2O3

Iron oxide proportion

Rest

Residual complementing the sum to 1

Mask

Indicator of grid point outside the boundary of the country. NOTE: to use it with setMask() you will need to invert it using !as.logical(TellusTI$Mask)

This is a migrated version of the data set to a regular grid, ideal for illustrating and testing multiple point geostatistical algorithms.

References

https://www.bgs.ac.uk/gsni/tellus/index.html

Examples

## Not run: 
 getwd()
 getTellus(TI=TRUE, cleanup=TRUE)
 dir(pattern="Tellus")

## End(Not run)

Apply Functions Over Array or DataFrameStack Margins

Description

Returns a vector or array or list of values obtained by applying a function to the margins of an array or matrix. Method gmApply.default() is just a wrapper on base::apply(). Method gmApply() reimplements the functionality with future access to parallel computing and appropriate default values for the MARGIN. ALWAYS use named arguments here!

Usage

gmApply(X, ...)

## Default S3 method:
gmApply(X, MARGIN, FUN, ...)

## S3 method for class 'DataFrameStack'
gmApply(X, MARGIN = stackDim(X), FUN, ..., .parallel = FALSE)

Arguments

Value

In principle, if MARGIN==stackDim(X) (the default), the oputput is a list with the result of using FUN on each element of the stack. If FUN returns a matrix or a data.frame assimilable to one element of the stack, a transformation of this output to a DataFrameStack is attempted.

For X non-DataFrameStack or MARGIN!=stackDim(X) see base::apply().

Methods (by class)

• default: wrapper around base::apply()

• DataFrameStack: Apply Functions Over DataFrameStack Margins

Examples

dm = list(point=1:100, var=LETTERS[1:2], rep=paste("r",1:5, sep=""))
ar = array(rnorm(1000), dim=c(100,2,5), dimnames = dm)
dfs = DataFrameStack(ar, stackDim="rep")
gmApply(dfs, FUN=colMeans)
rs = gmApply(dfs, FUN=function(x) x+1)
class(rs)
getStackElement(rs,1)
getStackElement(dfs,1)

parameters for Spatial Gaussian methods of any kind

Description

abstract class, containing any parameter specification for a spatial algorithm for interpolation, simulation or validation making use of Gaussian assumptions

parameters for Gaussian Simulation methods

Description

abstract class, containing any parameter specification of a spatial simulation algorithm exploiting a Gaussian two-point model structure

parameters for Multiple-Point Statistics methods

Description

abstract class, containing any parameter specification of a spatial multipoint algorithm

Neighbourhood description

Description

abstract class, containing any specification of a spatial neighbourhood

Parameter specification for a spatial simulation algorithm

Description

abstract class, containing any parameter specification for a spatial simulation algorithm

General description of a spatial data container

Description

abstract class, containing any specification of a spatial data container

Details

setClassUnion(name="gmTrainingImage", members=c("SpatialGridDataFrame", "SpatialPixelsDataFrame") )

Parameter specification for any spatial method

Description

abstract class, containing any parameter specification for any spatial method. Members of this class are gmNeighbourhoodSpecification gmMPSParameters and gmValidationStrategy.

Conditional spatial model data container

Description

This class is devised to contain a conditional spatial model, with: some conditioning data (a sp::SpatialPointsDataFrame()), an unconditional geospatial model (a structure with e.g. a training image; or the information defining a Gaussian random field); and eventually some extra method parameters. The class extends sp::SpatialPointsDataFrame() and has therefore its slots, plus model (for the unconditional model) and parameters (for the extra method information)

Usage

## S4 method for signature 'gmSpatialModel'
variogram(object, methodPars = NULL, ...)

## S4 method for signature 'gmSpatialModel'
logratioVariogram(data, ..., azimuth = 0, azimuth.tol = 180/length(azimuth))

## S4 method for signature 'gmSpatialModel'
as.gstat(object, ...)

Arguments

Value

You will seldom create the spatial model directly. Use instead the creators ⁠make.gm*⁠ linked below

Methods (by generic)

• variogram: Compute a variogram, see variogram_gmSpatialModel() for details

• logratioVariogram: S4 wrapper method around logratioVariogram() for gmSpatialModel objects

• as.gstat: convert from gmSpatialModel to gstat; see as.gstat() for details

Slots

data

a data.frame (or class extending it) containing the conditional data

coords

a matrix or dataframe of 2-3 columns containing the sampling locations of the conditional data

coords.nrs

see sp::SpatialPointsDataFrame()

bbox

see sp::SpatialPointsDataFrame()

proj4string

see sp::SpatialPointsDataFrame()

model

gmUnconditionalSpatialModel. Some unconditional geospatial model. It can be NULL.

parameters

gmSpatialMethodParameters. Some method parameters. It can be NULL

See Also

Other gmSpatialModel: Predict(), as.gmSpatialModel(), make.gmCompositionalGaussianSpatialModel(), make.gmCompositionalMPSSpatialModel(), make.gmMultivariateGaussianSpatialModel()

Examples

data("jura", package="gstat")
library(sp)
X = jura.pred[,1:2]
Zc = jura.pred[,7:13]
spdf = sp::SpatialPointsDataFrame(coords=X, data=Zc)
new("gmSpatialModel", spdf)
make.gmCompositionalGaussianSpatialModel(data=Zc, coords=X, V="alr")

MPS training image class

Description

abstract class, containing any specification of a multiple-point training image. This must be analogous to sp::SpatialGridDataFrame(), and implement a method for sp::getGridTopology() and to coerce it to sp::SpatialGridDataFrame.

General description of a spatial model

Description

abstract class, containing any specification of an unconditional spatial model

Validation strategy description

Description

abstract class, containing any specification of a validation strategy for spatial models

Cokriging of all sorts, internal function

Description

internal function to compute cokriging (simple, ordinary, universal, with trend)

Usage

gsi.Cokriging(
  Xout,
  Xin,
  Zin,
  vgram,
  Fin = rep(1, nrow(Xin)),
  Fout = rep(1, nrow(Xout)),
  krigVar = FALSE,
  tol = 1e-15
)

Arguments

Value

A (M,D)-matrix of predictions, eventually with an attribute "krigVar" containing the output defined by argument krigVar

Internal function, conditional turning bands realisations

Description

Internal function to compute conditional turning bands simulations

Usage

gsi.CondTurningBands(
  Xout,
  Xin,
  Zin,
  vgram,
  nbands = 400,
  tol = 1e-15,
  nsim = NULL
)

Arguments

Value

an array with (npoint, nvar, nsim)-elements, being npoint=nrow(X) and nvar = nr of variables in vgram

Workhorse function for direct sampling

Description

This function implements in R the direct sampling algorithm

Usage

gsi.DS(
  n,
  f,
  t,
  n_realiz,
  dim_TI,
  dim_SimGrid,
  TI_input,
  SimGrid_input,
  ivars_TI = 3:ncol(TI_input),
  SimGrid_mask = ncol(SimGrid_input),
  invertMask = TRUE
)

Arguments

Value

A sp::SpatialPixelsDataFrame() or sp::SpatialGridDataFrame(), depending on whether the whole grid is simulated. The '@data' slot of these objects contains a DataFrameStack() with the stacking dimension running through the realisations. It is safer to use this functionality through the interface make.gmCompositionalMPSSpatialModel(), then request a direct simulation with DSpars() and finally run it with predict_gmSpatialModel.

Author(s)

Hassan Talebi (copyright holder), Raimon Tolosana-Delgado

Examples

## training image:
x = 1:10
y = 1:7
xy_TI = expand.grid(x=x, y=y)
TI_input = cbind(xy_TI, t(apply(xy_TI, 1, function(x) c(sum(x), abs(x[2]-x[1]))+rnorm(2, sd=0.01))))
colnames(TI_input) = c("x", "y", "V1", "V2")
o1 = image_cokriged(TI_input, ivar="V1")
o2 = image_cokriged(TI_input, ivar="V2")
## simulation grid:
SimGrid = TI_input
SimGrid$mask = with(SimGrid, x==1 | x==10 | y==1 | y==7)
tk = SimGrid$mask
tk[sample(70, 50)] = TRUE 
SimGrid[tk,3:4]=NA
image_cokriged(SimGrid, ivar="V1", breaks=o1$breaks, col=o1$col)
image_cokriged(SimGrid, ivar="V2", breaks=o2$breaks, col=o2$col)
image_cokriged(SimGrid, ivar="mask", breaks=c(-0.0001, 0.5, 1.001))
## Not run: 
res = gsi.DS(n=5, f=0.75, t=0.05, n_realiz=2, dim_TI=c(10,7),  dim_SimGrid=c(10,7), 
       TI_input=as.matrix(TI_input), SimGrid_input=as.matrix(SimGrid), 
       ivars_TI = c("V1", "V2"), SimGrid_mask="mask", invertMask=TRUE)
image_cokriged(cbind(xy_TI, getStackElement(res,1)), ivar="V1", breaks=o1$breaks, col=o1$col)
image_cokriged(cbind(xy_TI, getStackElement(res,2)), ivar="V1", breaks=o1$breaks, col=o1$col)
image_cokriged(cbind(xy_TI, getStackElement(res,1)), ivar="V2", breaks=o2$breaks, col=o2$col)
image_cokriged(cbind(xy_TI, getStackElement(res,2)), ivar="V2", breaks=o2$breaks, col=o2$col)

## End(Not run)

Empirical variogram or covariance function in 2D

Description

compute the empirical variogram or covariance function in a 2D case study

Usage

gsi.EVario2D(
  X,
  Z,
  Ff = rep(1, nrow(X)),
  maxdist = max(dist(X[sample(nrow(X), min(nrow(X), 1000)), ]))/2,
  lagNr = 15,
  lags = seq(from = 0, to = maxdist, length.out = lagNr + 1),
  azimuthNr = 4,
  azimuths = seq(from = 0, to = 180, length.out = azimuthNr + 1)[1:azimuthNr],
  maxbreadth = Inf,
  minpairs = 10,
  cov = FALSE
)

Arguments

Value

An empirical variogram for the provided data. NOTE: avoid using directly gsi.* functions! They represent either internal functions, or preliminary, not fully-tested functions. Use variogram instead.

See Also

Other gmEVario functions: as.gmEVario.gstatVariogram(), gsi.EVario3D(), ndirections(), plot.gmEVario(), variogramModelPlot()

Examples

library(gstat)
data("jura", package = "gstat")
X = as.matrix(jura.pred[,1:2])
Z = as.matrix(jura.pred[,c("Zn","Cd","Pb")])
vge = gsi.EVario2D(X,Z)
dim(vge)
dimnames(vge)
class(vge["gamma",1])
dim(vge["gamma",1][[1]])
vge["npairs",1]
vge["lags",1]

Empirical variogram or covariance function in 3D

Description

compute the empirical variogram or covariance function in a 3D case study

Usage

gsi.EVario3D(
  X,
  Z,
  Ff = rep(1, nrow(X)),
  maxdist = max(dist(X[sample(nrow(X), min(nrow(X), 1000)), ]))/2,
  lagNr = 15,
  lags = seq(from = 0, to = maxdist, length.out = lagNr + 1),
  dirvecs = t(c(1, 0, 0)),
  angtol = 90,
  maxbreadth = Inf,
  minpairs = 10,
  cov = FALSE
)

Arguments

Value

An empirical variogram for the provided data. NOTE: avoid using directly gsi.* functions! They represent either internal functions, or preliminary, not fully-tested functions. Use variogram instead.

See Also

Other gmEVario functions: as.gmEVario.gstatVariogram(), gsi.EVario2D(), ndirections(), plot.gmEVario(), variogramModelPlot()

Internal function, unconditional turning bands realisations

Description

Internal function to compute unconditional turning bands simulations

Usage

gsi.TurningBands(X, vgram, nBands, nsim = NULL)

Arguments

Value

an array with (npoint, nvar, nsim)-elements, being npoint=nrow(X) and nvar = nr of variables in vgram

Compute covariance matrix oout of locations

Description

internal function to compute the variogram model for all combinations of two

Usage

gsi.calcCgram(X, Y, vgram, ijEqual = FALSE)

Arguments

Value

a matrix of elements of the covariance, in NxM blocks of DxD elements (being N and M resp. the nr. of rows of X and Y, and D the nr of variables)

Reorganisation of cokriged compositions

Description

Produce compositional predictions out of a gstat::gstat() prediction

Usage

gsi.gstatCokriging2compo(COKresult, ...)

## Default S3 method:
gsi.gstatCokriging2compo(COKresult, ...)

## S3 method for class 'data.frame'
gsi.gstatCokriging2compo(
  COKresult,
  V = NULL,
  orignames = NULL,
  tol = 1e-12,
  nscore = FALSE,
  gg = NULL,
  ...
)

## Default S3 method:
gsi.gstatCokriging2rmult(COKresult, ...)

## S3 method for class 'data.frame'
gsi.gstatCokriging2rmult(COKresult, nscore = FALSE, gg = NULL, ...)

Arguments

Value

an (N,D)-object of class c("spatialGridAcomp","acomp") with the predictions, together with an extra attribute "krigVar" containing the cokriging covariance matrices in an (N, D, D)-array; here N=number of interpolated locations, D=number of original components of the composition

Methods (by class)

• default: Reorganisation of cokriged compositions

• data.frame: Reorganisation of cokriged compositions

• default: Reorganisation of cokriged multivariate data

• data.frame: Reorganisation of cokriged multivariate data

See Also

image_cokriged.spatialGridRmult() for an example

extract information about the original data, if available

Description

originally implemented in package:compositions

Usage

gsi.orig(x, y = NULL)

Arguments

Value

The original untransformed data (with a compositional class), resp the TRANSPOSED INVERSE matrix, i.e. the one to recover the original components.

Create a matrix of logcontrasts and name prefix

Description

Given a representation specification for compositions, this function creates the matrix of logcontrasts and provides a suitable prefix name for naming variables.

Usage

gsi.produceV(
  V = NULL,
  D = nrow(V),
  orignames = rownames(V),
  giveInv = FALSE,
  prefix = NULL
)

Arguments

Value

A list with at least two elements

1. V containing the final matrix of logcontrasts

2. prefix containing the final prefix for names of transformed variables

3. W eventually, the (transposed, generalised) inverse of V, if giveInv=TRUE

Examples

gsi.produceV("alr", D=3)
gsi.produceV("ilr", D=3, orignames = c("Ca", "K", "Na"))
gsi.produceV("alr", D=3, orignames = c("Ca", "K", "Na"), giveInv = TRUE)

Check presence of missings check presence of missings in a data.frame

Description

Check presence of missings check presence of missings in a data.frame

Usage

## S3 method for class 'data.frame'
has.missings(x, mc = NULL, ...)

Arguments

Value

A single true/false response about the presence of any missing value on the whole data set

Examples

library(compositions)
data(Windarling)
has.missings(Windarling)
head(Windarling)
Windarling[1,1] = NA
head(Windarling)
has.missings(Windarling)

Plot variogram maps for anisotropic logratio variograms

Description

Image method to obtain variogram maps for anisotropic logratio variograms

Usage

## S3 method for class 'logratioVariogramAnisotropy'
image(
  x,
  jointColor = FALSE,
  breaks = NULL,
  probs = seq(0, 1, 0.1),
  col = spectralcolors,
  ...
)

Arguments

Value

This function is called for its effect of producing a figure. Additionally, the graphical parameters active prior to calling this function are returned invisibly.

Examples

data("jura", package="gstat")
X = jura.pred[,1:2]
Zc = compositions::acomp(jura.pred[,7:9])
vg = logratioVariogram(data=Zc, loc=X, azimuth=c(0:18)*10, 
                           azimuth.tol=22.5)
image(vg)
image(vg, jointColor=TRUE)

Image method for mask objects

Description

Plot a mask of a 2D grid; see constructMask() for an example

Usage

## S3 method for class 'mask'
image(x, col = c(NA, 2), ...)

Arguments

Value

nothing

Plot an image of gridded data

Description

Plot an image of one variable (possibly, one realisation) of output of cokriging or cosimulation functions.

Usage

image_cokriged(x, ...)

## Default S3 method:
image_cokriged(
  x,
  ivar = 3,
  breaks = quantile(as.data.frame(x)[, ivar], probs = c(0:10)/10, na.rm = TRUE),
  col = spectralcolors(length(breaks) - 1),
  legendPropSpace = 0.2,
  legendPos = "top",
  main = ifelse(is.character(ivar), ivar, colnames(x)[ivar]),
  ...
)

## S3 method for class 'spatialGridRmult'
image_cokriged(
  x,
  ivar = 1,
  isim = NULL,
  breaks = 10,
  mask = attr(x, "mask"),
  col = spectralcolors(length(breaks) - 1),
  legendPropSpace = 0.2,
  legendPos = "top",
  main = ifelse(is.character(ivar), ivar, dimnames(x)[[length(dimnames(x))]][ivar]),
  ...
)

## S3 method for class 'spatialGridAcomp'
image_cokriged(
  x,
  ivar = 1,
  isim = NULL,
  breaks = 10,
  mask = attr(x, "mask"),
  col = spectralcolors(length(breaks) - 1),
  legendPropSpace = 0.2,
  legendPos = "top",
  main = ifelse(is.character(ivar), ivar, dimnames(x)[[length(dimnames(x))]][ivar]),
  ...
)

Arguments

Value

Invisibly, a list with elements breaks and col containing the breaks and hexadecimal colors finally used for the z-values of the image. Particularly useful for plotting other plotting elements on the same color scale.

Methods (by class)

• default: Plot an image of gridded data

• spatialGridRmult: method for spatialGridRmult objects

• spatialGridAcomp: method for spatialGridAcomp objects

Examples

## Not run: 
getTellus(cleanup=TRUE, TI=TRUE)
load("Tellus_TI.RData")
head(Tellus_TI)
coords = as.matrix(Tellus_TI[,1:2])
compo = compositions::acomp(Tellus_TI[,3:7])
dt = spatialGridAcomp(coords=coords, compo=compo)
image_cokriged(dt, ivar="MgO") # equi-spaced
image_cokriged(dt, ivar="MgO", breaks = NULL) # equi-probable

## End(Not run)

Check for any anisotropy class

Description

Check that an object contains a valid specification of anisotropy, in any form

Usage

is.anisotropySpecification(x)

Arguments

Value

a logical, TRUE if the object is an anisotropy specification; FALSE otherwise

See Also

Other anisotropy: AnisotropyRangeMatrix(), AnisotropyScaling(), anis_GSLIBpar2A(), as.AnisotropyRangeMatrix(), as.AnisotropyScaling()

Examples

a =  anis2D_par2A(0.5, 30)
a
is.anisotropySpecification(a)

Check for anisotropy of a theoretical variogram

Description

Checks for anisotropy of a theoretical variogram or covariance function model

Usage

is.isotropic(x, tol = 1e-10, ...)

Arguments

Value

Generic function. Returns of boolean answering the question of the name, or NA if object x does not contain a known theoretical variogram

Length, and number of columns or rows

Description

Provide number of structures, and nr of variables of an LMC of class gmCgram

Usage

## S3 method for class 'gmCgram'
length(x)

Arguments

Value

length returns the number of structures (nugget not counted), while ncol and nrow return these values for the nugget (assuming that they will be also valid for the sill).

See Also

Other gmCgram functions: [.gmCgram(), [[.gmCgram(), as.function.gmCgram(), as.gmCgram.variogramModelList(), ndirections(), plot.gmCgram(), variogramModelPlot()

Examples

utils::data("variogramModels")
v1 = setCgram(type=vg.Gau, sill=diag(3)+0.5, anisRanges = 2*diag(c(3,0.5)))
v2 = setCgram(type=vg.Exp, sill=0.3*diag(3), anisRanges = 0.5*diag(2))
vm = v1+v2
length(vm)
ncol(vm)
nrow(vm)

Empirical logratio variogram calculation

Description

Calculation of an empirical logratio variogram (aka variation-variogram)

Usage

logratioVariogram(data, ...)

Arguments

Value

The output of this function depends on the number of azimuth values provided: if it is one single value (or if you explicitly call logratioVariogram.default) the result is a list of class "logratioVariogram" with the following elements

• vg: A nbins x D x D array containing the logratio variograms

• h: A nbins x D x D array containing the mean distance the value is computed on.

• n: A nbins x D x D array containing the number of nonmissing pairs used for the corresponding value. If azimuth is a vector of directions, then the result is a matrix-like list of "logratioVariogram" objects. Each element of the mother list (or column of the matrix) is the variogram of one direction. The output has a compound class c("logratioVariogramAnisotropy", "logratioVariogram").

Examples

data("jura", package="gstat")
X = jura.pred[,1:2]
Zc = compositions::acomp(jura.pred[,7:13])
vg = logratioVariogram(data=Zc, loc=X)
class(vg)
summary(vg)
vg = logratioVariogram(data=Zc, loc=X, azimuth=c(0,90))
class(vg)
summary(vg)

Logratio variogram of a compositional data

Description

gmGeostats reimplementation of the compositions::logratioVariogram function

Usage

## S4 method for signature 'acomp'
logratioVariogram(
  data = comp,
  loc,
  ...,
  azimuth = 0,
  azimuth.tol = 180/length(azimuth),
  comp = NULL
)

Arguments

Value

a "logratioVariogram" object, or a "logratioVariogramAnisotropy" object if you provide more than one azimuth. See logratioVariogram() for details and examples.

Variogram method for gmSpatialModel objects

Description

Compute the empirical variogram of the conditioning data contained in a gmSpatialModel object

Usage

logratioVariogram_gmSpatialModel(
  data,
  ...,
  azimuth = 0,
  azimuth.tol = 180/length(azimuth)
)

variogram_gmSpatialModel(object, methodPars = NULL, ...)

Arguments

Value

Currently the function is just a convenience wrapper on the variogram calculation functionalities of package "gstat", and returns objects of class "gstatVariogram". Check the help of gstat::variogram for further information. In the near future, methods will be created, which will depend on the properties of the two arguments provided, object and methodPars.

Functions

• logratioVariogram_gmSpatialModel: logratio variogram method (see logratioVariogram() for details)

Construct a Gaussian gmSpatialModel for regionalized compositions

Description

Construct a regionalized compositional data container to be used for Gaussian-based geostatistics: variogram modelling, cokriging and simulation.

Usage

make.gmCompositionalGaussianSpatialModel(
  data,
  coords = attr(data, "coords"),
  V = "ilr",
  prefix = NULL,
  model = NULL,
  beta = model$beta,
  formula = model$formula,
  ng = NULL,
  nmax = ng$nmax,
  nmin = ng$nmin,
  omax = ng$omax,
  maxdist = ng$maxdist,
  force = ng$force
)

Arguments

Value

A "gmSpatialModel" object with all information provided appropriately structured. See gmSpatialModel.

See Also

SequentialSimulation(), TurningBands() or CholeskyDecomposition() for specifying the exact simulation method and its parameters, predict_gmSpatialModel for running predictions or simulations

Other gmSpatialModel: Predict(), as.gmSpatialModel(), gmSpatialModel-class, make.gmCompositionalMPSSpatialModel(), make.gmMultivariateGaussianSpatialModel()

Examples

data("jura", package="gstat")
X = jura.pred[1:20,1:2]
Zc = compositions::acomp(jura.pred[1:20,7:13])
make.gmCompositionalGaussianSpatialModel(data=Zc, coords=X, V="alr")

Construct a Multi-Point gmSpatialModel for regionalized compositions

Description

Construct a regionalized compositional data container to be used for multipoint geostatistics.

Usage

make.gmCompositionalMPSSpatialModel(
  data,
  coords = attr(data, "coords"),
  V = "ilr",
  prefix = NULL,
  model = NULL
)

Arguments

Value

A "gmSpatialModel" object with all information provided appropriately structured. See gmSpatialModel.

See Also

DirectSamplingParameters() for specifying a direct simulation method parameters, predict_gmSpatialModel for running the simulation

Other gmSpatialModel: Predict(), as.gmSpatialModel(), gmSpatialModel-class, make.gmCompositionalGaussianSpatialModel(), make.gmMultivariateGaussianSpatialModel()

Construct a Gaussian gmSpatialModel for regionalized multivariate data

Description

Construct a regionalized multivariate data container to be used for Gaussian-based geostatistics: variogram modelling, cokriging and simulation.

Usage

make.gmMultivariateGaussianSpatialModel(
  data,
  coords = attr(data, "coords"),
  model = NULL,
  beta = model$beta,
  formula = model$formula,
  ng = NULL,
  nmax = ng$nmax,
  nmin = ng$nmin,
  omax = ng$omax,
  maxdist = ng$maxdist,
  force = ng$force
)

Arguments

Value

A "gmSpatialModel" object with all information provided appropriately structured. See gmSpatialModel.

See Also

SequentialSimulation(), TurningBands() or CholeskyDecomposition() for specifying the exact simulation method and its parameters, predict_gmSpatialModel for running predictions or simulations

Other gmSpatialModel: Predict(), as.gmSpatialModel(), gmSpatialModel-class, make.gmCompositionalGaussianSpatialModel(), make.gmCompositionalMPSSpatialModel()

Examples

data("jura", package="gstat")
X = jura.pred[,1:2]
Zc = jura.pred[,7:13]
make.gmMultivariateGaussianSpatialModel(data=Zc, coords=X)

Mean accuracy

Description

Mean method for accuracy curves, after Deutsch (1997)

Usage

## S3 method for class 'accuracy'
mean(x, ...)

Arguments

Value

The value of the mean accuracy after Deutsch (1997); details can be found in precision().

See Also

Other accuracy functions: accuracy(), plot.accuracy(), precision(), validate(), xvErrorMeasures.default(), xvErrorMeasures()

Average measures of spatial decorrelation

Description

Compute average measres of spatial decorrelation out of the result of a call to spatialDecorrelation()

Usage

## S3 method for class 'spatialDecorrelationMeasure'
mean(x, ...)

Arguments

Value

a mean for each measure available on x, i.e. this can be a vector. The output is named.

See Also

spatialDecorrelation() for an example

Number of directions of an empirical variogram

Description

Returns the number of directions at which an empirical variogram was computed

Usage

ndirections(x)

## Default S3 method:
ndirections(x)

## S3 method for class 'azimuth'
ndirections(x)

## S3 method for class 'azimuthInterval'
ndirections(x)

## S3 method for class 'logratioVariogramAnisotropy'
ndirections(x)

## S3 method for class 'logratioVariogram'
ndirections(x)

## S3 method for class 'gmEVario'
ndirections(x)

## S3 method for class 'gstatVariogram'
ndirections(x)

Arguments

Value

Generic function. It provides the number of directions at which an empirical variogram was computed

Methods (by class)

• default: generic method

• azimuth: method for objects of class "azimuth" (vectors of single angles)

• azimuthInterval: method for objects of class "azimuthInterval" (data.frames of intervals for angles)

• logratioVariogramAnisotropy: method for empirical logratio variograms with anisotropy

• logratioVariogram: method for empirical logratio variograms without anisotropy

• gmEVario: method for empirical gmGeostats variograms

• gstatVariogram: method for empirical gstat variograms

See Also

logratioVariogram(), gstat::variogram()

Other gmEVario functions: as.gmEVario.gstatVariogram(), gsi.EVario2D(), gsi.EVario3D(), plot.gmEVario(), variogramModelPlot()

Other gmCgram functions: [.gmCgram(), [[.gmCgram(), as.function.gmCgram(), as.gmCgram.variogramModelList(), length.gmCgram(), plot.gmCgram(), variogramModelPlot()

Test for lack of spatial correlation

Description

Permutation test for checking lack of spatial correlation.

Usage

noSpatCorr.test(Z, ...)

## S3 method for class 'data.frame'
noSpatCorr.test(Z, X, ...)

## Default S3 method:
noSpatCorr.test(Z, ...)

## S3 method for class 'matrix'
noSpatCorr.test(
  Z,
  X,
  R = 299,
  maxlag0 = 0.1 * max(as.matrix(dist(X))),
  minlagInf = 0.25 * max(as.matrix(dist(X))),
  ...
)

Arguments

Value

Produces a test of lack of spatial correlation by means of permutations. The test statistic is based on the smallest eigenvalue of the generalised eigenvalues of the matrices of covariance for short range and for long range.

Methods (by class)

• data.frame: Test for lack of spatial correlation

• default: Test for lack of spatial correlation, works only for Spatial objects with a "data" slot

• matrix: Test for lack of spatial correlation

Examples

data("jura", package="gstat")
X = jura.pred[, 1:2]
Z = data.frame(compositions::ilr(jura.pred[,-(1:6)]))
## Not run:  
noSpatCorr.test(Z=Z, X=X)
# now destroy the spatial structure reshuffling the coordinates:
ip = sample(nrow(X))
noSpatCorr.test(Z=Z, X=X[ip,]) 

## End(Not run)

Multiple maps Matrix of maps showing different combinations of components of a composition, user defined

Description

Multiple maps Matrix of maps showing different combinations of components of a composition, user defined

Usage

pairsmap(data, ...)

## S3 method for class 'SpatialPointsDataFrame'
pairsmap(data, ...)

## Default S3 method:
pairsmap(
  data,
  loc,
  colscale = rev(rainbow(10, start = 0, end = 4/6)),
  cexrange = c(0.1, 2),
  scale = rank,
  commonscale = FALSE,
  mfrow = c(floor(sqrt(ncol(data))), ceiling(ncol(data)/floor(sqrt(ncol(data))))),
  foregroundcolor = "black",
  closeplot = TRUE,
  ...
)

Arguments

Value

The function is primarily called for producing a matrix of plots of each component of a multivariate data set, such as a compositional data set. Each plot represents a map whose symbols are colored and sized according to a z-scale controlled by one of the variables of the data set. It can be used virtually with any geometry, any kind of data (compositional, positive, raw, etc). This function returns invisibly the graphical parameters that were active prior to calling this function. This allows the user to add further stuff to the plots (mostly, using par(mfg=c(i,j)) to plot on the diagram (i,j)), or else freeze the plot (by wrapping the call to pwlrmap on a call to par).

Methods (by class)

• SpatialPointsDataFrame: Multiple maps

• default: Multiple maps

Examples

data("Windarling")
library("compositions")
coords = as.matrix(Windarling[,c("Easting","Northing")])
compo = Windarling[,c("Fe","Al2O3","SiO2", "Mn", "P")]
compo$Rest = 100-rowSums(compo)
compo = acomp(compo)
pairsmap(data=clr(compo), loc=coords) # clr
pairsmap(data=compo, loc=coords) # closed data

Plot method for accuracy curves

Description

Plot an accuracy curve out of the outcome of accuracy().

Usage

## S3 method for class 'accuracy'
plot(
  x,
  xlim = c(0, 1),
  ylim = c(0, 1),
  xaxs = "i",
  yaxs = "i",
  type = "o",
  col = "red",
  asp = 1,
  xlab = "confidence",
  ylab = "coverage",
  pty = "s",
  main = "accuracy plot",
  colref = col[1],
  ...
)

Arguments

Value

Nothing, called to plot the accuracy curve \{(p_i, \pi_i), i=1,2, \ldots, I\}, where \{p_i\} are a sequence of nominal confidence of prediction intervals and each \pi_i is the actual coverage of an interval with nominal confidence p_i.

See Also

Other accuracy functions: accuracy(), mean.accuracy(), precision(), validate(), xvErrorMeasures.default(), xvErrorMeasures()

Draw cuves for covariance/variogram models

Description

Represent a gmCgram object as a matrix of lines in several plots

Usage

## S3 method for class 'gmCgram'
plot(
  x,
  xlim.up = NULL,
  xlim.lo = NULL,
  vdir.up = NULL,
  vdir.lo = NULL,
  xlength = 200,
  varnames = colnames(x$nugget),
  add = FALSE,
  commonAxis = TRUE,
  cov = TRUE,
  closeplot = TRUE,
  ...
)

Arguments

Value

This function is called for its side effect of producing a plot: the plot will be open to further changes if you provide closeplot=FALSE. Additionally, the function invisibly returns the graphical parameters that were active before starting the plot. Hence, if you want to freeze a plot and not add anymore to it, you can do par(plot(x, closeplot=FALSE, ...)), or plot(x, closeplot=TRUE, ...). If you want to further add stuff to it, better just call plot(x, closeplot=FALSE,...). The difference is only relevant when working with the screen graphical device.

See Also

Other gmCgram functions: [.gmCgram(), [[.gmCgram(), as.function.gmCgram(), as.gmCgram.variogramModelList(), length.gmCgram(), ndirections(), variogramModelPlot()

Examples

utils::data("variogramModels")
v1 = setCgram(type=vg.Gau, sill=diag(3)-0.5, anisRanges = 2*diag(c(3,0.5)))
v2 = setCgram(type=vg.Exp, sill=0.3*diag(3), anisRanges = 0.5*diag(2))
vm = v1+v2
plot(vm)
plot(vm, cov=FALSE)

Plot empirical variograms

Description

Flexible plot of an empirical variogram of class gmEVario

Usage

## S3 method for class 'gmEVario'
plot(
  x,
  xlim.up = NULL,
  xlim.lo = NULL,
  vdir.up = NULL,
  vdir.lo = NULL,
  varnames = dimnames(x$gamma)[[2]],
  type = "o",
  add = FALSE,
  commonAxis = TRUE,
  cov = attr(x, "type") == "covariance",
  closeplot = TRUE,
  ...
)

Arguments

Value

invisibly, the graphical parameters active before calling the function. This is useful for freezing the plot if you provided closeplot=FALSE.

How to use arguments vdir.lo and vdir.up? Each empirical variogram x has been computed along certain distances, recorded in its attributes and retrievable with command ndirections.

See Also

Other gmEVario functions: as.gmEVario.gstatVariogram(), gsi.EVario2D(), gsi.EVario3D(), ndirections(), variogramModelPlot()

Examples

library(gstat)
data("jura", package = "gstat")
X = as.matrix(jura.pred[,1:2])
Z = as.matrix(jura.pred[,c("Zn","Cd","Pb")])
vge = gsi.EVario2D(X,Z)
plot(vge)
plot(vge, pch=22, lty=1, bg="grey")

Plot variogram lines of empirical directional logratio variograms

Description

Plots an "logratioVariogramAnisotropy" object in a series of panels, with each direction represented as a broken line.

Usage

## S3 method for class 'logratioVariogramAnisotropy'
plot(
  x,
  azimuths = colnames(x),
  col = rev(rainbow(length(azimuths))),
  type = "o",
  V = NULL,
  lty = 1,
  pch = 1:length(azimuths),
  model = NULL,
  figsp = 0,
  closeplot = TRUE,
  ...
)

Arguments

Value

Nothing. The function is called to create a plot.

Examples

data("jura", package="gstat")
X = jura.pred[,1:2]
Zc = compositions::acomp(jura.pred[,7:9])
vg = logratioVariogram(data=Zc, loc=X, azimuth=c(0:3)*45, 
                         azimuth.tol=22.5)
plot(vg) 

Plotting method for swarmPlot objects

Description

Produces a plot for objects obtained by means of swarmPlot(). These are lists of two-element list, respectively giving the (x,y) coordinates of a series of points defining a curve. They are typically obtained out of a DataFrameStack(), and each curve represents a a certain relevant summary of one of the elements of the stack. All parameters of graphics::plot.default() are available, plus the following

Usage

## S3 method for class 'swarmPlot'
plot(
  x,
  type = "l",
  col = "grey",
  xlim = NULL,
  ylim = NULL,
  xlab = "x",
  ylab = "y",
  main = NULL,
  asp = NA,
  sub = NULL,
  log = "",
  ann = par("ann"),
  axes = TRUE,
  frame.plot = axes,
  bg = NA,
  pch = 1,
  cex = 1,
  lty = 1,
  lwd = 1,
  add = FALSE,
  ...
)

Arguments

Value

Nothing. The function is called to make a plot.

Examples

dm = list(point=1:100, var=LETTERS[1:2], rep=paste("r",1:5, sep=""))
ar = array(rnorm(1000), dim=c(100,2,5), dimnames = dm)
dfs = DataFrameStack(ar, stackDim="rep")
rs = swarmPlot(dfs, PLOTFUN=function(x) density(x[,1]), .plotargs=FALSE)
plot(rs, col="yellow", lwd=3)

Precision calculations

Description

Precision calculations

Usage

precision(x, ...)

## S3 method for class 'accuracy'
precision(x, ...)

Arguments

Value

output depends on input and meaning of the function (the term precision is highly polysemic)

Methods (by class)

• accuracy: Compute precision and goodness for accuracy curves, after Deutsch (1997), using the accuracy curve obtained with accuracy(). This returns a named vector with two values, one for precision and one for goodness.

  Mean accuracy, precision and goodness were defined by Deutsch (1997) for an accuracy curve \{(p_i, \pi_i), i=1,2, \ldots, I\}, where \{p_i\} are a sequence of nominal confidence of prediction intervals and each \pi_i is the actual coverage of an interval with nominal confidence p_i. Out of these values, the mean accuracy (see mean.accuracy()) is computed as

  A = \int_{0}^{1} I\{(\pi_i-p_i)>0\} dp,

  where the indicator I\{(\pi_i-p_i)>0\} is 1 if the condition is satisfied and 0 otherwise. Out of it, the area above the 1:1 bisector and under the accuracy curve is the precision P = 1-2\int_{0}^{1} (\pi_i-p_i)\cdot I\{(\pi_i-p_i)>0\} dp, which only takes into account those points of the accuracy curve where \pi_i>p_i. To consider the whole curve, goodness can be used

  G = 1-\int_{0}^{1} (\pi_i-p_i)\cdot (3\cdot I\{(\pi_i-p_i)>0\}-2) dp.

See Also

Other accuracy functions: accuracy(), mean.accuracy(), plot.accuracy(), validate(), xvErrorMeasures.default(), xvErrorMeasures()

Compute model variogram values Evaluate the variogram model provided at some lag vectors

Description

Compute model variogram values

Evaluate the variogram model provided at some lag vectors

Usage

## S3 method for class 'LMCAnisCompo'
predict(object, newdata, ...)

Arguments

Value

an array of dimension (nr of lags, D, D) with D the number of variables in the model.

Examples

data("jura", package="gstat")
Zc = compositions::acomp(jura.pred[,7:9])
lrmd = LMCAnisCompo(Zc, models=c("nugget", "sph"), azimuths=c(0,45))
predict(lrmd, outer(0:5, c(0,1)))

Predict method for generalised diagonalisation objects

Description

Predict method for generalised diagonalisation objects

Usage

## S3 method for class 'genDiag'
predict(object, newdata = NULL, ...)

Arguments

Value

A data set or compositional object of the nature of the original data used for creating the genDiag object.

See Also

Other generalised Diagonalisations: Maf(), coloredBiplot.genDiag()

Examples

data("jura", package="gstat")
juracomp = compositions::acomp(jura.pred[, -(1:6)]) 
lrvg = logratioVariogram(data=juracomp, loc=jura.pred[,1:2])
mf = Maf(juracomp, vg=lrvg)
mf
biplot(mf)
predict(mf) 
unclass(predict(mf)) - unclass(juracomp) # predict recovers the original composition

Print method for mask objects

Description

Print method for mask objects. See constructMask() for examples. If you want to see the whole content of the mask, then use unclass(...)

Usage

## S3 method for class 'mask'
print(x, ...)

Arguments

Value

the summary of number of nodes inside/outside the mask

See Also

Other masking functions: constructMask(), getMask(), setMask(), unmask()

Compositional maps, pairwise logratios Matrix of maps showing different combinations of components of a composition, in pairwise logratios

Description

Compositional maps, pairwise logratios Matrix of maps showing different combinations of components of a composition, in pairwise logratios

Usage

pwlrmap(
  loc,
  comp,
  colscale = rev(rainbow(10, start = 0, end = 4/6)),
  cexrange = c(0.1, 2),
  scale = rank,
  commonscale = FALSE,
  foregroundcolor = "black",
  closeplot = TRUE
)

Arguments

Value

The function is primarily called for producing a matrix of (D,D) plots of the D-part compositional samples, where at each plot we represent a map whose symbols are colored and sized according to a z-scale controlled by a different logratio. For each plot, this is the logratio of the row variable by the column variable. However, in case that closeplot=FALSE, this function returns invisibly the graphical parameters that were active prior to calling this function. This allows the user to add further stuff to the plots (mostly, using par(mfg=c(i,j)) to plot on the diagram (i,j)), or manually freeze the plot (by wrapping the call to pwlrmap on a call to par).

Examples

data("Windarling")
coords = as.matrix(Windarling[,c("Easting","Northing")])
compo = Windarling[,c("Fe","Al2O3","SiO2", "Mn", "P")]
compo$Rest = 100-rowSums(compo)
compo = compositions::acomp(compo) 
# in quantiles (default, ranking controls color and size)
pwlrmap(coords, compo) 

# in logratios (I=identity)
pwlrmap(coords, compo, scale=I)
# in ratios (i.e., apply exp)
pwlrmap(coords, compo, scale=exp)  
# use only color, no change in symbol size
pwlrmap(coords, compo, cexrange=c(1,1)) 
# change all
pwlrmap(coords, compo, commonscale=TRUE, cexrange=c(1.2,1.2), 
                    colscale=rev(rainbow(40, start=0, end=4/6))) 

Generate D-variate variogram models

Description

Function to set up D-variate variogram models based on model type, the variogram parameters sill and nugget and a matrix describing the anisotropy of the range.

Usage

setCgram(type, nugget = sill * 0, sill, anisRanges, extraPar = 0)

Arguments

Details

The argument type must be an integer indicating the model to be used. Some constants are available to make reading code more understandable. That is, you can either write 1, vg.sph, vg.Sph or vg.Spherical, they will all work and produce a spherical model. The same applies for the following models: vg.Gauss = vg.Gau = vg.gau = 0; ⁠vg.Exponential = vg.Exp = vg. exp = 2⁠. These constants are available after calling data("variogramModels"). No other model is currently available, but this data object will be regularly updated. The constant vector gsi.validModels contains all currently valid models.

Argument anisRange expects a matrix $M$ such that

h^2 = (\mathbf{x}_i-\mathbf{x}_j)\cdot M^{-1}\cdot (\mathbf{x}_i-\mathbf{x}_j)^t

is the (square of) the lag distance to be fed into the correlation function.

Value

an object of class "gmCgram" containing the linear model of corregionalization of the nugget and the structure given.

Examples

utils::data("variogramModels") # shortcut for all model constants
v1 = setCgram(type=vg.Gau, sill=diag(2), anisRanges = 3*diag(c(3,1)))
v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
vm = v1+v2
plot(vm)

Set or get the ordering of a grid

Description

Specify or retrieve the ordering in which a grid is stored in a vector (or matrix).

Usage

setGridOrder(x, refpoint, cycle)

setGridOrder_sp(x, G = 2)

setGridOrder_array(x, G = 2)

gridOrder_sp(G = 2)

gridOrder_gstat(G = 2)

gridOrder_array(G = 2)

gridOrder_GSLib(G = 2)

Arguments

Details

A "gridOrder" attribute is a list consisting of two named elements:

refpoint

one of "topleft", "bottomleft", "topright" or "bottomright" in 2D, or also of "topleftsurf", "bottomleftsurf", "toprightsurf", "bottomrightsurf", "topleftdeep","bottomleftdeep", "toprightdeep" or "bottomrightdeep" in 3D ("deep" is accessory, i.e. "topleft"== "topleftdeep"), indicating the location on the grid of the first point of the object x

cycle

a permutation of 1:G indicating in which order run the dymensions, from faster to slower

Thus, a conventional ordering of a (nX*nY)-element vector into a matrix to plot with graphics::image() corresponds to an refpoint="bottomleft" and cycle=1:2, i.e. start with the lower left corner and run first by rows (eastwards), then by columns (northwards). This is constructed by gridOrder_array(G), and can be directly set to an object x by setGridOrder_array(x,G); gridOrder_GSLib is an alias for gridOrder_array.

The grids from package "sp" (and many other in R), on the contrary follow the convention refpoint="topleft" and cycle=1:2, i.e. start with the upper left corner and run first by rows (eastwards), then by columns (southwards). This is constructed by gridOrder_sp(G), and can be directly set to an object x by setGridOrder_sp(x,G); gridOrder_gstat is an alias for gridOrder_sp.

Value

setGridOrder(x,...) returns the object x with the grid order description attached as an attribute "gridOrder"; getGridOrder(x) retrieves this attribute and returns it.

Functions

• setGridOrder_sp: Set or get the ordering of a grid

• setGridOrder_array: Set or get the ordering of a grid

• gridOrder_sp: Set or get the ordering of a grid

• gridOrder_gstat: Set or get the ordering of a grid

• gridOrder_array: Set or get the ordering of a grid

• gridOrder_GSLib: Set or get the ordering of a grid

See Also

sortDataInGrid() for ways of reordering a grid

Examples

gt = sp::GridTopology(cellcentre.offset=c(1,11), cellsize=c(1,1), cells.dim=c(5,3))
sp::coordinates(sp::SpatialGrid(grid=gt))
gridOrder_sp(2)

Set a mask on an object

Description

Set a mask on an object See constructMask() for examples on how to construct masks.

Usage

setMask(x, ...)

## Default S3 method:
setMask(x, mask, coordinates = 1:2, ...)

## S3 method for class 'data.frame'
setMask(x, mask, coordinates = 1:2, ...)

## S3 method for class 'DataFrameStack'
setMask(x, mask, coordinates = attr(x, "coordinates"), ...)

## S3 method for class 'SpatialGrid'
setMask(x, mask, ...)

## S3 method for class 'GridTopology'
setMask(x, mask, ...)

## S3 method for class 'SpatialPoints'
setMask(x, mask, ...)

Arguments

Value

The object x appropriately masked (for the setter methods).

Methods (by class)

• default: Set a mask on an object

• data.frame: Set a mask on a data.frame object

• DataFrameStack: Set a mask on a DataFrameStack object

• SpatialGrid: Set a mask on a SpatialGrid object

• GridTopology: Set a mask on a GridTopology object

• SpatialPoints: Set a mask on a SpatialPoints object

See Also

Other masking functions: constructMask(), getMask(), print.mask(), unmask()

Reorder data in a grid

Description

Reorder the data in a compact grid, changing between ordering specifications

Usage

sortDataInGrid(
  x,
  grid = attr(x, "grid"),
  orderIn = attr(x, "gridOrder"),
  orderOut = list(refpoint = "bottomleft", cycle = 1:2)
)

Arguments

Value

the data from x (typically a matrix), but reordered as orderOut

See Also

setGridOrder() for ways of specifying the grid ordering

Examples

## Not run: 
getTellus(cleanup=TRUE, TI=TRUE)
load("Tellus_TI.RData")
coords = as.matrix(Tellus_TI[,1:2])
compo = compositions::acomp(Tellus_TI[,3:7])
dt = spatialGridAcomp(coords=coords, compo=compo)
image_cokriged(dt, ivar="MgO", breaks = NULL) 
x = sort(unique(coords[,1]))
y = sort(unique(coords[,2]))
x0 = c(min(x), min(y))
Ax = c(mean(diff(x)), mean(diff(y)))
n = c(length(x), length(y))
gr = sp::GridTopology(cellcentre.offset=x0, cellsize=Ax, cells.dim=n)
dt0 = sortDataInGrid(Tellus_TI, grid=gr, orderIn=gridOrder_array(2), 
                    orderOut=list(refpoint="bottomright", cycle=2:1))
coords = as.matrix(dt0[,1:2])
compo = compositions::acomp(dt0[,3:7])
spatialGridAcomp(coords=coords, compo=compo)

## End(Not run)

Compute diagonalisation measures

Description

Compute one or more diagonalisation measures out of an empirical multivariate variogram.

Usage

spatialDecorrelation(vgemp, ...)

## S3 method for class 'gstatVariogram'
spatialDecorrelation(
  vgemp,
  vgemp0 = NULL,
  method = "add",
  quadratic = method[1] != "rdd",
  ...
)

## S3 method for class 'logratioVariogram'
spatialDecorrelation(vgemp, vgemp0 = NULL, method = "add", ...)

## S3 method for class 'gmEVario'
spatialDecorrelation(vgemp, vgemp0 = NULL, method = "add", ...)

Arguments

Details

The three measures provided are

absolute deviation from diagonality ("add")

defined as the sum of all off-diagonal elements of the variogram, possibly squared ($p=2$ if quadratic=TRUE the default; otherwise $p=1$)

\zeta(h)=\sum_{k=1}^n\sum_{j\neq k}^n \gamma_{k,j}^p(h)

relative deviation from diagonality ("rdd")

comparing the absolute sum of off-diagonal elements with the sum of the diagonal elements of the variogram, each possibly squared ($p=2$ if quadratic=TRUE; otherwise $p=1$ the default)

\tau(h)=\frac{\sum_{k=1}^n\sum_{j \neq k}^n |\gamma_{k,j}(h)|^p}{\sum_{k=1}^n|\gamma_{k,k}(h)|^p}

spatial diagonalisation efficiency ("sde")

is the only one requiring vgemp0, because it compares an initial state with a diagonalised state of the variogram system

\kappa(h)=1- \frac{\sum_{k=1}^n\sum_{j \neq k}^n |\gamma_{k,j}(h)|^p}{\sum_{k=1}^n\sum_{j \neq k}^n |\gamma_{(0)k,j}(h)|^p }

The value of $p$ is controlled by the first value of method. That is, the results with method=c("rdd", "add") are not the same as those obtained with method=c("add", "rdd"), as in the first case $p=1$ and in the second case $p=2$.

Value

an object of a similar nature to vgemp, but where the desired quantities are reported for each lag. This can then be plotted or averages be computed.

Methods (by class)

• gstatVariogram: Compute diagonalisation measures

• logratioVariogram: Compute diagonalisation measures

• gmEVario: Compute diagonalisation measures

Examples

data("jura", package="gstat")
X = jura.pred[, 1:2]
Z = jura.pred[,-(1:6)]
gm1 = make.gmCompositionalGaussianSpatialModel(data=Z, coords=X, V="alr")
vg1 = variogram(as.gstat(gm1)) 
(r1 = spatialDecorrelation(vg1, method=c("add", "rdd")))
plot(r1)
mean(r1)
require("compositions")
pc = princomp(acomp(Z))
v = pc$loadings
colnames(v)=paste("pc", 1:ncol(v), sep="")
gm2 = make.gmCompositionalGaussianSpatialModel(data=Z, coords=X, V=v, prefix="pc")
vg2 = variogram(as.gstat(gm2)) 
(r2 = spatialDecorrelation(vg2, method=c("add", "rdd")))
plot(r2)
mean(r2)
(r21 = spatialDecorrelation(vg2, vg1, method="sde") )
plot(r21)
mean(r21)

Construct a regionalized composition / reorder compositional simulations

Description

Connect some coordinates to a composition (of hard data, of predictions or of simulations); currently, the coordinates are stored in an attribute and the dataset is given a complex S3 class. This functionality will change in the future, to make use of package "sp" classes.

Usage

spatialGridAcomp(coords, compo, dimcomp = 2, dimsim = NA)

Arguments

Value

A (potentially transposed/aperm-ed) matrix of class c("spatialGridAcomp","acomp") with the coordinates in an extra attribute "coords".

See Also

image_cokriged.spatialGridAcomp() for an example; gsi.gstatCokriging2compo() to restructure the output from gstat::predict.gstat() confortably

Construct a regionalized multivariate data

Description

Connect some coordinates to a multivariate data set (of hard data, of predictions or of simulations); currently, the coordinates are stored in an attribute and the dataset is given a complex S3 class. This functionality will change in the future, to make use of package "sp" classes.

Usage

spatialGridRmult(coords, data, dimcomp = 2, dimsim = NA)

Arguments

Value

A (potentially transposed/aperm-ed) matrix of class c("spatialGridAcomp","acomp") with the coordinates in an extra attribute "coords".

See Also

image_cokriged.spatialGridRmult() for an example; gsi.gstatCokriging2rmult() to restructure the output from gstat::predict.gstat() confortably

Spectral colors palette based on the RColorBrewer::brewer.pal(11,"Spectral")

Description

Spectral colors palette based on the RColorBrewer::brewer.pal(11,"Spectral")

Usage

spectralcolors(n)

Arguments

Value

a palette, i.e. a list of colors, from dark blue to dark red over pale yellow.

Examples

(cls=spectralcolors(20))

Spherifying transform Compute a transformation that spherifies a certain data set

Description

Spherifying transform Compute a transformation that spherifies a certain data set

Usage

sphTrans(Y, ...)

## Default S3 method:
sphTrans(Y, weights = NULL, p = 1:ncol(Y), ...)

Arguments

Value

a function with arguments (x, inv=FALSE), where x will be the data to apply the transformation to, and inv=FALSE will indicate if the direct or the inverse transformation is desired. This function applied to the same data returns a translated, rotated and scaled, so that the new scores are centered, have variance 1, and no correlation.

Methods (by class)

• default: Spherifying transform

Author(s)

K. Gerald van den Boogaart, Raimon Tolosana-Delgado

See Also

ana, anaBackward, sphTrans

Examples

library(compositions)
data("jura", package="gstat")
Y = acomp(jura.pred[,c(10,12,13)])
oldpar = par(mfrow = c(1,1))
plot(Y)
sph = sphTrans(Y)
class(sph)
z = sph(Y)
plot(z)
par(oldpar)
cor(cbind(z, ilr(Y)))
colMeans(cbind(z, ilr(Y)))

Get/set name/index of (non)stacking dimensions

Description

Return (or set) the name or index of either the stacking dimension, or else of the non-stacking dimension (typically, the dimension runing through the variables)

Usage

stackDim(x, ...)

## S3 method for class 'DataFrameStack'
stackDim(x, ...)

noStackDim(x, ...)

## Default S3 method:
noStackDim(x, ...)

stackDim(x) <- value

## Default S3 replacement method:
stackDim(x) <- value

Arguments

Value

the index or the name of the asked dimension.

Functions

• stackDim.DataFrameStack: Get/set name/index of (non)stacking dimensions

• noStackDim: Get/set name/index of (non)stacking dimensions

• noStackDim.default: Get/set name/index of (non)stacking dimensions

• stackDim<-: Get/set name/index of (non)stacking dimensions

• stackDim<-.default: Get/set name/index of (non)stacking dimensions

Examples

ar = array(1:30, dim = c(5,2,3), dimnames=list(obs=1:5, vars=c("A","B"), rep=1:3))
dfs = DataFrameStack(ar, stackDim="rep")
dfs
stackDim(dfs)
noStackDim(dfs)
getStackElement(dfs, 1)
stackDim(dfs) <- "vars"
getStackElement(dfs, 1)

Get name/index of the stacking dimension of a Spatial object

Description

Get name/index of the stacking dimension of a Spatial object which data slot is of class DataFrameStack()

Usage

## S4 method for signature 'Spatial'
stackDim(x)

Arguments

Value

see stackDim() for details

See Also

stackDim()

Plot a swarm of calculated output through a DataFrameStack

Description

Take a DataFrameStack() and apply a certain plotting function to each elements of the stack. The result (typically a curve per each stack element), may then be plotted all together

Usage

swarmPlot(
  X,
  MARGIN = stackDim(X),
  PLOTFUN,
  ...,
  .plotargs = list(type = "l"),
  .parallelBackend = NA
)

Arguments

Value

Invisibly, this function returns a list of the evaluation of PLOTFUN on each element of the stack. If .plotargs other than FALSE, then the function calls plot.swarmPlot() to produce a plot.

Examples

dm = list(point=1:100, var=LETTERS[1:2], rep=paste("r",1:5, sep=""))
ar = array(rnorm(1000), dim=c(100,2,5), dimnames = dm)
dfs = DataFrameStack(ar, stackDim="rep")
swarmPlot(dfs, PLOTFUN=function(x) density(x[,1]))

Swath plots

Description

Plots of data vs. one spatial coordinate, with local average spline curve

Usage

swath(data, ...)

## Default S3 method:
swath(
  data,
  loc,
  pch = ".",
  withLoess = TRUE,
  commonScale = TRUE,
  xlab = deparse(substitute(loc)),
  ...,
  mfrow
)

## S3 method for class 'acomp'
swath(
  data,
  loc,
  pch = ".",
  withLoess = TRUE,
  commonScale = NA,
  xlab = deparse(substitute(loc)),
  ...
)

## S3 method for class 'ccomp'
swath(
  data,
  loc,
  pch = ".",
  withLoess = TRUE,
  commonScale = NA,
  xlab = deparse(substitute(loc)),
  ...
)

## S3 method for class 'rcomp'
swath(
  data,
  loc,
  pch = ".",
  withLoess = TRUE,
  commonScale = NA,
  xlab = deparse(substitute(loc)),
  ...
)

Arguments

Value

Nothing, as the function is primarily called to produce a plot

Methods (by class)

• default: swath plot default method

• acomp: Swath plots for acomp objects

• ccomp: Swath plots for ccomp objects

• rcomp: Swath plots for rcomp objects

Examples

data("Windarling")
library("sp")
compo = compositions::acomp(Windarling[,c("Fe","Al2O3","SiO2", "Mn", "P")])
Northing = Windarling$Northing
swath(compo, Northing)
wind.spdf = sp::SpatialPointsDataFrame(sp::SpatialPoints(Windarling[,6:7]), 
     data=compo)
swath(wind.spdf, loc=Northing)

Unmask a masked object

Description

Unmask a masked object, i.e. recover the original grid and extend potential data containers associated to it with NAs. See examples in constructMask()

Usage

unmask(x, ...)

## S3 method for class 'data.frame'
unmask(
  x,
  mask = attr(x, "mask"),
  fullgrid = attr(mask, "fullgrid"),
  forceCheck = is(fullgrid, "GridTopology"),
  ...
)

## S3 method for class 'DataFrameStack'
unmask(
  x,
  mask = attr(x, "mask"),
  fullgrid = attr(mask, "fullgrid"),
  forceCheck = is(fullgrid, "GridTopology"),
  ...
)

## S3 method for class 'SpatialPixels'
unmask(
  x,
  mask = NULL,
  fullgrid = attr(mask, "fullgrid"),
  forceCheck = FALSE,
  ...
)

## S3 method for class 'SpatialPoints'
unmask(
  x,
  mask = attr(x@data, "mask"),
  fullgrid = attr(mask, "fullgrid"),
  forceCheck = FALSE,
  ...
)

Arguments

Value

The original grid data and extend potential data containers associated to it with NAs. See examples in constructMask(). The nature of the output depends on the nature of x: a "data.frame" produced a "data.frame"; a "unmask.DataFrameStack" produces a "unmask.DataFrameStack"; a "SpatialPoints" produces a "SpatialPoints"; and finally a "SpatialPixels" produces either a "SpatialPixels" or a "SpatialGrid" (if it is full). Note that only in the case that is(x,"SpatialPixels")=TRUE is mask required, for the other methods all arguments have reasonable defaults.

Functions

• unmask: Unmask a masked object

• unmask.DataFrameStack: Unmask a masked object

• unmask.SpatialPixels: Unmask a masked object

• unmask.SpatialPoints: Unmask a masked object

See Also

Other masking functions: constructMask(), getMask(), print.mask(), setMask()

Validate a spatial model

Description

Validate a spatial model by predicting some values. Typically this will be a validation set, or else some subset of the conditing data.

Usage

validate(object, strategy, ...)

## S3 method for class 'LeaveOneOut'
validate(object, strategy, ...)

## S3 method for class 'NfoldCrossValidation'
validate(object, strategy, ...)

Arguments

Value

A data frame of predictions (possibly with kriging variances and covariances, or equivalent uncertainty measures) for each element of the validation set

Methods (by class)

• LeaveOneOut: Validate a spatial model

• NfoldCrossValidation: Validate a spatial model

See Also

Other validation functions: LeaveOneOut, NfoldCrossValidation

Other accuracy functions: accuracy(), mean.accuracy(), plot.accuracy(), precision(), xvErrorMeasures.default(), xvErrorMeasures()

Examples

data("Windarling")
X = Windarling[,c("Easting","Northing")]
Z = compositions::acomp(Windarling[,c(9:12,16)])
gm = make.gmCompositionalGaussianSpatialModel(data=Z, coords=X)
vg = variogram(gm)
md = gstat::vgm(range=30, model="Sph", nugget=1, psill=1)
gs = fit_lmc(v=vg, g=gm, model=md) 
## Not run:  v1 = validate(gs, strategy=LeaveOneOut()) # quite slow 
vs2 = NfoldCrossValidation(nfolds=sample(1:10, nrow(X), replace=TRUE))
vs2
## Not run:  v2 = validate(gs, strategy=vs2) # quite slow 

Quick plotting of empirical and theoretical variograms Quick and dirty plotting of empirical variograms/covariances with or without their models

Description

Quick plotting of empirical and theoretical variograms Quick and dirty plotting of empirical variograms/covariances with or without their models

Usage

variogramModelPlot(vg, ...)

## S3 method for class 'gmEVario'
variogramModelPlot(
  vg,
  model = NULL,
  col = rev(rainbow(ndirections(vg))),
  commonAxis = FALSE,
  newfig = TRUE,
  closeplot = TRUE,
  ...
)

Arguments

Value

The function is primarily called for producing a plot. However, it invisibly returns the graphical parameters active before the call occurred. This is useful for constructing complex diagrams, by adding layers of info. If you want to "freeze" your plot, embed your call in another call to par, e.g. par(variogramModelPlot(...)); if you want to leave the plot open for further changes give the extra argument closeplot=FALSE.

Functions

• variogramModelPlot: Quick plotting of empirical and theoretical variograms

See Also

logratioVariogram()

Other variogramModelPlot: variogramModelPlot.gstatVariogram(), variogramModelPlot.logratioVariogram()

Other gmEVario functions: as.gmEVario.gstatVariogram(), gsi.EVario2D(), gsi.EVario3D(), ndirections(), plot.gmEVario()

Other gmCgram functions: [.gmCgram(), [[.gmCgram(), as.function.gmCgram(), as.gmCgram.variogramModelList(), length.gmCgram(), ndirections(), plot.gmCgram()

Examples

utils::data("variogramModels")
v1 = setCgram(type=vg.Gau, sill=diag(3)+0.5, anisRanges = 5e-1*diag(c(3,0.5)))
v2 = setCgram(type=vg.Exp, sill=0.3*diag(3), anisRanges = 5e-2*diag(2))
vm = v1+v2
plot(vm, closeplot=TRUE)
library(gstat)
data("jura", package = "gstat")
X = as.matrix(jura.pred[,1:2])
Z = as.matrix(jura.pred[,c("Zn","Cd","Pb")])
vge = gsi.EVario2D(X,Z)
variogramModelPlot(vge, vm)

Quick plotting of empirical and theoretical variograms Quick and dirty plotting of empirical variograms/covariances with or without their models

Description

Quick plotting of empirical and theoretical variograms Quick and dirty plotting of empirical variograms/covariances with or without their models

Usage

## S3 method for class 'gstatVariogram'
variogramModelPlot(
  vg,
  model = NULL,
  col = rev(rainbow(1 + length(unique(vg$dir.hor)))),
  commonAxis = FALSE,
  newfig = TRUE,
  closeplot = TRUE,
  ...
)

Arguments

Value

The function is primarily called for producing a plot. However, it invisibly returns the graphical parameters active before the call occurred. This is useful for constructing complex diagrams, by giving argument closeplot=FALSE and then adding layers of information. If you want to "freeze" your plot, either give closeplot=TRUE or embed your call in another call to par, e.g. par(variogramModelPlot(...)).

See Also

gstat::plot.gstatVariogram()

Other variogramModelPlot: variogramModelPlot.logratioVariogram(), variogramModelPlot()

Examples

data("jura", package="gstat")
X = jura.pred[,1:2]
Zc = jura.pred[,7:13]
gg = make.gmCompositionalGaussianSpatialModel(Zc, X, V="alr", formula = ~1)
vg = variogram(gg)
md = gstat::vgm(model="Sph", psill=1, nugget=1, range=1.5)
gg = fit_lmc(v=vg, g=gg, model=md)
variogramModelPlot(vg, model=gg)

Quick plotting of empirical and theoretical logratio variograms Quick and dirty plotting of empirical logratio variograms with or without their models

Description

Quick plotting of empirical and theoretical logratio variograms Quick and dirty plotting of empirical logratio variograms with or without their models

Usage

## S3 method for class 'logratioVariogram'
variogramModelPlot(vg, model = NULL, col = rev(rainbow(ndirections(vg))), ...)

Arguments

Value

The function is primarily called for producing a plot. However, it invisibly returns the graphical parameters active before the call occurred. This is useful for constructing complex diagrams, by adding layers of info. If you want to "freeze" your plot, embed your call in another call to par, e.g. par(variogramModelPlot(...)).

See Also

Other variogramModelPlot: variogramModelPlot.gstatVariogram(), variogramModelPlot()

Write a regionalized data set in GSLIB format

Description

Write a regionalized data set in plain text GSLIB format

Usage

write.GSLib(x, file, header = basename(file))

Arguments

Value

The status of closing the file, see close for details, although this is seldom problematic. This function is basically called for its side-effect of writing a data set in the simplified Geo-EAS format that is used in GSLIB.

See Also

http://www.gslib.com/gslib_help/format.html

Examples

data("jura", package="gstat")
## Not run: write.GSLib(jura.pred, file="jurapred.txt")

Cross-validation errror measures

Description

Compute one or more error measures from cross-validation output

Usage

xvErrorMeasures(x, ...)

## S3 method for class 'data.frame'
xvErrorMeasures(
  x,
  observed = x$observed,
  output = "MSDR1",
  univariate = length(dim(observed)) == 0,
  ...
)

## S3 method for class 'DataFrameStack'
xvErrorMeasures(
  x,
  observed,
  output = "ME",
  univariate = length(dim(observed)) == 0,
  ...
)

Arguments

Details

"ME" stands for mean error (average of the differences between true values and predicted values), "MSE" stands for mean square error (average of the square differences between true values and predicted values), and "MSDR" for mean squared deviation ratio (average of the square between true values and predicted values each normalized by its kriging variance). These quantities are classically used in evaluating output results of validation exercises of one single variable. For multivariate cases, "ME" (a vector) and "MSE" (a scalar) work as well, while two different definitions of a multivariate mean squared deviation ratio can be given:

• "MSDR1" is the average Mahalanobis square error (see accuracy() for explanations)

• "MSDR2" is the average univariate "MSDR" over all variables.

Value

If just some of c("ME","MSE","MSDR","MSDR1","MSDR2") are requested, the output is a named vector with the desired quantities. If only "Mahalanobis" is requested, the output is a vector of Mahalanobis square errors. If you mix up things and ask for "Mahalanobis" and some of the quantities mentioned above, the result will be a named list with the requested quantities. (NOTE: some options are not available for x a "DataFrameStack")

Functions

• xvErrorMeasures: Cross-validation errror measures

• xvErrorMeasures.DataFrameStack: Cross-validation errror measures

See Also

Other accuracy functions: accuracy(), mean.accuracy(), plot.accuracy(), precision(), validate(), xvErrorMeasures.default()

Cross-validation errror measures

Description

Compute one or more error measures from cross-validation output

Usage

## Default S3 method:
xvErrorMeasures(x, krigVar, observed, output = "MSDR1", ...)

Arguments

Details

"ME" stands for mean error (average of the differences between true values and predicted values), "MSE" stands for mean square error (average of the square differences between true values and predicted values), and "MSDR" for mean squared deviation ratio (average of the square between true values and predicted values each normalized by its kriging variance). These quantities are classically used in evaluating output results of validation exercises of one single variable. For multivariate cases, see xvErrorMeasures.data.frame().

Value

If just some of c("ME","MSE","MSDR") are requested, the output is a named vector with the desired quantities. If only "Mahalanobis" is requested, the output is a vector of Mahalanobis square errors. If you mix up things and ask for "Mahalanobis" and some of the quantities mentioned above, the result will be a named list with the requested quantities.

See Also

Other accuracy functions: accuracy(), mean.accuracy(), plot.accuracy(), precision(), validate(), xvErrorMeasures()