From c6149a3e7b54ac0b870a21eba2a4f745569a546e Mon Sep 17 00:00:00 2001 From: Solveig Pospiech <s.pospiech@hzdr.de> Date: Sun, 28 Nov 2021 12:11:15 +0100 Subject: [PATCH] setCgram: small improvements in the help-page --- R/variograms.R | 412 +++++++++++++++++++++++++------------------------ 1 file changed, 207 insertions(+), 205 deletions(-) diff --git a/R/variograms.R b/R/variograms.R index 1f68ea5..256bf6f 100644 --- a/R/variograms.R +++ b/R/variograms.R @@ -1,34 +1,36 @@ #### theoretical variogram --------------------- -#' gmGeostats Variogram models -#' set up a D-variate variogram models +#' Generate D-variate variogram models #' -#' @param type constant indicating the model of correlation function -#' @param nugget (DxD)-matrix for the nugget effect +#' @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. +#' +#' @param type model of correlation function. The function expects a constant, e.g. the internal constants 'vg.Gau' for Gaussian model or 'vg.Exp'. for exponential models. See examples for usage. +#' @param nugget (DxD)-matrix for the nugget effect. Default is a muted nugget (0). #' @param sill (DxD)-matrix for the partial sills of the correlation function #' @param anisRanges 2x2 or 3x3 matrix of ranges (see details) #' @param extraPar for certain correlation functions, extra parameters (smoothness, period, etc) #' #' @return an object of class "gmCgram" containing the linear model of corregionalization -#' of the nugget and the structure given. -#' @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` or `vg.Spherical`, they will all work and produce -#' a spherical model. The same applies for the following models: `vg.Gauss=0`; -#' `vg.Exp=vg.Exponential=2`. These constants are available after calling -#' `data("variogramModels")`. -#' No other model is currently available, but this data object will be +#' of the nugget and the structure given. +#' @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 +#' +#' Argument `anisRange` expects a matrix $M$ such that #' \deqn{ #' 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. -#' @family gmCgram +#' @family gmCgram #' @export -#' @aliases vg.Exp vg.exp vg.Exponential vg.Gau vg.gauss +#' @aliases vg.Exp vg.exp vg.Exponential vg.Gau vg.gauss #' vg.Gauss vg.Sph vg.sph vg.Spherical gsi.validModels #' #' @examples @@ -39,7 +41,7 @@ #' plot(vm) setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ utils::data("variogramModels") - stopifnot(all(dim(sill)==dim(nugget)), + stopifnot(all(dim(sill)==dim(nugget)), ncol(sill)==nrow(sill), type %in% gsi.validModels) dim(sill) = c(1,dim(sill)) @@ -57,18 +59,18 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ #' Subsetting of gmCgram variogram structures -#' +#' #' Extraction or combination of nested structures of a gmCgram object #' #' @param x `gmCgram` variogram object #' @param i indices of the structures that are desired to be kept (0=nugget) or removed (see details) #' @param ... extra arguments for generic functionality #' -#' @return a `gmCgram` variogram object with the desired structures only. -#' @details This function can be used to: extract the nugget (i=0), extract some +#' @return a `gmCgram` variogram object with the desired structures only. +#' @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 +#' 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. @@ -87,24 +89,24 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ nSo = length(i) nSi = dim(x$M)[1] if(i==0){# case only nugget wanted - out = with(x, list(type=type[i], - data=data[i], - nugget=nugget, - sill=structure(0*sill[1,,], dim=c(nSo, nD, nD)), + out = with(x, list(type=type[i], + data=data[i], + nugget=nugget, + sill=structure(0*sill[1,,], dim=c(nSo, nD, nD)), M=structure(M[1,,], dim=c(nSo, nG, nG)))) }else if( (0 %in% i) & !any(i<0)){ j = i[i>0] # case some structures wanted, including nugget - out = with(x, list(type=type[j], - data=data[j], - nugget=nugget, - sill=structure(sill[j,,], dim=c(nSo-1, nD, nD)), + out = with(x, list(type=type[j], + data=data[j], + nugget=nugget, + sill=structure(sill[j,,], dim=c(nSo-1, nD, nD)), M=structure(M[j,,], dim=c(nSo-1, nG, nG)))) }else if(all(i>0) | all(i<0)){# case some structured (un)wanted, nugget surely not wanted if(all(i<0)) nSo = nSi-nSo - out = with(x, list(type=type[i], - data=data[i], - nugget=nugget*0, - sill=structure(sill[i,,], dim=c(nSo, nD, nD)), + out = with(x, list(type=type[i], + data=data[i], + nugget=nugget*0, + sill=structure(sill[i,,], dim=c(nSo, nD, nD)), M=structure(M[i,,], dim=c(nSo, nG, nG)))) }else{# unsolvable case stop("index set i cannot merge 0 and negative numbers") @@ -115,7 +117,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ #' Combination of gmCgram variogram structures -#' +#' #' combination of nested structures of a gmCgram object #' #' @param x `gmCgram` variogram object @@ -130,7 +132,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ #' vm = v1+v2 "+.gmCgram" <- function(x,y) { y = as.gmCgram(y) - stopifnot(class(y)=="gmCgram", + stopifnot(class(y)=="gmCgram", dim(x$sill)[-1]==dim(y$sill)[-1], dim(x$M)[-1]==dim(y$M)[-1]) myfun = function(A,B){ @@ -155,7 +157,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ #' Subsetting of gmCgram variogram structures -#' +#' #' Extraction of some variables of a gmCgram object #' #' @param x \code{gmCgram} variogram object @@ -164,15 +166,15 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ #' is specified, \code{j} will be taken as equal to \code{i}!) #' @param ... extra arguments for generic functionality #' -#' @return a \code{gmCgram} variogram object with the desired variables only. -#' @details This function can be used to extract the model for a a subset of variables. -#' If only \code{i} is specified, \code{j} will be taken as equal to \code{i}. +#' @return a \code{gmCgram} variogram object with the desired variables only. +#' @details This function can be used to extract the model for a a subset of variables. +#' If only \code{i} is specified, \code{j} will be taken as equal to \code{i}. #' If you want to select all \code{i}'s for certain \code{j}'s or vice versa, give -#' \code{i=1:dim(x$nugget)[1]} and \code{j=} your desired indices, respectively +#' \code{i=1:dim(x$nugget)[1]} and \code{j=} your desired indices, respectively #' \code{j=1:dim(x$nugget)[2]} and \code{i=} your desired indices; replace \code{x} by the -#' object you are giving. If \code{i!=j}, the output will be a \code{c("gmXCgram","gmCgram")} +#' object you are giving. If \code{i!=j}, the output will be a \code{c("gmXCgram","gmCgram")} #' object, otherwise it will be a regular class \code{"gmCgram"} object. -#' If you want to extract "slots" or +#' 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. #' @export @@ -188,10 +190,10 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ nDi = length(i) nDj = length(j) nS = dim(x$M)[1] - out = with(x, list(type=type, - data=data, - nugget=structure(nugget[i,j, drop=F], dim=c(nDi, nDj)), - sill=structure(sill[,i,j, drop=F], dim=c(nS, nDi, nDj)), + out = with(x, list(type=type, + data=data, + nugget=structure(nugget[i,j, drop=F], dim=c(nDi, nDj)), + sill=structure(sill[,i,j, drop=F], dim=c(nS, nDi, nDj)), M=M)) class(out) = class(x) if(!all(i==j)) class(out) = unique(c("gmXCgram", class(out))) @@ -199,8 +201,8 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ } #' Length, and number of columns or rows -#' -#' Provide number of structures, and nr of variables of an LMC of class gmCgram +#' +#' Provide number of structures, and nr of variables of an LMC of class gmCgram #' #' @param x gmCgram object #' @@ -234,14 +236,14 @@ nrow.gmCgram = function(x) nrow(x$nugget) #' Convert a gmCgram object to an (evaluable) function -#' +#' #' Evaluate a gmCgram on some h values, or convert the gmCgram object into an evaluable function #' #' @param x a gmCgram object #' @param ... extra arguments for generic functionality #' #' @return a \code{function} that can be evaluated normally, with an argument \code{X} -#' and possibly another argument \code{Y}; both must have the same number of columns +#' and possibly another argument \code{Y}; both must have the same number of columns #' than the geographic dimension of the variogram (i.e. \code{dim(x$M)[3]}). #' @export #' @method as.function gmCgram @@ -280,55 +282,55 @@ predict.gmCgram = function(object, newdata, ...){ } -#' Convert theoretical structural functions to gmCgram format -#' -#' Convert covariance function or variogram models to the format gmCgram +#' Convert theoretical structural functions to gmCgram format +#' +#' Convert covariance function or variogram models to the format gmCgram #' of package gmGeostats -#' +#' #' @param m model to be converted #' @param ... further parameters #' -#' @return the covariance/variogram model, recasted to class \code{gmCgram}. +#' @return the covariance/variogram model, recasted to class \code{gmCgram}. #' This is a generic function. Methods exist for objects of class #' \code{LMCAnisCompo} (for compositional data) and \code{variogramModelList} -#' (as provided by package \code{gstat}). +#' (as provided by package \code{gstat}). #' @export #' @family gmCgram functions -as.gmCgram <- function(m, ...) UseMethod("as.gmCgram",m) +as.gmCgram <- function(m, ...) UseMethod("as.gmCgram",m) -#' @describeIn as.gmCgram Convert theoretical structural functions to gmCgram format +#' @describeIn as.gmCgram Convert theoretical structural functions to gmCgram format #' @method as.gmCgram default #' @export -as.gmCgram.default <- function(m,...) m +as.gmCgram.default <- function(m,...) m #' Draw cuves for covariance/variogram models -#' +#' #' Represent a gmCgram object as a matrix of lines in several plots #' #' @param x object to draw, of class gmCgram // curently only valid for symmetric functions #' @param xlim.up range of lag values to use in plots of the upper triangle #' @param xlim.lo range of lag values to use in plots of the lower triangle -#' @param vdir.up geograohic directions to represent in the upper triangle +#' @param vdir.up geograohic directions to represent in the upper triangle #' @param vdir.lo geograohic directions to represent in the lower triangle #' @param xlength number of discretization points to use for the curves (defaults to 200) #' @param varnames string vector, variable names to use in the labelling of axes #' @param add logical, should a new plot be created or stuff be added to an existing one? #' @param commonAxis logical, is a common Y axis for all plots in a row desired? -#' @param cov logical, should the covariance function (=TRUE) or the variogram (=FALSE) be plotted? -#' @param closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)? +#' @param cov logical, should the covariance function (=TRUE) or the variogram (=FALSE) be plotted? +#' @param closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)? #' defaults to TRUE #' @param ... further graphical parameters for the plotting function #' #' @return 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, +#' 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 \code{par(plot(x, closeplot=FALSE, ...))}, #' or \code{plot(x, closeplot=TRUE, ...)}. #' If you want to further add stuff to it, better just call \code{plot(x, closeplot=FALSE,...)}. The difference -#' is only relevant when working with the screen graphical device. -#' @export +#' is only relevant when working with the screen graphical device. +#' @export #' @method plot gmCgram #' @family gmCgram functions #' @examples @@ -338,7 +340,7 @@ as.gmCgram.default <- function(m,...) m #' vm = v1+v2 #' plot(vm) #' plot(vm, cov=FALSE) -plot.gmCgram = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= NULL, xlength=200, varnames = colnames(x$nugget), +plot.gmCgram = function(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, ...){ Dg = dim(x$M)[2] Ns = dim(x$M)[1] @@ -386,12 +388,12 @@ plot.gmCgram = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= N if(!is.null(xlim.lo)){ xseq.lo = seq(from=xlim.lo[1], to=xlim.lo[2], length.out=xlength) }else{ xseq.lo=NULL} - + opar = par() opar = par_remove_readonly(opar) - + if(closeplot) on.exit(par(opar)) - + getVdens = function(vdir, xseq){ if(is.null(vdir)|is.null(xseq)) return(NULL) Vdens = sapply(1:nrow(vdir), function(k){ @@ -407,16 +409,16 @@ plot.gmCgram = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= N dim(C0) = c(Dv,Dv) Vdens = sweep(-Vdens, c(1,3), C0, "+") ## this must be corrected when we allow non-symmetric covariances } - return(Vdens) + return(Vdens) } Vdens.up = getVdens(vdir.up, xseq.up) Vdens.lo = getVdens(vdir.lo, xseq.lo) - + myplot = function(...) matplot(type="l",ylab="", xlab="",xaxt="n", ...) if(add) myplot = function(...) matlines(...) if(!add){ par(mfrow=c(Dv+1,Dv+1), mar=c(2,3,0,0), oma=c(1,4,1,1), xpd=NA) - myplot(c(0,0), c(0,0), pch="", ann=FALSE, bty="n", yaxt="n") + myplot(c(0,0), c(0,0), pch="", ann=FALSE, bty="n", yaxt="n") } for(i in 1:Dv){ for(j in 1:Dv){ @@ -441,8 +443,8 @@ plot.gmCgram = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= N mtext(text=varnames[i], side = 1, line=3) mtext(text=varnames[i], side = 2, line=3) } - } - }} + } + }} mtext(text="lag distance", side=1, outer = TRUE, line=0) mtext(text=c("semivariogram","covariance")[cov+1], side=2, outer = TRUE, line=2) invisible(opar) @@ -452,11 +454,11 @@ plot.gmCgram = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= N #' Check for anisotropy of a theoretical variogram -#' +#' #' Checks for anisotropy of a theoretical variogram or covariance function model #' @param x variogram or covariance model object -#' @param tol tolerance for +#' @param tol tolerance for #' @param ... extra arguments for generic functionality #' #' @return Generic function. Returns of boolean answering the question of the name, @@ -511,7 +513,7 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){ #' Variogram method for gmSpatialModel objects -#' +#' #' Compute the empirical variogram of the conditioning data contained in a [gmSpatialModel-class] object #' #' @param object a gmSpatialModel object containing spatial data. @@ -519,7 +521,7 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){ #' @param ... further parameters to [gstat::variogram()] #' #' @return Currently the function is just a convenience wrapper on -#' the variogram calculation functionalities of package "gstat", +#' the variogram calculation functionalities of package "gstat", #' and returns objects of class "\code{gstatVariogram}". Check the #' help of \code{gstat::variogram} for further information. #' In the near future, methods will be created, which will depend on @@ -534,15 +536,15 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){ # Variogram calculations -# +# # Compute empricial variograms out of a spatial data object -# +# # @param object spatial data container # @param ... further parameters for variogram calculation -# -# @return depending on the input data, different kinds of empirical variograms +# +# @return depending on the input data, different kinds of empirical variograms # will be produced. See appropriate method descriptions. -# +# # @importFrom sp variogram # @export ##variogram <- function(object, ...) UseMethod("variogram", object) @@ -558,25 +560,25 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){ #' Empirical variogram or covariance function in 2D -#' -#' compute the empirical variogram or covariance function in a 2D case study -#' +#' +#' compute the empirical variogram or covariance function in a 2D case study +#' #' @param X matrix of Nx2 columns with the geographic coordinates #' @param Z matrix or data.frame of data with dimension (N,Dv) -#' @param Ff for variogram, matrix of basis functions with nrow(Ff)=N (can be a N-vector of 1s); +#' @param Ff for variogram, matrix of basis functions with nrow(Ff)=N (can be a N-vector of 1s); #' for covariance function, a (N,Dv)-matrix or a Dv-vector giving the mean values #' @param maxdist maximum lag distance to consider #' @param lagNr number of lags to consider -#' @param lags if maxdist and lagNr are not specified, either: (a) a matrix of 2 columns giving -#' minimal and maximal lag distance defining the lag classes to consider, or (b) a vector of lag breaks +#' @param lags if maxdist and lagNr are not specified, either: (a) a matrix of 2 columns giving +#' minimal and maximal lag distance defining the lag classes to consider, or (b) a vector of lag breaks #' @param azimuthNr number of azimuths to consider -#' @param azimuths if azimuthNr is not specified, either: (a) a matrix of 2 columns giving +#' @param azimuths if azimuthNr is not specified, either: (a) a matrix of 2 columns giving #' minimal and maximal azimuth defining the azimuth classes to consider, or (b) a vector of azimuth breaks #' @param maxbreadth maximal breadth (in lag units) orthogonal to the lag direction #' @param minpairs minimal number of pairs falling in each class to consider the calculation sufficient; defaults to 10 #' @param cov boolean, is covariance (TRUE) or variogram (FALSE) desired? defaults to variogram #' -#' @return An empirical variogram for the provided data. NOTE: avoid using directly gsi.* functions! They +#' @return 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 \code{\link{variogram}} instead. #' @export #' @family gmEVario functions @@ -594,7 +596,7 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){ #' vge["npairs",1] #' vge["lags",1] gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), - maxdist= max(dist(X[sample(nrow(X),min(nrow(X),1000)),]))/2, + 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){ @@ -611,7 +613,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), } # expand the information given into a set of columns stating conditions if(length(dim(lags))==0){ - lags = data.frame(minlag=lags[-length(lags)], maxlag=lags[-1]) + lags = data.frame(minlag=lags[-length(lags)], maxlag=lags[-1]) if(maxbreadth!=Inf) lags[,"maxbreadth"]=maxbreadth }else if(dim(lags)==2){ lags = data.frame(lags) @@ -621,21 +623,21 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), if(length(dim(azimuths))==0){ tol = (azimuths[2]-azimuths[1])/2 if(is.na(tol)) tol=180 - azimuths = data.frame(minaz=azimuths-tol, maxaz=azimuths+tol) + azimuths = data.frame(minaz=azimuths-tol, maxaz=azimuths+tol) }else if(dim(azimuths)==2){ azimuths = data.frame(azimuths) colnames(azimuths) = c("minaz","maxaz") }else stop("azimuths can be either a vector of lags or a data.frame, see ?gsi.EVario2D") - + # compute pairs of locations ij = expand.grid(1:nrow(X), 1:nrow(X))# indices XX = X[ij[,1],]-X[ij[,2],] # locations - XXabs = gmApply(XX, 1, function(x) sqrt(sum(x^2))) + XXabs = gmApply(XX, 1, function(x) sqrt(sum(x^2))) XXaz = gmApply(XX, 1, function(x) pi/2-atan2(x[2],x[1])) +2*pi XXaz = XXaz %% pi # residual to 180?? - + # compute residuals and pairs of variables appropriate to the structural function - if(cov){ + if(cov){ if(all(dim(Z)==dim(Ff))){ Z = as.matrix(Z-Ff) }else if(Dv==length(c(unlist(Ff)))){ @@ -648,8 +650,8 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), Z = lm(as.matrix(Z)~as.matrix(Ff)+0)$residuals ## ideally this should be a GLS fit ZZ = Z[ij[,1],]-Z[ij[,2],] } - - + + # output ## ATTENTION: needs to be changed to return a structure (3,Na)-matrix of objects, # like logratioVariogramAnisotropy @@ -667,7 +669,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), if(ncol(lags)>2){ tk_b = outer(xxabs * abs(sin((xxaz-(azs[i,2]-azs[i,1])))), lags[,3], "<=") tk_h = tk_h & tk_b - } + } n[,i] = colSums(tk_h) for(j in 1:Nh){ if(!is.na(n[j,i]) && n[j,i]>minpairs){ @@ -675,7 +677,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), vg[j,,,i] = gmApply((ZZ[tk_a,,][tk_h[,j],,]), c(2,3),"sum")/(n[j,i]) }else{ # semi-variogram aux = rowSums( - gmApply(X=ZZ[tk_a,][tk_h[,j],], 1, function(x)outer(x,x,"*")) + gmApply(X=ZZ[tk_a,][tk_h[,j],], 1, function(x)outer(x,x,"*")) )/(2*n[j,i]) vg[j,,,i] = aux } @@ -685,7 +687,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), } return(list(gamma=vg[,,,i], lags=gsi.lagClass(lags), npairs =n[,i])) }) - + # output attr(res, "directions") = gsi.azimuthInterval(azimuths) # attr(res, "lags") = gsi.lagClass(lags) @@ -695,34 +697,34 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), } #' Empirical variogram or covariance function in 3D -#' -#' compute the empirical variogram or covariance function in a 3D case study -#' +#' +#' compute the empirical variogram or covariance function in a 3D case study +#' #' @param X matrix of Nx3 columns with the geographic coordinates #' @param Z matrix or data.frame of data with dimension (N,Dv) -#' @param Ff for variogram, matrix of basis functions with nrow(Ff)=N -#' (can be a N-vector of 1s; should include the vector of 1s!!); +#' @param Ff for variogram, matrix of basis functions with nrow(Ff)=N +#' (can be a N-vector of 1s; should include the vector of 1s!!); #' for covariance function, a (N,Dv)-matrix or a Dv-vector giving the mean values #' @param maxdist maximum lag distance to consider #' @param lagNr number of lags to consider -#' @param lags if maxdist and lagNr are not specified, either: (a) a matrix of 2 columns giving -#' minimal and maximal lag distance defining the lag classes to consider, or (b) a vector of lag breaks -#' @param dirvecs matrix which rows are the director vectors along which variograms will be computed (these will be normalized!) -#' @param angtol scalar, angular tolerance applied (in degrees; defaults to 90??, ie. isotropic) +#' @param lags if maxdist and lagNr are not specified, either: (a) a matrix of 2 columns giving +#' minimal and maximal lag distance defining the lag classes to consider, or (b) a vector of lag breaks +#' @param dirvecs matrix which rows are the director vectors along which variograms will be computed (these will be normalized!) +#' @param angtol scalar, angular tolerance applied (in degrees; defaults to 90??, ie. isotropic) #' @param maxbreadth maximal breadth (in lag units) orthogonal to the lag direction (defaults to `Inf`, ie. not used) #' @param minpairs minimal number of pairs falling in each class to consider the calculation sufficient; defaults to 10 #' @param cov boolean, is covariance (TRUE) or variogram (FALSE) desired? defaults to variogram #' -#' @return An empirical variogram for the provided data. NOTE: avoid using directly gsi.* functions! They +#' @return 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 \code{\link{variogram}} instead. -#' -#' +#' +#' #' @export #' @family gmEVario functions #' # @examples # dt <- readr::read_csv("~/Geochem_sum_imp2.csv") -# +# # X = as.matrix(dt[,c("X","Y","Z")]) # Z = as.matrix(compositions::clr(dt[,c("Cu","Zn","Pb", "As", "Cd", "In", "Other")])) # dirvecs = rbind( c(1,0,0), c(1,1,0), c(0,1,0), c(-1,1,0), c(0,0,1)) @@ -736,11 +738,11 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), # vge["lags",1] # plot.gmEVario(vge, varnames = colnames(Z), commonAxis=FALSE, # vdir.up = 1:4, vdir.lo=5, xlim.up=c(0,50), xlim.lo=c(0,15) -# ) +# ) ## ie. the 4 planar directions on the upper triangle, ## the downhole direction in the lower triangle gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), - maxdist= max(dist(X[sample(nrow(X),min(nrow(X),1000)),]))/2, + 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){ @@ -757,15 +759,15 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ## expand the information given into a set of columns stating conditions # prepare lags if(length(dim(lags))==0){ - lags = data.frame(minlag=lags[-length(lags)], maxlag=lags[-1]) - # if(maxbreadth!=Inf) + lags = data.frame(minlag=lags[-length(lags)], maxlag=lags[-1]) + # if(maxbreadth!=Inf) lags[,"maxbreadth"] = maxbreadth }else if(dim(lags)==2){ lags = data.frame(lags) colnames(lags) = c("minlag","maxlag","maxbreadth")[1:ncol(lags)] }else stop("lags can be either a vector of lags or a data.frame, see ?gsi.EVario3D") lagsSq = lags^2 - + # prepare directions with tolerance if(length(dim(dirvecs))==0){ if(length(dirvecs)!=3) stop("dirvecs must be a matrix with 3 columns!") @@ -775,15 +777,15 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), if(length(angtol) %in% c(1, nrow(dirvecs))){ dirvecs = cbind(dirvecs, tol=cos(angtol*pi/180) ) }else stop("angtol can be either a scalar or a vector of length=nrow(dirvecs), see ?gsi.EVario3D") - + # compute pairs of locations ij = expand.grid(1:nrow(X), 1:nrow(X))# indices XX = X[ij[,1],]-X[ij[,2],] # locations XXmodSq = gmApply(XX, 1, function(x) sum(x^2)) XXmod = sqrt(XXmodSq) - + # compute residuals and pairs of variables appropriate to the structural function - if(cov){ + if(cov){ if(all(dim(Z)==dim(Ff))){ Z = as.matrix(Z-Ff) }else if(Dv==length(c(unlist(Ff)))){ @@ -796,7 +798,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), Z = lm(as.matrix(Z)~as.matrix(Ff)+0)$residuals ## ideally this should be a GLS fit ZZ = Z[ij[,1],]-Z[ij[,2],] } - + # output ## ATTENTION: needs to be changed to return a structure (3,Na)-matrix of objects, # like logratioVariogramAnisotropy @@ -816,7 +818,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), residualsSq = XXmodSq[tk_a]-projections[tk_a]^2 tk_b = outer(residualsSq, lagsSq[,3], "<=") tk_h = tk_h & tk_b - } + } n[,i] = colSums(tk_h) for(j in 1:Nh){ if(!is.na(n[j,i]) && n[j,i]>minpairs){ @@ -824,7 +826,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), vg[j,,,i] = gmApply((ZZ[tk_a,,][tk_h[,j],,]), c(2,3),"sum")/(n[j,i]) }else{ # semi-variogram aux = rowSums( - gmApply(X=ZZ[tk_a,][tk_h[,j],], 1, function(x)outer(x,x,"*")) + gmApply(X=ZZ[tk_a,][tk_h[,j],], 1, function(x)outer(x,x,"*")) )/(2*n[j,i]) vg[j,,,i] = aux } @@ -834,7 +836,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), } return(list(gamma=vg[,,,i], lags=gmGeostats:::gsi.lagClass(lags), npairs =n[,i])) }) - + # output attr(res, "directions") = gmGeostats:::gsi.directorVector(dirvecs[,1:3]) # attr(res, "lags") = gsi.lagClass(lags) @@ -848,13 +850,13 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), #' Plot empirical variograms -#' +#' #' Flexible plot of an empirical variogram of class gmEVario -#' +#' #' @param x object to print, of class gmEVario #' @param xlim.up range of X values to be used for the diagrams of the upper triangle #' @param xlim.lo range of X values to be used for the diagrams of the lower triangle -#' @param vdir.up in case of anisotropic variograms, indices of the directions to be plotted +#' @param vdir.up in case of anisotropic variograms, indices of the directions to be plotted #' on the upper triangle #' @param vdir.lo ..., indices of the directions to be plotted on the lower triangle #' @param varnames variable names to be used @@ -862,16 +864,16 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), #' @param add boolean, add stuff to an existing diagram? #' @param commonAxis boolean, should vertical axes be shared by all plots in a row? #' @param cov boolean, is this a covariance? (if FALSE, it is a variogram) -#' @param closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)? +#' @param closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)? #' defaults to TRUE #' @param ... further parameters to \code{\link{matplot}} #' -#' @return invisibly, the graphical parameters active before calling the function. +#' @return 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 \code{x} has been +#' +#' How to use arguments `vdir.lo` and `vdir.up`? Each empirical variogram \code{x} has been #' computed along certain distances, recorded in its attributes and retrievable with command -#' \code{\link{ndirections}}. +#' \code{\link{ndirections}}. #' @export #' @family gmEVario functions #' @method plot gmEVario @@ -884,8 +886,8 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), #' vge = gsi.EVario2D(X,Z) #' plot(vge) #' plot(vge, pch=22, lty=1, bg="grey") -plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= NULL, - varnames = dimnames(x$gamma)[[2]], type="o", +plot.gmEVario = function(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, ...){ Dv = dim(x[1,1][[1]])[2] @@ -912,12 +914,12 @@ plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= xlim.lo = c(0, maxdist) } } - + opar = par() opar = par_remove_readonly(opar) - + if(closeplot) on.exit(par(opar)) - + myplot = function(...) matplot(type=type, ylab="", xlab="",xaxt="n", ...) if(add) myplot = function(...) matpoints(type=type, ...) if(!add){ @@ -931,7 +933,7 @@ plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= ylim = range(sapply(vdir.lo, function(kk) x["gamma",kk][[1]][,i,j]), na.rm = TRUE) if(commonAxis) ylim=range(sapply(vdir.lo, function(kk) x["gamma",kk][[1]][,,j]), na.rm=TRUE) myplot( - sapply(vdir.lo, function(kk) gsi.midValues.lagClass(x["lags",kk][[1]])), + sapply(vdir.lo, function(kk) gsi.midValues.lagClass(x["lags",kk][[1]])), sapply(vdir.lo, function(kk) x["gamma",kk][[1]][,i,j]), ylim=ylim, ...) if(i==j){ axis(side = 3) @@ -944,15 +946,15 @@ plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= ylim = range(sapply(vdir.up, function(kk) x["gamma",kk][[1]][,i,j]), na.rm = TRUE) if(commonAxis) ylim=range(sapply(vdir.up, function(kk) x["gamma",kk][[1]][,,j]), na.rm=TRUE) myplot( - sapply(vdir.up, function(kk) gsi.midValues.lagClass(x["lags",kk][[1]])), + sapply(vdir.up, function(kk) gsi.midValues.lagClass(x["lags",kk][[1]])), sapply(vdir.up, function(kk) x["gamma",kk][[1]][,i,j]), ylim=ylim, ...) if(i==j){ axis(side = 1) mtext(text=varnames[i], side = 1, line=3) mtext(text=varnames[i], side = 2, line=3) } - } - }} + } + }} mtext(text="lag distance", side=1, outer = TRUE, line=0) mtext(text=c("semivariogram","covariance")[cov+1], side=2, outer = TRUE, line=2) attr(opar, "vdir.up")=vdir.up @@ -965,19 +967,19 @@ plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= -#' Convert empirical structural function to gmEVario format -#' -#' Convert empirical covariance functions or variograms to the format gmEVario +#' Convert empirical structural function to gmEVario format +#' +#' Convert empirical covariance functions or variograms to the format gmEVario #' of package gmGeostats -#' +#' #' @param vgemp variogram/covariance function to be converted #' @param ... further parameters #' -#' @return the empirical covariance function or variogram, recasted to class -#' \code{gmEVario}. This is a generic function. Methods exist for objects of -#' class \code{logratioVariogram}\code{logratioVariogramAnisotropy} +#' @return the empirical covariance function or variogram, recasted to class +#' \code{gmEVario}. This is a generic function. Methods exist for objects of +#' class \code{logratioVariogram}\code{logratioVariogramAnisotropy} #' (for compositional data) and \code{gstatVariogram} -#' (from package \code{gstat}). +#' (from package \code{gstat}). #' @export #' @aliases as.gmEVario.gstatVariogram as.gmEVario.logratioVariogram #' as.gmEVario.logratioVariogramAnisotropy @@ -995,26 +997,26 @@ variogramModelPlot <- function(vg, ...) UseMethod("variogramModelPlot", vg) #' Quick plotting of empirical and theoretical variograms -#' Quick and dirty plotting of empirical variograms/covariances with or without their models +#' Quick and dirty plotting of empirical variograms/covariances with or without their models #' @param vg empirical variogram or covariance function #' @param model optional, theoretical variogram or covariance function #' @param col colors to use for the several directional variograms #' @param commonAxis boolean, should all plots in a row share the same vertical axis? -#' @param newfig boolean, should a new figure be created? otherwise user should ensure the device space is appropriately managed -#' @param closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)? +#' @param newfig boolean, should a new figure be created? otherwise user should ensure the device space is appropriately managed +#' @param closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)? #' defaults to TRUE #' @param ... further parameters to underlying plot or matplot functions #' -#' @return 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 +#' @return 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 \code{\link{par}}, e.g. \code{par(variogramModelPlot(...))}; if you #' want to leave the plot open for further changes give the extra argument `closeplot=FALSE`. #' @export -#' @family variogramModelPlot -#' @family gmEVario functions -#' @family gmCgram functions +#' @family variogramModelPlot +#' @family gmEVario functions +#' @family gmCgram functions #' @seealso [logratioVariogram()] #' @method variogramModelPlot gmEVario #' @@ -1030,15 +1032,15 @@ variogramModelPlot <- function(vg, ...) UseMethod("variogramModelPlot", vg) #' Z = as.matrix(jura.pred[,c("Zn","Cd","Pb")]) #' vge = gsi.EVario2D(X,Z) #' variogramModelPlot(vge, vm) -#' -#' +#' +#' variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogramModelList object containing a variogram model fitted to vg col = rev(rainbow(ndirections(vg))), commonAxis = FALSE, newfig = TRUE, closeplot=TRUE, ...){ opar = plot(vg, commonAxis=commonAxis, add=!newfig, col=col, closeplot=is.null(model), ...) opar = par_remove_readonly(opar) - + if(closeplot) on.exit(par(opar)) vdir.lo = attr(opar, "vdir.lo") vdir.up = attr(opar, "vdir.up") @@ -1046,8 +1048,8 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra return(invisible(opar)) # OTHERWISE: add the curves for the model aux = as.directorVector(attr(vg, "directions")) - if(!is.null(vdir.lo)) vdir.lo = aux[vdir.lo,] - if(!is.null(vdir.up)) vdir.up = aux[vdir.up,] + if(!is.null(vdir.lo)) vdir.lo = aux[vdir.lo,] + if(!is.null(vdir.up)) vdir.up = aux[vdir.up,] opar = plot(as.gmCgram(model), vdir.up= vdir.up, vdir.lo= vdir.lo, add=TRUE, cov =FALSE, ...) #f = as.function(as.gmCgram(gg)) #Dv = dim(vg$gamma)[2] @@ -1062,7 +1064,7 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra # hd = seq(from=xlim[1], to=xlim[2], length.out = 200) # Y = outer(hd, 1:nrow(dirs.lo), function(h,k){ # f(rep(0,Dg), dirs.lo[k,]*h)[i,j] - # }) + # }) # matlines(hd, Y, col=col[vdir.lo], ...) # } # if((i<=j)&(!is.null(vdir.up))){ @@ -1072,26 +1074,26 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra # hd = seq(from=xlim[1], to=xlim[2], length.out = 200) # Y = outer(hd, 1:nrow(dirs.up), function(h,k){ # f(rep(0,Dg), dirs.up[k,]*h)[i,j] - # }) + # }) # matlines(hd, Y, col=col[vdir.up], ...) - # } - # }} + # } + # }} invisible(opar) -} +} #' Number of directions of an empirical variogram -#' +#' #' Returns the number of directions at which an empirical variogram was computed -#' +#' #' @param x empirical variogram object #' -#' @return Generic function. It provides the +#' @return Generic function. It provides the #' number of directions at which an empirical variogram was computed #' @export -#' @family gmEVario functions -#' @family gmCgram functions +#' @family gmEVario functions +#' @family gmCgram functions #' @seealso [logratioVariogram()], [gstat::variogram()] ndirections <- function(x){ UseMethod("ndirections", x) } @@ -1107,7 +1109,7 @@ ndirections.azimuth = function(x) length(x) #' @method ndirections azimuthInterval #' @export ndirections.azimuthInterval = function(x) length(x[[1]]) -#' @describeIn ndirections method for empirical logratio variograms with anisotropy +#' @describeIn ndirections method for empirical logratio variograms with anisotropy #' @method ndirections logratioVariogramAnisotropy #' @export ndirections.logratioVariogramAnisotropy = function(x) ndirections(attr(x,"directions")) @@ -1124,7 +1126,7 @@ ndirections.gmEVario = function(x) ndirections(attr(x,"directions")) #' @export ndirections.gstatVariogram = function(x){ length(unique(paste(x$dir.hor, x$dir.ver))) -} +} @@ -1135,12 +1137,12 @@ ndirections.gstatVariogram = function(x){ gsi.azimuth = function(x){ class(x) = c("azimuth","directionClass") return(x) -} +} gsi.azimuthInterval = function(x){ class(x) = c("azimuthInterval","directionClass") return(x) -} +} @@ -1155,14 +1157,14 @@ print.directionClass = function(x, complete=TRUE, ...){ if(complete){ print(unclass(x), ...) } -} +} as.directorVector <- function(x){ UseMethod("as.directorVector",x) } -#' @method as.directorVector default +#' @method as.directorVector default as.directorVector.default = function(x,...) x #' @method as.directorVector azimuth @@ -1185,7 +1187,7 @@ as.directorVector.azimuthInterval = function(x, D=2){ gsi.lagdists = function(x){ class(x) = c("lagdist","lagClass") return(x) -} +} gsi.lagClass = function(x){ if(length(dim(x))!=2) stop("provided lag classes are not matrix-like!") @@ -1209,7 +1211,7 @@ gsi.midValues.default = function(x) x #' @method gsi.midValues lagClass gsi.midValues.lagClass <- function(x){ (x[[1]]+x[[2]])/2 } - + #' @method gsi.midValues azimuthIntervals gsi.midValues.azimuthInterval <- function(x){ (x[[1]]+x[[2]])/2 } @@ -1218,49 +1220,49 @@ gsi.midValues.azimuthInterval <- function(x){ (x[[1]]+x[[2]])/2 } ## theoretical structural functions -# S3 -> S4 classes +# S3 -> S4 classes # cat("creating variogram model classes\n") # abstract classes #' @title Structural function model specification -#' @description Abstract class, containing any specification of a variogram (or covariance) model. +#' @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 "gmCgram" (see [setCgram()] for an example), and (possibly) coercion to #' class "variogramModel" or "variogramModelList" (see [gstat::vgm()]) #' @export #' @include compositionsCompatibility.R #' @include gstatCompatibility.R #' @include preparations.R -setClassUnion(name="ModelStructuralFunctionSpecification", +setClassUnion(name="ModelStructuralFunctionSpecification", members=c("NULL","gmCgram", "LMCAnisCompo", "variogramModelList", "variogramModel")) -# +# # #### container class -------------- # # An S4 class to represent a Gaussian random field specification # # -# # @slot structure ModelStructuralFunctionSpecification. Variogram or +# # @slot structure ModelStructuralFunctionSpecification. Variogram or # # (generalised) covariance function specification, typically an object # # obtained from a call to functions such as \code{\link{setCgram}}, # # \code{\link{LMCAnisCompo}} or \code{gstat::vgm}. # # @slot formula formula specifying the structure -# # of dependence of the mean of the random field w.r.to spatial coordinates -# # and/or covariables; typically it will have no left-hand-side term; +# # of dependence of the mean of the random field w.r.to spatial coordinates +# # and/or covariables; typically it will have no left-hand-side term; # # @slot beta numeric, a vector with as many coefficients as terms the formula # # above requires for a full specification of the trend; if unknown, these can -# # be NAs, as many as needed. +# # be NAs, as many as needed. # # # # @return A object with the slots populated as given # # @export # # @seealso [gmSpatialModel-class], and the `make.gm*` functions referenced there -# setClass("gmGaussianModel", +# setClass("gmGaussianModel", # slots = list(structure = "ModelStructuralFunctionSpecification", # formula="formula", # beta = "structure") # ) -# +# # #setMethod("initialize", signature="gmGaussianModel", # # def=function(.Object, structure, formula, beta){ # # .Object@formula = formula @@ -1274,21 +1276,21 @@ setClassUnion(name="ModelStructuralFunctionSpecification", #### container class -------------- # An S4 class to represent a Gaussian random field specification # -# @slot structure ModelStructuralFunctionSpecification. Variogram or +# @slot structure ModelStructuralFunctionSpecification. Variogram or # (generalised) covariance function specification, typically an object # obtained from a call to functions such as \code{\link{setCgram}}, # \code{\link{LMCAnisCompo}} or \code{gstat::vgm}. # @slot formula formula specifying the structure -# of dependence of the mean of the random field w.r.to spatial coordinates -# and/or covariables; typically it will have no left-hand-side term; +# of dependence of the mean of the random field w.r.to spatial coordinates +# and/or covariables; typically it will have no left-hand-side term; # @slot beta numeric, a vector with as many coefficients as terms the formula # above requires for a full specification of the trend; if unknown, these can -# be NAs, as many as needed. +# be NAs, as many as needed. # # @return A object with the slots populated as given # @export # @seealso [gmSpatialModel-class], and the `make.gm*` functions referenced there -setClass("gmGaussianModel", +setClass("gmGaussianModel", slots = list(structure = "ModelStructuralFunctionSpecification", formula="formula", beta = "numeric") @@ -1302,9 +1304,9 @@ setClass("gmGaussianModel", # abstract classes #' @title Empirical structural function specification -#' @description Abstract class, containing any specification of an empirical variogram +#' @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 "gmEVario" (see [gsi.EVario2D()] for an example), and (possibly) coercion to #' class "gstatVariogram" (see [gstat::variogram()]) #' @export #' @include compositionsCompatibility.R -- GitLab