Skip to content
Snippets Groups Projects
Commit c6149a3e authored by Pospiech, Solveig (FWGB) - 133453's avatar Pospiech, Solveig (FWGB) - 133453
Browse files

setCgram: small improvements in the help-page

parent 26ccbe0f
No related branches found
No related tags found
No related merge requests found
#### theoretical variogram --------------------- #### theoretical variogram ---------------------
#' gmGeostats Variogram models #' Generate D-variate variogram models
#' set up a D-variate variogram models
#' #'
#' @param type constant indicating the model of correlation function #' @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 nugget (DxD)-matrix for the nugget effect #'
#' @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 sill (DxD)-matrix for the partial sills of the correlation function
#' @param anisRanges 2x2 or 3x3 matrix of ranges (see details) #' @param anisRanges 2x2 or 3x3 matrix of ranges (see details)
#' @param extraPar for certain correlation functions, extra parameters (smoothness, period, etc) #' @param extraPar for certain correlation functions, extra parameters (smoothness, period, etc)
#' #'
#' @return an object of class "gmCgram" containing the linear model of corregionalization #' @return an object of class "gmCgram" containing the linear model of corregionalization
#' of the nugget and the structure given. #' of the nugget and the structure given.
#' @details The argument `type` must be an integer indicating the model to be used. #' @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 #' 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 #' 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=0`; #' a spherical model. The same applies for the following models:
#' `vg.Exp=vg.Exponential=2`. These constants are available after calling #' `vg.Gauss = vg.Gau = vg.gau = 0`;
#' `data("variogramModels")`. #' `vg.Exponential = vg.Exp = vg. exp = 2`.
#' No other model is currently available, but this data object will be #' These constants are available after calling `data("variogramModels")`.
#' No other model is currently available, but this data object will be
#' regularly updated. #' regularly updated.
#' The constant vector `gsi.validModels` contains all currently valid models. #' 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{ #' \deqn{
#' h^2 = (\mathbf{x}_i-\mathbf{x}_j)\cdot M^{-1}\cdot (\mathbf{x}_i-\mathbf{x}_j)^t #' 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. #' is the (square of) the lag distance to be fed into the correlation function.
#' @family gmCgram #' @family gmCgram
#' @export #' @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 #' vg.Gauss vg.Sph vg.sph vg.Spherical gsi.validModels
#' #'
#' @examples #' @examples
...@@ -39,7 +41,7 @@ ...@@ -39,7 +41,7 @@
#' plot(vm) #' plot(vm)
setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
utils::data("variogramModels") utils::data("variogramModels")
stopifnot(all(dim(sill)==dim(nugget)), stopifnot(all(dim(sill)==dim(nugget)),
ncol(sill)==nrow(sill), ncol(sill)==nrow(sill),
type %in% gsi.validModels) type %in% gsi.validModels)
dim(sill) = c(1,dim(sill)) dim(sill) = c(1,dim(sill))
...@@ -57,18 +59,18 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ ...@@ -57,18 +59,18 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
#' Subsetting of gmCgram variogram structures #' Subsetting of gmCgram variogram structures
#' #'
#' Extraction or combination of nested structures of a gmCgram object #' Extraction or combination of nested structures of a gmCgram object
#' #'
#' @param x `gmCgram` variogram 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 i indices of the structures that are desired to be kept (0=nugget) or removed (see details)
#' @param ... extra arguments for generic functionality #' @param ... extra arguments for generic functionality
#' #'
#' @return a `gmCgram` variogram object with the desired structures only. #' @return a `gmCgram` variogram object with the desired structures only.
#' @details This function can be used to: extract the nugget (i=0), extract some #' @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), #' 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; #' 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 #' 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 #' "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) #' variogram matrices, use the `[`-notation. The contrary operation (adding structures together)
#' is obtained by summing (+) two `gmCgram` objects. #' is obtained by summing (+) two `gmCgram` objects.
...@@ -87,24 +89,24 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ ...@@ -87,24 +89,24 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
nSo = length(i) nSo = length(i)
nSi = dim(x$M)[1] nSi = dim(x$M)[1]
if(i==0){# case only nugget wanted if(i==0){# case only nugget wanted
out = with(x, list(type=type[i], out = with(x, list(type=type[i],
data=data[i], data=data[i],
nugget=nugget, nugget=nugget,
sill=structure(0*sill[1,,], dim=c(nSo, nD, nD)), sill=structure(0*sill[1,,], dim=c(nSo, nD, nD)),
M=structure(M[1,,], dim=c(nSo, nG, nG)))) M=structure(M[1,,], dim=c(nSo, nG, nG))))
}else if( (0 %in% i) & !any(i<0)){ }else if( (0 %in% i) & !any(i<0)){
j = i[i>0] # case some structures wanted, including nugget j = i[i>0] # case some structures wanted, including nugget
out = with(x, list(type=type[j], out = with(x, list(type=type[j],
data=data[j], data=data[j],
nugget=nugget, nugget=nugget,
sill=structure(sill[j,,], dim=c(nSo-1, nD, nD)), sill=structure(sill[j,,], dim=c(nSo-1, nD, nD)),
M=structure(M[j,,], dim=c(nSo-1, nG, nG)))) 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 }else if(all(i>0) | all(i<0)){# case some structured (un)wanted, nugget surely not wanted
if(all(i<0)) nSo = nSi-nSo if(all(i<0)) nSo = nSi-nSo
out = with(x, list(type=type[i], out = with(x, list(type=type[i],
data=data[i], data=data[i],
nugget=nugget*0, nugget=nugget*0,
sill=structure(sill[i,,], dim=c(nSo, nD, nD)), sill=structure(sill[i,,], dim=c(nSo, nD, nD)),
M=structure(M[i,,], dim=c(nSo, nG, nG)))) M=structure(M[i,,], dim=c(nSo, nG, nG))))
}else{# unsolvable case }else{# unsolvable case
stop("index set i cannot merge 0 and negative numbers") stop("index set i cannot merge 0 and negative numbers")
...@@ -115,7 +117,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ ...@@ -115,7 +117,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
#' Combination of gmCgram variogram structures #' Combination of gmCgram variogram structures
#' #'
#' combination of nested structures of a gmCgram object #' combination of nested structures of a gmCgram object
#' #'
#' @param x `gmCgram` variogram object #' @param x `gmCgram` variogram object
...@@ -130,7 +132,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ ...@@ -130,7 +132,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
#' vm = v1+v2 #' vm = v1+v2
"+.gmCgram" <- function(x,y) { "+.gmCgram" <- function(x,y) {
y = as.gmCgram(y) y = as.gmCgram(y)
stopifnot(class(y)=="gmCgram", stopifnot(class(y)=="gmCgram",
dim(x$sill)[-1]==dim(y$sill)[-1], dim(x$sill)[-1]==dim(y$sill)[-1],
dim(x$M)[-1]==dim(y$M)[-1]) dim(x$M)[-1]==dim(y$M)[-1])
myfun = function(A,B){ myfun = function(A,B){
...@@ -155,7 +157,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ ...@@ -155,7 +157,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
#' Subsetting of gmCgram variogram structures #' Subsetting of gmCgram variogram structures
#' #'
#' Extraction of some variables of a gmCgram object #' Extraction of some variables of a gmCgram object
#' #'
#' @param x \code{gmCgram} variogram object #' @param x \code{gmCgram} variogram object
...@@ -164,15 +166,15 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ ...@@ -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}!) #' is specified, \code{j} will be taken as equal to \code{i}!)
#' @param ... extra arguments for generic functionality #' @param ... extra arguments for generic functionality
#' #'
#' @return a \code{gmCgram} variogram object with the desired variables only. #' @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. #' @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 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 #' 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 #' \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. #' 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 #' "elements" of the variogram, use the $-notation. If you want to extract variables of the
#' variogram matrices, use the `[`-notation. #' variogram matrices, use the `[`-notation.
#' @export #' @export
...@@ -188,10 +190,10 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){ ...@@ -188,10 +190,10 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
nDi = length(i) nDi = length(i)
nDj = length(j) nDj = length(j)
nS = dim(x$M)[1] nS = dim(x$M)[1]
out = with(x, list(type=type, out = with(x, list(type=type,
data=data, data=data,
nugget=structure(nugget[i,j, drop=F], dim=c(nDi, nDj)), nugget=structure(nugget[i,j, drop=F], dim=c(nDi, nDj)),
sill=structure(sill[,i,j, drop=F], dim=c(nS, nDi, nDj)), sill=structure(sill[,i,j, drop=F], dim=c(nS, nDi, nDj)),
M=M)) M=M))
class(out) = class(x) class(out) = class(x)
if(!all(i==j)) class(out) = unique(c("gmXCgram", class(out))) 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){ ...@@ -199,8 +201,8 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
} }
#' Length, and number of columns or rows #' 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 #' @param x gmCgram object
#' #'
...@@ -234,14 +236,14 @@ nrow.gmCgram = function(x) nrow(x$nugget) ...@@ -234,14 +236,14 @@ nrow.gmCgram = function(x) nrow(x$nugget)
#' Convert a gmCgram object to an (evaluable) function #' Convert a gmCgram object to an (evaluable) function
#' #'
#' Evaluate a gmCgram on some h values, or convert the gmCgram object into 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 x a gmCgram object
#' @param ... extra arguments for generic functionality #' @param ... extra arguments for generic functionality
#' #'
#' @return a \code{function} that can be evaluated normally, with an argument \code{X} #' @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]}). #' than the geographic dimension of the variogram (i.e. \code{dim(x$M)[3]}).
#' @export #' @export
#' @method as.function gmCgram #' @method as.function gmCgram
...@@ -280,55 +282,55 @@ predict.gmCgram = function(object, newdata, ...){ ...@@ -280,55 +282,55 @@ predict.gmCgram = function(object, newdata, ...){
} }
#' Convert theoretical structural functions to gmCgram format #' Convert theoretical structural functions to gmCgram format
#' #'
#' Convert covariance function or variogram models to the format gmCgram #' Convert covariance function or variogram models to the format gmCgram
#' of package gmGeostats #' of package gmGeostats
#' #'
#' @param m model to be converted #' @param m model to be converted
#' @param ... further parameters #' @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 #' This is a generic function. Methods exist for objects of class
#' \code{LMCAnisCompo} (for compositional data) and \code{variogramModelList} #' \code{LMCAnisCompo} (for compositional data) and \code{variogramModelList}
#' (as provided by package \code{gstat}). #' (as provided by package \code{gstat}).
#' @export #' @export
#' @family gmCgram functions #' @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 #' @method as.gmCgram default
#' @export #' @export
as.gmCgram.default <- function(m,...) m as.gmCgram.default <- function(m,...) m
#' Draw cuves for covariance/variogram models #' Draw cuves for covariance/variogram models
#' #'
#' Represent a gmCgram object as a matrix of lines in several plots #' 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 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.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 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 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 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 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 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 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 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 closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)?
#' defaults to TRUE #' defaults to TRUE
#' @param ... further graphical parameters for the plotting function #' @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 #' @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 #' 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, ...))}, #' 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, ...)}. #' 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 #' 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. #' is only relevant when working with the screen graphical device.
#' @export #' @export
#' @method plot gmCgram #' @method plot gmCgram
#' @family gmCgram functions #' @family gmCgram functions
#' @examples #' @examples
...@@ -338,7 +340,7 @@ as.gmCgram.default <- function(m,...) m ...@@ -338,7 +340,7 @@ as.gmCgram.default <- function(m,...) m
#' vm = v1+v2 #' vm = v1+v2
#' plot(vm) #' plot(vm)
#' plot(vm, cov=FALSE) #' 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, ...){ add=FALSE, commonAxis=TRUE, cov =TRUE, closeplot=TRUE, ...){
Dg = dim(x$M)[2] Dg = dim(x$M)[2]
Ns = dim(x$M)[1] 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 ...@@ -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)){ if(!is.null(xlim.lo)){
xseq.lo = seq(from=xlim.lo[1], to=xlim.lo[2], length.out=xlength) xseq.lo = seq(from=xlim.lo[1], to=xlim.lo[2], length.out=xlength)
}else{ xseq.lo=NULL} }else{ xseq.lo=NULL}
opar = par() opar = par()
opar = par_remove_readonly(opar) opar = par_remove_readonly(opar)
if(closeplot) on.exit(par(opar)) if(closeplot) on.exit(par(opar))
getVdens = function(vdir, xseq){ getVdens = function(vdir, xseq){
if(is.null(vdir)|is.null(xseq)) return(NULL) if(is.null(vdir)|is.null(xseq)) return(NULL)
Vdens = sapply(1:nrow(vdir), function(k){ 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 ...@@ -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) dim(C0) = c(Dv,Dv)
Vdens = sweep(-Vdens, c(1,3), C0, "+") ## this must be corrected when we allow non-symmetric covariances 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.up = getVdens(vdir.up, xseq.up)
Vdens.lo = getVdens(vdir.lo, xseq.lo) Vdens.lo = getVdens(vdir.lo, xseq.lo)
myplot = function(...) matplot(type="l",ylab="", xlab="",xaxt="n", ...) myplot = function(...) matplot(type="l",ylab="", xlab="",xaxt="n", ...)
if(add) myplot = function(...) matlines(...) if(add) myplot = function(...) matlines(...)
if(!add){ if(!add){
par(mfrow=c(Dv+1,Dv+1), mar=c(2,3,0,0), oma=c(1,4,1,1), xpd=NA) 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(i in 1:Dv){
for(j 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 ...@@ -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 = 1, line=3)
mtext(text=varnames[i], side = 2, line=3) mtext(text=varnames[i], side = 2, line=3)
} }
} }
}} }}
mtext(text="lag distance", side=1, outer = TRUE, line=0) mtext(text="lag distance", side=1, outer = TRUE, line=0)
mtext(text=c("semivariogram","covariance")[cov+1], side=2, outer = TRUE, line=2) mtext(text=c("semivariogram","covariance")[cov+1], side=2, outer = TRUE, line=2)
invisible(opar) invisible(opar)
...@@ -452,11 +454,11 @@ plot.gmCgram = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= N ...@@ -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 #' Check for anisotropy of a theoretical variogram
#' #'
#' Checks for anisotropy of a theoretical variogram or covariance function model #' Checks for anisotropy of a theoretical variogram or covariance function model
#' @param x variogram or covariance model object #' @param x variogram or covariance model object
#' @param tol tolerance for #' @param tol tolerance for
#' @param ... extra arguments for generic functionality #' @param ... extra arguments for generic functionality
#' #'
#' @return Generic function. Returns of boolean answering the question of the name, #' @return Generic function. Returns of boolean answering the question of the name,
...@@ -511,7 +513,7 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){ ...@@ -511,7 +513,7 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){
#' Variogram method for gmSpatialModel objects #' Variogram method for gmSpatialModel objects
#' #'
#' Compute the empirical variogram of the conditioning data contained in a [gmSpatialModel-class] object #' Compute the empirical variogram of the conditioning data contained in a [gmSpatialModel-class] object
#' #'
#' @param object a gmSpatialModel object containing spatial data. #' @param object a gmSpatialModel object containing spatial data.
...@@ -519,7 +521,7 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){ ...@@ -519,7 +521,7 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){
#' @param ... further parameters to [gstat::variogram()] #' @param ... further parameters to [gstat::variogram()]
#' #'
#' @return Currently the function is just a convenience wrapper on #' @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 #' and returns objects of class "\code{gstatVariogram}". Check the
#' help of \code{gstat::variogram} for further information. #' help of \code{gstat::variogram} for further information.
#' In the near future, methods will be created, which will depend on #' In the near future, methods will be created, which will depend on
...@@ -534,15 +536,15 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){ ...@@ -534,15 +536,15 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){
# Variogram calculations # Variogram calculations
# #
# Compute empricial variograms out of a spatial data object # Compute empricial variograms out of a spatial data object
# #
# @param object spatial data container # @param object spatial data container
# @param ... further parameters for variogram calculation # @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. # will be produced. See appropriate method descriptions.
# #
# @importFrom sp variogram # @importFrom sp variogram
# @export # @export
##variogram <- function(object, ...) UseMethod("variogram", object) ##variogram <- function(object, ...) UseMethod("variogram", object)
...@@ -558,25 +560,25 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){ ...@@ -558,25 +560,25 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){
#' Empirical variogram or covariance function in 2D #' 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 X matrix of Nx2 columns with the geographic coordinates
#' @param Z matrix or data.frame of data with dimension (N,Dv) #' @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 #' for covariance function, a (N,Dv)-matrix or a Dv-vector giving the mean values
#' @param maxdist maximum lag distance to consider #' @param maxdist maximum lag distance to consider
#' @param lagNr number of lags 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 #' @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 #' 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 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 #' 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 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 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 #' @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. #' represent either internal functions, or preliminary, not fully-tested functions. Use \code{\link{variogram}} instead.
#' @export #' @export
#' @family gmEVario functions #' @family gmEVario functions
...@@ -594,7 +596,7 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){ ...@@ -594,7 +596,7 @@ variogram_gmSpatialModel <- function(object, methodPars=NULL, ...){
#' vge["npairs",1] #' vge["npairs",1]
#' vge["lags",1] #' vge["lags",1]
gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), 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), 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], azimuthNr=4, azimuths = seq(from=0, to=180, length.out=azimuthNr+1)[1:azimuthNr],
maxbreadth=Inf, minpairs=10, cov=FALSE){ maxbreadth=Inf, minpairs=10, cov=FALSE){
...@@ -611,7 +613,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -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 # expand the information given into a set of columns stating conditions
if(length(dim(lags))==0){ 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 if(maxbreadth!=Inf) lags[,"maxbreadth"]=maxbreadth
}else if(dim(lags)==2){ }else if(dim(lags)==2){
lags = data.frame(lags) lags = data.frame(lags)
...@@ -621,21 +623,21 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -621,21 +623,21 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)),
if(length(dim(azimuths))==0){ if(length(dim(azimuths))==0){
tol = (azimuths[2]-azimuths[1])/2 tol = (azimuths[2]-azimuths[1])/2
if(is.na(tol)) tol=180 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){ }else if(dim(azimuths)==2){
azimuths = data.frame(azimuths) azimuths = data.frame(azimuths)
colnames(azimuths) = c("minaz","maxaz") colnames(azimuths) = c("minaz","maxaz")
}else stop("azimuths can be either a vector of lags or a data.frame, see ?gsi.EVario2D") }else stop("azimuths can be either a vector of lags or a data.frame, see ?gsi.EVario2D")
# compute pairs of locations # compute pairs of locations
ij = expand.grid(1:nrow(X), 1:nrow(X))# indices ij = expand.grid(1:nrow(X), 1:nrow(X))# indices
XX = X[ij[,1],]-X[ij[,2],] # locations 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 = gmApply(XX, 1, function(x) pi/2-atan2(x[2],x[1])) +2*pi
XXaz = XXaz %% pi # residual to 180?? XXaz = XXaz %% pi # residual to 180??
# compute residuals and pairs of variables appropriate to the structural function # compute residuals and pairs of variables appropriate to the structural function
if(cov){ if(cov){
if(all(dim(Z)==dim(Ff))){ if(all(dim(Z)==dim(Ff))){
Z = as.matrix(Z-Ff) Z = as.matrix(Z-Ff)
}else if(Dv==length(c(unlist(Ff)))){ }else if(Dv==length(c(unlist(Ff)))){
...@@ -648,8 +650,8 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -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 Z = lm(as.matrix(Z)~as.matrix(Ff)+0)$residuals ## ideally this should be a GLS fit
ZZ = Z[ij[,1],]-Z[ij[,2],] ZZ = Z[ij[,1],]-Z[ij[,2],]
} }
# output # output
## ATTENTION: needs to be changed to return a structure (3,Na)-matrix of objects, ## ATTENTION: needs to be changed to return a structure (3,Na)-matrix of objects,
# like logratioVariogramAnisotropy # like logratioVariogramAnisotropy
...@@ -667,7 +669,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -667,7 +669,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)),
if(ncol(lags)>2){ if(ncol(lags)>2){
tk_b = outer(xxabs * abs(sin((xxaz-(azs[i,2]-azs[i,1])))), lags[,3], "<=") tk_b = outer(xxabs * abs(sin((xxaz-(azs[i,2]-azs[i,1])))), lags[,3], "<=")
tk_h = tk_h & tk_b tk_h = tk_h & tk_b
} }
n[,i] = colSums(tk_h) n[,i] = colSums(tk_h)
for(j in 1:Nh){ for(j in 1:Nh){
if(!is.na(n[j,i]) && n[j,i]>minpairs){ if(!is.na(n[j,i]) && n[j,i]>minpairs){
...@@ -675,7 +677,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -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]) vg[j,,,i] = gmApply((ZZ[tk_a,,][tk_h[,j],,]), c(2,3),"sum")/(n[j,i])
}else{ # semi-variogram }else{ # semi-variogram
aux = rowSums( 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]) )/(2*n[j,i])
vg[j,,,i] = aux vg[j,,,i] = aux
} }
...@@ -685,7 +687,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -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])) return(list(gamma=vg[,,,i], lags=gsi.lagClass(lags), npairs =n[,i]))
}) })
# output # output
attr(res, "directions") = gsi.azimuthInterval(azimuths) attr(res, "directions") = gsi.azimuthInterval(azimuths)
# attr(res, "lags") = gsi.lagClass(lags) # attr(res, "lags") = gsi.lagClass(lags)
...@@ -695,34 +697,34 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -695,34 +697,34 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)),
} }
#' Empirical variogram or covariance function in 3D #' 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 X matrix of Nx3 columns with the geographic coordinates
#' @param Z matrix or data.frame of data with dimension (N,Dv) #' @param Z matrix or data.frame of data with dimension (N,Dv)
#' @param Ff for variogram, matrix of basis functions with nrow(Ff)=N #' @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!!); #' (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 #' for covariance function, a (N,Dv)-matrix or a Dv-vector giving the mean values
#' @param maxdist maximum lag distance to consider #' @param maxdist maximum lag distance to consider
#' @param lagNr number of lags 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 #' @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 #' 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 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 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 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 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 #' @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. #' represent either internal functions, or preliminary, not fully-tested functions. Use \code{\link{variogram}} instead.
#' #'
#' #'
#' @export #' @export
#' @family gmEVario functions #' @family gmEVario functions
#' #'
# @examples # @examples
# dt <- readr::read_csv("~/Geochem_sum_imp2.csv") # dt <- readr::read_csv("~/Geochem_sum_imp2.csv")
# #
# X = as.matrix(dt[,c("X","Y","Z")]) # X = as.matrix(dt[,c("X","Y","Z")])
# Z = as.matrix(compositions::clr(dt[,c("Cu","Zn","Pb", "As", "Cd", "In", "Other")])) # 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)) # 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)), ...@@ -736,11 +738,11 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)),
# vge["lags",1] # vge["lags",1]
# plot.gmEVario(vge, varnames = colnames(Z), commonAxis=FALSE, # 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) # 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, ## ie. the 4 planar directions on the upper triangle,
## the downhole direction in the lower triangle ## the downhole direction in the lower triangle
gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), 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), lagNr = 15, lags = seq(from=0, to=maxdist, length.out=lagNr+1),
dirvecs=t(c(1,0,0)), angtol=90, dirvecs=t(c(1,0,0)), angtol=90,
maxbreadth=Inf, minpairs=10, cov=FALSE){ maxbreadth=Inf, minpairs=10, cov=FALSE){
...@@ -757,15 +759,15 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -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 ## expand the information given into a set of columns stating conditions
# prepare lags # prepare lags
if(length(dim(lags))==0){ 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) # if(maxbreadth!=Inf)
lags[,"maxbreadth"] = maxbreadth lags[,"maxbreadth"] = maxbreadth
}else if(dim(lags)==2){ }else if(dim(lags)==2){
lags = data.frame(lags) lags = data.frame(lags)
colnames(lags) = c("minlag","maxlag","maxbreadth")[1:ncol(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") }else stop("lags can be either a vector of lags or a data.frame, see ?gsi.EVario3D")
lagsSq = lags^2 lagsSq = lags^2
# prepare directions with tolerance # prepare directions with tolerance
if(length(dim(dirvecs))==0){ if(length(dim(dirvecs))==0){
if(length(dirvecs)!=3) stop("dirvecs must be a matrix with 3 columns!") 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)), ...@@ -775,15 +777,15 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)),
if(length(angtol) %in% c(1, nrow(dirvecs))){ if(length(angtol) %in% c(1, nrow(dirvecs))){
dirvecs = cbind(dirvecs, tol=cos(angtol*pi/180) ) 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") }else stop("angtol can be either a scalar or a vector of length=nrow(dirvecs), see ?gsi.EVario3D")
# compute pairs of locations # compute pairs of locations
ij = expand.grid(1:nrow(X), 1:nrow(X))# indices ij = expand.grid(1:nrow(X), 1:nrow(X))# indices
XX = X[ij[,1],]-X[ij[,2],] # locations XX = X[ij[,1],]-X[ij[,2],] # locations
XXmodSq = gmApply(XX, 1, function(x) sum(x^2)) XXmodSq = gmApply(XX, 1, function(x) sum(x^2))
XXmod = sqrt(XXmodSq) XXmod = sqrt(XXmodSq)
# compute residuals and pairs of variables appropriate to the structural function # compute residuals and pairs of variables appropriate to the structural function
if(cov){ if(cov){
if(all(dim(Z)==dim(Ff))){ if(all(dim(Z)==dim(Ff))){
Z = as.matrix(Z-Ff) Z = as.matrix(Z-Ff)
}else if(Dv==length(c(unlist(Ff)))){ }else if(Dv==length(c(unlist(Ff)))){
...@@ -796,7 +798,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -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 Z = lm(as.matrix(Z)~as.matrix(Ff)+0)$residuals ## ideally this should be a GLS fit
ZZ = Z[ij[,1],]-Z[ij[,2],] ZZ = Z[ij[,1],]-Z[ij[,2],]
} }
# output # output
## ATTENTION: needs to be changed to return a structure (3,Na)-matrix of objects, ## ATTENTION: needs to be changed to return a structure (3,Na)-matrix of objects,
# like logratioVariogramAnisotropy # like logratioVariogramAnisotropy
...@@ -816,7 +818,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -816,7 +818,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)),
residualsSq = XXmodSq[tk_a]-projections[tk_a]^2 residualsSq = XXmodSq[tk_a]-projections[tk_a]^2
tk_b = outer(residualsSq, lagsSq[,3], "<=") tk_b = outer(residualsSq, lagsSq[,3], "<=")
tk_h = tk_h & tk_b tk_h = tk_h & tk_b
} }
n[,i] = colSums(tk_h) n[,i] = colSums(tk_h)
for(j in 1:Nh){ for(j in 1:Nh){
if(!is.na(n[j,i]) && n[j,i]>minpairs){ if(!is.na(n[j,i]) && n[j,i]>minpairs){
...@@ -824,7 +826,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -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]) vg[j,,,i] = gmApply((ZZ[tk_a,,][tk_h[,j],,]), c(2,3),"sum")/(n[j,i])
}else{ # semi-variogram }else{ # semi-variogram
aux = rowSums( 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]) )/(2*n[j,i])
vg[j,,,i] = aux vg[j,,,i] = aux
} }
...@@ -834,7 +836,7 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -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])) return(list(gamma=vg[,,,i], lags=gmGeostats:::gsi.lagClass(lags), npairs =n[,i]))
}) })
# output # output
attr(res, "directions") = gmGeostats:::gsi.directorVector(dirvecs[,1:3]) attr(res, "directions") = gmGeostats:::gsi.directorVector(dirvecs[,1:3])
# attr(res, "lags") = gsi.lagClass(lags) # attr(res, "lags") = gsi.lagClass(lags)
...@@ -848,13 +850,13 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -848,13 +850,13 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)),
#' Plot empirical variograms #' Plot empirical variograms
#' #'
#' Flexible plot of an empirical variogram of class gmEVario #' Flexible plot of an empirical variogram of class gmEVario
#' #'
#' @param x object to print, 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.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 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 #' on the upper triangle
#' @param vdir.lo ..., indices of the directions to be plotted on the lower triangle #' @param vdir.lo ..., indices of the directions to be plotted on the lower triangle
#' @param varnames variable names to be used #' @param varnames variable names to be used
...@@ -862,16 +864,16 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -862,16 +864,16 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)),
#' @param add boolean, add stuff to an existing diagram? #' @param add boolean, add stuff to an existing diagram?
#' @param commonAxis boolean, should vertical axes be shared by all plots in a row? #' @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 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 #' defaults to TRUE
#' @param ... further parameters to \code{\link{matplot}} #' @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`. #' 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 #' computed along certain distances, recorded in its attributes and retrievable with command
#' \code{\link{ndirections}}. #' \code{\link{ndirections}}.
#' @export #' @export
#' @family gmEVario functions #' @family gmEVario functions
#' @method plot gmEVario #' @method plot gmEVario
...@@ -884,8 +886,8 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)), ...@@ -884,8 +886,8 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)),
#' vge = gsi.EVario2D(X,Z) #' vge = gsi.EVario2D(X,Z)
#' plot(vge) #' plot(vge)
#' plot(vge, pch=22, lty=1, bg="grey") #' plot(vge, pch=22, lty=1, bg="grey")
plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= NULL, plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= NULL,
varnames = dimnames(x$gamma)[[2]], type="o", varnames = dimnames(x$gamma)[[2]], type="o",
add=FALSE, commonAxis=TRUE, cov =attr(x,"type")=="covariance", add=FALSE, commonAxis=TRUE, cov =attr(x,"type")=="covariance",
closeplot=TRUE, ...){ closeplot=TRUE, ...){
Dv = dim(x[1,1][[1]])[2] 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= ...@@ -912,12 +914,12 @@ plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo=
xlim.lo = c(0, maxdist) xlim.lo = c(0, maxdist)
} }
} }
opar = par() opar = par()
opar = par_remove_readonly(opar) opar = par_remove_readonly(opar)
if(closeplot) on.exit(par(opar)) if(closeplot) on.exit(par(opar))
myplot = function(...) matplot(type=type, ylab="", xlab="",xaxt="n", ...) myplot = function(...) matplot(type=type, ylab="", xlab="",xaxt="n", ...)
if(add) myplot = function(...) matpoints(type=type, ...) if(add) myplot = function(...) matpoints(type=type, ...)
if(!add){ if(!add){
...@@ -931,7 +933,7 @@ plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= ...@@ -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) 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) if(commonAxis) ylim=range(sapply(vdir.lo, function(kk) x["gamma",kk][[1]][,,j]), na.rm=TRUE)
myplot( 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, ...) sapply(vdir.lo, function(kk) x["gamma",kk][[1]][,i,j]), ylim=ylim, ...)
if(i==j){ if(i==j){
axis(side = 3) axis(side = 3)
...@@ -944,15 +946,15 @@ plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= ...@@ -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) 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) if(commonAxis) ylim=range(sapply(vdir.up, function(kk) x["gamma",kk][[1]][,,j]), na.rm=TRUE)
myplot( 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, ...) sapply(vdir.up, function(kk) x["gamma",kk][[1]][,i,j]), ylim=ylim, ...)
if(i==j){ if(i==j){
axis(side = 1) axis(side = 1)
mtext(text=varnames[i], side = 1, line=3) mtext(text=varnames[i], side = 1, line=3)
mtext(text=varnames[i], side = 2, line=3) mtext(text=varnames[i], side = 2, line=3)
} }
} }
}} }}
mtext(text="lag distance", side=1, outer = TRUE, line=0) mtext(text="lag distance", side=1, outer = TRUE, line=0)
mtext(text=c("semivariogram","covariance")[cov+1], side=2, outer = TRUE, line=2) mtext(text=c("semivariogram","covariance")[cov+1], side=2, outer = TRUE, line=2)
attr(opar, "vdir.up")=vdir.up 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= ...@@ -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 structural function to gmEVario format
#' #'
#' Convert empirical covariance functions or variograms to the format gmEVario #' Convert empirical covariance functions or variograms to the format gmEVario
#' of package gmGeostats #' of package gmGeostats
#' #'
#' @param vgemp variogram/covariance function to be converted #' @param vgemp variogram/covariance function to be converted
#' @param ... further parameters #' @param ... further parameters
#' #'
#' @return the empirical covariance function or variogram, recasted to class #' @return the empirical covariance function or variogram, recasted to class
#' \code{gmEVario}. This is a generic function. Methods exist for objects of #' \code{gmEVario}. This is a generic function. Methods exist for objects of
#' class \code{logratioVariogram}\code{logratioVariogramAnisotropy} #' class \code{logratioVariogram}\code{logratioVariogramAnisotropy}
#' (for compositional data) and \code{gstatVariogram} #' (for compositional data) and \code{gstatVariogram}
#' (from package \code{gstat}). #' (from package \code{gstat}).
#' @export #' @export
#' @aliases as.gmEVario.gstatVariogram as.gmEVario.logratioVariogram #' @aliases as.gmEVario.gstatVariogram as.gmEVario.logratioVariogram
#' as.gmEVario.logratioVariogramAnisotropy #' as.gmEVario.logratioVariogramAnisotropy
...@@ -995,26 +997,26 @@ variogramModelPlot <- function(vg, ...) UseMethod("variogramModelPlot", vg) ...@@ -995,26 +997,26 @@ variogramModelPlot <- function(vg, ...) UseMethod("variogramModelPlot", vg)
#' Quick plotting of empirical and theoretical variograms #' 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 vg empirical variogram or covariance function
#' @param model optional, theoretical variogram or covariance function #' @param model optional, theoretical variogram or covariance function
#' @param col colors to use for the several directional variograms #' @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 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 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 closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)?
#' defaults to TRUE #' defaults to TRUE
#' @param ... further parameters to underlying plot or matplot functions #' @param ... further parameters to underlying plot or matplot functions
#' #'
#' @return The function is primarily called for producing a plot. However, it #' @return The function is primarily called for producing a plot. However, it
#' invisibly returns the graphical parameters active before the call #' invisibly returns the graphical parameters active before the call
#' occurred. This is useful for constructing complex diagrams, by adding layers #' 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 #' 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 #' 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`. #' want to leave the plot open for further changes give the extra argument `closeplot=FALSE`.
#' @export #' @export
#' @family variogramModelPlot #' @family variogramModelPlot
#' @family gmEVario functions #' @family gmEVario functions
#' @family gmCgram functions #' @family gmCgram functions
#' @seealso [logratioVariogram()] #' @seealso [logratioVariogram()]
#' @method variogramModelPlot gmEVario #' @method variogramModelPlot gmEVario
#' #'
...@@ -1030,15 +1032,15 @@ variogramModelPlot <- function(vg, ...) UseMethod("variogramModelPlot", vg) ...@@ -1030,15 +1032,15 @@ variogramModelPlot <- function(vg, ...) UseMethod("variogramModelPlot", vg)
#' Z = as.matrix(jura.pred[,c("Zn","Cd","Pb")]) #' Z = as.matrix(jura.pred[,c("Zn","Cd","Pb")])
#' vge = gsi.EVario2D(X,Z) #' vge = gsi.EVario2D(X,Z)
#' variogramModelPlot(vge, vm) #' variogramModelPlot(vge, vm)
#' #'
#' #'
variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogramModelList object containing a variogram model fitted to vg variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogramModelList object containing a variogram model fitted to vg
col = rev(rainbow(ndirections(vg))), col = rev(rainbow(ndirections(vg))),
commonAxis = FALSE, commonAxis = FALSE,
newfig = TRUE, closeplot=TRUE, ...){ newfig = TRUE, closeplot=TRUE, ...){
opar = plot(vg, commonAxis=commonAxis, add=!newfig, col=col, closeplot=is.null(model), ...) opar = plot(vg, commonAxis=commonAxis, add=!newfig, col=col, closeplot=is.null(model), ...)
opar = par_remove_readonly(opar) opar = par_remove_readonly(opar)
if(closeplot) on.exit(par(opar)) if(closeplot) on.exit(par(opar))
vdir.lo = attr(opar, "vdir.lo") vdir.lo = attr(opar, "vdir.lo")
vdir.up = attr(opar, "vdir.up") vdir.up = attr(opar, "vdir.up")
...@@ -1046,8 +1048,8 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra ...@@ -1046,8 +1048,8 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra
return(invisible(opar)) return(invisible(opar))
# OTHERWISE: add the curves for the model # OTHERWISE: add the curves for the model
aux = as.directorVector(attr(vg, "directions")) aux = as.directorVector(attr(vg, "directions"))
if(!is.null(vdir.lo)) vdir.lo = aux[vdir.lo,] if(!is.null(vdir.lo)) vdir.lo = aux[vdir.lo,]
if(!is.null(vdir.up)) vdir.up = aux[vdir.up,] 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, ...) opar = plot(as.gmCgram(model), vdir.up= vdir.up, vdir.lo= vdir.lo, add=TRUE, cov =FALSE, ...)
#f = as.function(as.gmCgram(gg)) #f = as.function(as.gmCgram(gg))
#Dv = dim(vg$gamma)[2] #Dv = dim(vg$gamma)[2]
...@@ -1062,7 +1064,7 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra ...@@ -1062,7 +1064,7 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra
# hd = seq(from=xlim[1], to=xlim[2], length.out = 200) # hd = seq(from=xlim[1], to=xlim[2], length.out = 200)
# Y = outer(hd, 1:nrow(dirs.lo), function(h,k){ # Y = outer(hd, 1:nrow(dirs.lo), function(h,k){
# f(rep(0,Dg), dirs.lo[k,]*h)[i,j] # f(rep(0,Dg), dirs.lo[k,]*h)[i,j]
# }) # })
# matlines(hd, Y, col=col[vdir.lo], ...) # matlines(hd, Y, col=col[vdir.lo], ...)
# } # }
# if((i<=j)&(!is.null(vdir.up))){ # if((i<=j)&(!is.null(vdir.up))){
...@@ -1072,26 +1074,26 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra ...@@ -1072,26 +1074,26 @@ variogramModelPlot.gmEVario <- function(vg, model = NULL, # gstat or variogra
# hd = seq(from=xlim[1], to=xlim[2], length.out = 200) # hd = seq(from=xlim[1], to=xlim[2], length.out = 200)
# Y = outer(hd, 1:nrow(dirs.up), function(h,k){ # Y = outer(hd, 1:nrow(dirs.up), function(h,k){
# f(rep(0,Dg), dirs.up[k,]*h)[i,j] # f(rep(0,Dg), dirs.up[k,]*h)[i,j]
# }) # })
# matlines(hd, Y, col=col[vdir.up], ...) # matlines(hd, Y, col=col[vdir.up], ...)
# } # }
# }} # }}
invisible(opar) invisible(opar)
} }
#' Number of directions of an empirical variogram #' Number of directions of an empirical variogram
#' #'
#' Returns the number of directions at which an empirical variogram was computed #' Returns the number of directions at which an empirical variogram was computed
#' #'
#' @param x empirical variogram object #' @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 #' number of directions at which an empirical variogram was computed
#' @export #' @export
#' @family gmEVario functions #' @family gmEVario functions
#' @family gmCgram functions #' @family gmCgram functions
#' @seealso [logratioVariogram()], [gstat::variogram()] #' @seealso [logratioVariogram()], [gstat::variogram()]
ndirections <- function(x){ UseMethod("ndirections", x) } ndirections <- function(x){ UseMethod("ndirections", x) }
...@@ -1107,7 +1109,7 @@ ndirections.azimuth = function(x) length(x) ...@@ -1107,7 +1109,7 @@ ndirections.azimuth = function(x) length(x)
#' @method ndirections azimuthInterval #' @method ndirections azimuthInterval
#' @export #' @export
ndirections.azimuthInterval = function(x) length(x[[1]]) 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 #' @method ndirections logratioVariogramAnisotropy
#' @export #' @export
ndirections.logratioVariogramAnisotropy = function(x) ndirections(attr(x,"directions")) ndirections.logratioVariogramAnisotropy = function(x) ndirections(attr(x,"directions"))
...@@ -1124,7 +1126,7 @@ ndirections.gmEVario = function(x) ndirections(attr(x,"directions")) ...@@ -1124,7 +1126,7 @@ ndirections.gmEVario = function(x) ndirections(attr(x,"directions"))
#' @export #' @export
ndirections.gstatVariogram = function(x){ ndirections.gstatVariogram = function(x){
length(unique(paste(x$dir.hor, x$dir.ver))) length(unique(paste(x$dir.hor, x$dir.ver)))
} }
...@@ -1135,12 +1137,12 @@ ndirections.gstatVariogram = function(x){ ...@@ -1135,12 +1137,12 @@ ndirections.gstatVariogram = function(x){
gsi.azimuth = function(x){ gsi.azimuth = function(x){
class(x) = c("azimuth","directionClass") class(x) = c("azimuth","directionClass")
return(x) return(x)
} }
gsi.azimuthInterval = function(x){ gsi.azimuthInterval = function(x){
class(x) = c("azimuthInterval","directionClass") class(x) = c("azimuthInterval","directionClass")
return(x) return(x)
} }
...@@ -1155,14 +1157,14 @@ print.directionClass = function(x, complete=TRUE, ...){ ...@@ -1155,14 +1157,14 @@ print.directionClass = function(x, complete=TRUE, ...){
if(complete){ if(complete){
print(unclass(x), ...) print(unclass(x), ...)
} }
} }
as.directorVector <- function(x){ UseMethod("as.directorVector",x) } as.directorVector <- function(x){ UseMethod("as.directorVector",x) }
#' @method as.directorVector default #' @method as.directorVector default
as.directorVector.default = function(x,...) x as.directorVector.default = function(x,...) x
#' @method as.directorVector azimuth #' @method as.directorVector azimuth
...@@ -1185,7 +1187,7 @@ as.directorVector.azimuthInterval = function(x, D=2){ ...@@ -1185,7 +1187,7 @@ as.directorVector.azimuthInterval = function(x, D=2){
gsi.lagdists = function(x){ gsi.lagdists = function(x){
class(x) = c("lagdist","lagClass") class(x) = c("lagdist","lagClass")
return(x) return(x)
} }
gsi.lagClass = function(x){ gsi.lagClass = function(x){
if(length(dim(x))!=2) stop("provided lag classes are not matrix-like!") if(length(dim(x))!=2) stop("provided lag classes are not matrix-like!")
...@@ -1209,7 +1211,7 @@ gsi.midValues.default = function(x) x ...@@ -1209,7 +1211,7 @@ gsi.midValues.default = function(x) x
#' @method gsi.midValues lagClass #' @method gsi.midValues lagClass
gsi.midValues.lagClass <- function(x){ (x[[1]]+x[[2]])/2 } gsi.midValues.lagClass <- function(x){ (x[[1]]+x[[2]])/2 }
#' @method gsi.midValues azimuthIntervals #' @method gsi.midValues azimuthIntervals
gsi.midValues.azimuthInterval <- function(x){ (x[[1]]+x[[2]])/2 } gsi.midValues.azimuthInterval <- function(x){ (x[[1]]+x[[2]])/2 }
...@@ -1218,49 +1220,49 @@ 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 ## theoretical structural functions
# S3 -> S4 classes # S3 -> S4 classes
# cat("creating variogram model classes\n") # cat("creating variogram model classes\n")
# abstract classes # abstract classes
#' @title Structural function model specification #' @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 #' 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()]) #' class "variogramModel" or "variogramModelList" (see [gstat::vgm()])
#' @export #' @export
#' @include compositionsCompatibility.R #' @include compositionsCompatibility.R
#' @include gstatCompatibility.R #' @include gstatCompatibility.R
#' @include preparations.R #' @include preparations.R
setClassUnion(name="ModelStructuralFunctionSpecification", setClassUnion(name="ModelStructuralFunctionSpecification",
members=c("NULL","gmCgram", "LMCAnisCompo", "variogramModelList", "variogramModel")) members=c("NULL","gmCgram", "LMCAnisCompo", "variogramModelList", "variogramModel"))
# #
# #### container class -------------- # #### container class --------------
# # An S4 class to represent a Gaussian random field specification # # 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 # # (generalised) covariance function specification, typically an object
# # obtained from a call to functions such as \code{\link{setCgram}}, # # obtained from a call to functions such as \code{\link{setCgram}},
# # \code{\link{LMCAnisCompo}} or \code{gstat::vgm}. # # \code{\link{LMCAnisCompo}} or \code{gstat::vgm}.
# # @slot formula formula specifying the structure # # @slot formula formula specifying the structure
# # of dependence of the mean of the random field w.r.to spatial coordinates # # 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; # # 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 # # @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 # # 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 # # @return A object with the slots populated as given
# # @export # # @export
# # @seealso [gmSpatialModel-class], and the `make.gm*` functions referenced there # # @seealso [gmSpatialModel-class], and the `make.gm*` functions referenced there
# setClass("gmGaussianModel", # setClass("gmGaussianModel",
# slots = list(structure = "ModelStructuralFunctionSpecification", # slots = list(structure = "ModelStructuralFunctionSpecification",
# formula="formula", # formula="formula",
# beta = "structure") # beta = "structure")
# ) # )
# #
# #setMethod("initialize", signature="gmGaussianModel", # #setMethod("initialize", signature="gmGaussianModel",
# # def=function(.Object, structure, formula, beta){ # # def=function(.Object, structure, formula, beta){
# # .Object@formula = formula # # .Object@formula = formula
...@@ -1274,21 +1276,21 @@ setClassUnion(name="ModelStructuralFunctionSpecification", ...@@ -1274,21 +1276,21 @@ setClassUnion(name="ModelStructuralFunctionSpecification",
#### container class -------------- #### container class --------------
# An S4 class to represent a Gaussian random field specification # 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 # (generalised) covariance function specification, typically an object
# obtained from a call to functions such as \code{\link{setCgram}}, # obtained from a call to functions such as \code{\link{setCgram}},
# \code{\link{LMCAnisCompo}} or \code{gstat::vgm}. # \code{\link{LMCAnisCompo}} or \code{gstat::vgm}.
# @slot formula formula specifying the structure # @slot formula formula specifying the structure
# of dependence of the mean of the random field w.r.to spatial coordinates # 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; # 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 # @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 # 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 # @return A object with the slots populated as given
# @export # @export
# @seealso [gmSpatialModel-class], and the `make.gm*` functions referenced there # @seealso [gmSpatialModel-class], and the `make.gm*` functions referenced there
setClass("gmGaussianModel", setClass("gmGaussianModel",
slots = list(structure = "ModelStructuralFunctionSpecification", slots = list(structure = "ModelStructuralFunctionSpecification",
formula="formula", formula="formula",
beta = "numeric") beta = "numeric")
...@@ -1302,9 +1304,9 @@ setClass("gmGaussianModel", ...@@ -1302,9 +1304,9 @@ setClass("gmGaussianModel",
# abstract classes # abstract classes
#' @title Empirical structural function specification #' @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 #' (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()]) #' class "gstatVariogram" (see [gstat::variogram()])
#' @export #' @export
#' @include compositionsCompatibility.R #' @include compositionsCompatibility.R
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment