#### gstat easy/easier interface for multivariate data
setOldClass("gstat")
#' Fit an LMC to an empirical variogram
#'
#' Fit a linear model of coregionalisation to an empirical variogram
#'
#' @param v empirical variogram
#' @param ... further parameters
#' @export
#' @return Method fit_lmc.gstatVariogram is a wrapper around [gstat::fit.lmc()], that calls this function
#' and gives the resulting model its appropriate class (c("variogramModelList", "list")).
#' Method fit_lmc.default returns the fitted lmc (this function currently uses gstat as a
#' calculation machine, but this behavior can change in the future)
#' @aliases fit_lmc fit_lmc.default fit_lmc.logratioVariogramAnisotropy
#'
#' @examples
#' data("jura", package = "gstat")
#' X = jura.pred[,1:2]
#' Zc = jura.pred[,7:13]
#' gg = make.gmCompositionalGaussianSpatialModel(Zc, X, V="alr", formula = ~1)
#' vg = variogram(gg)
#' md = gstat::vgm(model="Sph", psill=1, nugget=1, range=1.5)
#' gg = fit_lmc(v=vg, g=gg, model=md)
#' variogramModelPlot(vg, model=gg)
fit_lmc <- function(v, ...) UseMethod("fit_lmc", v)
#' @describeIn fit_lmc wrapper around gstat::fit.lmc method
#' @param g spatial data object, containing the original data
#' @param model LMC or variogram model to fit
#' @export
#' @method fit_lmc gstatVariogram
fit_lmc.gstatVariogram <- function(v, g, model,...){
res = gstat::fit.lmc(as.gstatVariogram(v, ...), as.gstat(g, ...), as.variogramModel(model, ...))
class(res$model) = c("variogramModelList", "list")
return(res)
}
#' @describeIn fit_lmc flexible wrapper method for any class for which methods
#' for [as.gstatVariogram()], [as.gstat()] and [as.variogramModel()] exist.
#' In the future there may be direct specialised implementations not depending on
#' package gstat.
#' @export
#' @method fit_lmc default
fit_lmc.default <- function(v, g, model,...){
origclass = class(g)
res = fit_lmc(as.gstatVariogram(v), as.gstat(g), as.variogramModel(model), ...)$model
# activate in the future
res = as(res, origclass)
return(res)
}
#' @describeIn fit_lmc method for logratioVariogram wrapping compositions::fit.lmc.
#' In the future there may be direct specialised implementations,
#' including anisotropy (not yet possible).
#' @export
#' @method fit_lmc logratioVariogram
fit_lmc.logratioVariogram <- function(v, g, model,...){
res = compositions::fit.lmc(as.logratioVariogram(v), as.CompLinModCoReg(model), ...)
return(res)
}
#' Convert a regionalized data container to gstat
#'
#' Convert a regionalized data container to a "gstat" model object
#'
#' @param object regionalized data container
#' @param ... accessory parameters (currently not used)
#'
#' @return A regionalized data container of class "gstat",
#' eventually with variogram model included. See [gstat::gstat()] for more info.
#' @aliases as.gstat.default
#' @export
#'
#' @examples
#' data("jura", package = "gstat")
#' X = jura.pred[,1:2]
#' Zc = jura.pred[,7:13]
#' gg = make.gmCompositionalGaussianSpatialModel(Zc, X, V="alr", formula = ~1)
#' as.gstat(gg)
as.gstat <- function(object, ...) UseMethod("as.gstat", object)
#' @describeIn as.gstat default does nothing
#' @method as.gstat default
as.gstat.default <- function(object, ...){
return(object)
}
setGeneric("as.gstat", as.gstat)
# packs a regionalized composition and their geographic coordinates into a
# gstat object after an appropriate logratio representation
# coords: geographic coordinates (it works sure with data.frame)
# compo: an acomp object (NOT TRANSFORMED)
# V: can be either the matrix PSI (of the notes) or the strings "clr", "ilr" or "alr"
# nscore: should data be marginally transformed to normal scores?
# formulaterm: term for the formula argument of gstat (to control between UK and OK/SK)
# ...: further arguments to gstat (e.g. for controlling neighbourhood or specyfing a mean for SK)
compo2gstatLR = function(coords, compo, V=ilrBase(compo),
lrvgLMC=NULL, nscore=FALSE,
formulaterm = "~1", prefix=NULL, ...){
# prepare constants
V0 = V
D = ncol(compo)
o = gsi.produceV(V=V,D=D,orignames=colnames(compo),giveInv=FALSE, prefix=prefix)
prefix = o$prefix
V = o$V
# compute data (in lrs or in normal scores), set variable names
Zlr = compositions::idt(compo, V=V)
if(nscore){
source("nscore.R") # load the nscore.R functions
prefix= paste("NS",prefix,sep="")
Zlr = sapply(1:ncol(V), function(i){
rs = nscore(Zlr[,i])
aux = rs$nscore
attr(aux,"trn.table") = rs$trn.table # this ensures that the backtransformation is stored in the object
return(data.frame(aux))
})
Zlr = as.data.frame(Zlr)
}
if(is.null(colnames(Zlr))) colnames(Zlr) = paste(prefix, 1:(D-1), sep="")
# create gstat object
spatdescr = paste("~",c(paste(colnames(coords),collapse=" + ")), sep="")
gg = NULL
for(i in 1:(D-1)){
id = colnames(Zlr)[i]
frm = paste(id, formulaterm, sep="")
gg = gstat::gstat(gg, id=id, formula = stats::as.formula(frm), locations=stats::as.formula(spatdescr), data=data.frame(coords, Zlr), ...)
}
# if a logratio LMC was provided, convert it to gstat variogramModelList
# and attach it
if(!is.null(lrvgLMC) & !nscore){
# space for a future conversion of variation-variogram models to gstat-LR-variograms
gg$model = as.variogramModel(lrvgLMC, V=V0, prefix=prefix)
}
# return
return(gg)
}
## version for rmult
rmult2gstat = function(coords, data, V="cdt",
vgLMC=NULL, nscore=FALSE,
formulaterm = "~1", prefix=NULL, ...){
P = ncol(data)
if(nscore){
source("nscore.R") # load the nscore.R functions
prefix= paste("NS",prefix,sep="")
Z = sapply(1:P, function(i){
rs = nscore(data[,i])
aux = rs$nscore
attr(aux,"trn.table") = rs$trn.table # this ensures that the backtransformation is stored in the object
return(data.frame(aux))
})
Z = as.data.frame(Z)
}else{
Z = data
}
if(is.null(colnames(Z))) colnames(Z) = paste(prefix, 1:P, sep="")
# create gstat object
spatdescr = paste("~",c(paste(colnames(coords),collapse=" + ")), sep="")
gg = NULL
for(i in 1:P){
id = colnames(Z)[i]
frm = paste(id, formulaterm, sep="")
gg = gstat::gstat(gg, id=id, formula = stats::as.formula(frm), locations=stats::as.formula(spatdescr), data=data.frame(coords, Z), ...)
}
# if a logratio LMC was provided, convert it to gstat variogramModelList
# and attach it
if(!is.null(vgLMC) & !nscore){
# space for a future conversion of variation-variogram models to gstat-LR-variograms
gg$model = as.variogramModel(vgLMC, prefix=prefix)
}
# return
return(gg)
}
getGstatData = function(gg # gstat object
){
return(gg$data[[1]]$data@data)
}
#' Quick plotting of empirical and theoretical variograms
#' 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)?
#' 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 giving
#' argument `closeplot=FALSE` and then adding layers
#' of information. If you want to "freeze" your plot, either give `closeplot=TRUE` or
#' embed your call in another call to \code{\link{par}}, e.g. \code{par(variogramModelPlot(...))}.
#' @export
#' @importFrom gstat vgm
#' @seealso `gstat::plot.gstatVariogram()`
#' @method variogramModelPlot gstatVariogram
#' @family variogramModelPlot
#' @examples
#' data("jura", package="gstat")
#' X = jura.pred[,1:2]
#' Zc = jura.pred[,7:13]
#' gg = make.gmCompositionalGaussianSpatialModel(Zc, X, V="alr", formula = ~1)
#' vg = variogram(gg)
#' md = gstat::vgm(model="Sph", psill=1, nugget=1, range=1.5)
#' gg = fit_lmc(v=vg, g=gg, model=md)
#' variogramModelPlot(vg, model=gg)
variogramModelPlot.gstatVariogram =
function(vg, # gstatVariogram object
model = NULL, # gstat or variogramModelList object containing a variogram model fitted to vg
col = rev(rainbow(1+length(unique(vg$dir.hor)))),
commonAxis = FALSE,
newfig = TRUE,
closeplot = TRUE,
...
){
# capture graphical parameters
dotlist = list(...)
# check class of gg object and extract variable names
if(is.null(model)){
noms = levels(vg$id)
vrnames = sort(unique(unlist(strsplit(noms, split=".", fixed=TRUE))))
model = lapply(noms, function(i) gstat::vgm(model="Sph", psill=0, range=1, nugget=0))
names(model)=noms
}else{
if("variogramModelList" %in% class(model)){
vrnames = names(model)
vrnames = vrnames[-grep(".", vrnames, fixed=TRUE)]
}else if("gstat" %in% class(model)){
vrnames = names(model$data)
model = model$model
}else if("variogramModel" %in% class(model) ){
vrnames = levels(model$id)[1]
}else{
ggtr = tryCatch(as.variogramModel(model))
if(class(ggtr)=="try-error"){
stop("argument 'gg' must either be a variogramModel, a gstat object with a variogramModel, a variogramModelList object, or an object convertible to one")
}else{
return(variogramModelPlot(vg=vg, gg = ggtr, col = col, commonAxis = commonAxis, newfig = newfig, ...))
}
}
}
d = length(vrnames)
# plot empirical vario!
opar = par()
opar = par_remove_readonly(opar)
if(closeplot) on.exit(par(opar))
if(newfig) par(mfrow=c(d,d), mar=c(1,1,1,1)+0.5, oma=c(0,3,3,0))
for(i in 1:d){
for(j in 1:d){
if(i==j){
noms = vrnames[i]
}else{
noms = paste(vrnames[c(i,j)], vrnames[c(j,i)], sep=".")
}
# take the part of the empirical variogram corresponding to this (pair of) variable(s)
tk = vg$id %in% noms
empvar = vg[tk,]
# find relevant stats
azimuths = unique(empvar$dir.hor)
rgH = range(empvar$dist, na.rm=TRUE)
rgVG = range(empvar$gamma, na.rm=TRUE)
if(commonAxis){
tk.row = grep( vrnames[i],vg$id)
rgVG = range(vg[tk.row,"gamma"], na.rm=TRUE)
}
# take the corresponding model
modvar = model[which(names(model) %in% noms )][[1]]
# predict the several directions
myfun = function(azimuth){
dir = c(sin(azimuth*pi/180), cos((azimuth*pi/180)),0)
gstat::variogramLine(modvar, max(rgH), dir=dir)
}
preds = lapply(azimuths, myfun)
rgVG = range(0,rgVG, sapply(preds, function(X)X[,2]))
ldotlist = dotlist
if(!("xlim" %in% names(ldotlist))){
ldotlist$xlim = range(0,empvar$dist, na.rm=TRUE)
}
if(!("ylim" %in% names(ldotlist))){
ldotlist$ylim = range(ifelse(i==j,0,NA),rgVG, na.rm=TRUE)
}
if(!("pch" %in% names(ldotlist))){
ldotlist$pch = 19
}
if(!("lty" %in% names(ldotlist))){
ldotlist$lty = 2
}
if(!("lwd" %in% names(ldotlist))){
ldotlist$lwd = 2
}
if(!("xlab" %in% names(ldotlist))){
ldotlist$xlab = ""
}
if(!("ylab" %in% names(ldotlist))){
ldotlist$ylab = ""
}
ldotlist$col=col[as.integer(as.factor(empvar$dir.hor))]
ldotlist$x = empvar$dist
ldotlist$y = empvar$gamma
ldotlist$type="p"
do.call(plot, args=ldotlist)
sapply(1:length(azimuths), function(k){
intk = empvar$dir.hor==azimuths[k]
lines(gamma~dist, empvar[intk,], col=col[k], lty=ldotlist$lty)
lines(preds[[k]], col=col[k], lwd=ldotlist$lwd )
})
if(i==1) mtext(side=3, text = vrnames[j], line=2 )
if(j==1) mtext(side=2, text = vrnames[i], line=2 )
}
}
invisible(opar)
}
#### functions to change between LMC and empirical variograms from/to gstat -------
## as gmEVario
# @describeIn as.gmEVario
# @export
as.gmEVario.gstatVariogram = function(vgemp, ...) stop("not yet available")
## as.logratioVariogram (empirical) -------
# transforms a gstat empirical variogram into a logratioVariogram object (evtl. with anisotropy)
# @describeIn as.logratioVariogram
as.logratioVariogram.gstatVariogram = function(vgemp, # gstatVariogram object, emprical logratio variogram
V=NULL, # matrix or name of the logratio transformation used
tol=1e-12, # tolerance for generalized inverse (eventually for clr case)
orignames=NULL, # names of the original component
symmetrize=FALSE, # do you want a whole circle of directions?
...
){
# prepare dimensions, names and constants
DD = length(levels(vgemp$id))
D = (1+sqrt(1+8*DD))/2
if(is.null(orignames)) orignames = paste("v", 1:D, sep="")
if(length(orignames)!=D) stop("names provided not consistent with number of logratio variables. Did you forget the rest?")
o = gsi.produceV(V=V, D=D, orignames=orignames, giveInv=TRUE)
prefix = o$prefix
W = o$W
# separate each direction, if anisotropic vario (ATTENTION: 3D not yet supported)
vg4dir = split(vgemp, vgemp$dir.hor)
# function to build one logratioVariogram
buildOneLogratioVariogram = function(vg){
aux = split(vg, vg$id)
N = nrow(aux[[1]])
lrnames = unique(names(aux))
lrnames = sort(lrnames[-grep(".", lrnames, fixed=TRUE)])
h = array(aux[[1]][,2], dim=c(N, D, D), dimnames=list(NULL, orignames, orignames))
n = array(aux[[1]][,1], dim=c(N, D, D), dimnames=list(NULL, orignames, orignames))
v = array(0, dim=c(D-1, N, D-1), dimnames=list(lrnames, NULL, lrnames))
vvns = strsplit(names(aux), ".", fixed=TRUE)
for(i in 1:length(vvns)){
vns = vvns[[i]]
if(length(vns)==1){ # diagonal
v[vns, ,vns] = aux[[i]]$gamma
}else{#off-diagonal
v[vns[1], ,vns[2]] = aux[[i]]$gamma
v[vns[2], ,vns[1]] = aux[[i]]$gamma
}
}
d = D-1
dim(v) = c(N*d,d)
v = v %*% W
dim(v) = c(d, N*D)
v = t(W) %*% v
dim(v) = c(D,N,D)
# v = aperm(v, c(2,1,3))
v = gmApply(v,2,clrvar2variation)
dim(v) = c(D,D, N)
dimnames(v) = list(orignames, orignames, NULL)
v = aperm(v, c(3,1,2))
erg = structure(list(vg=v, h=h, n=n), class="logratioVariogram")
}
# create all variograms on all directions
res = sapply(vg4dir, buildOneLogratioVariogram)
# if a whole circle of directions is desired...
if(symmetrize){
cn = colnames(res)
res = cbind(res, res)
colnames(res) = c(cn, 180+as.double(cn))
}
# prepare and return the "logratioVariogramAnisotropy" object
hh = res["h",1][[1]][,1,1]
mndfh = mean(diff(hh))
dists = (hh[-1]+hh[-length(hh)])/2
attr(res, "dists") = c(0, dists, max(dists)+mndfh)
class(res)=c("logratioVariogramAnisotropy", "logratioVariogram")
return(res)
}
## as.gstatVariogram (empirical) -------
#' Represent an empirical variogram in "gstatVariogram" format
#'
#' Represent an empirical variogram in "gstatVariogram" format, from package "gstat"; see [gstat::variogram()]
#' for details.
#'
#' @param vgemp empirical variogram of any kind
#' @param ... further parameters (for generic functionality)
#'
#' @return The function returns an object of class "gstatVariogram" containing the empirical variogram provided.
#' See `gstat::variogram()` for details.
#' @export
#'
#' @examples
#' data("jura", package = "gstat")
#' X = jura.pred[,1:2]
#' Zc = compositions::acomp(jura.pred[,7:13])
#' lrvg = gmGeostats::logratioVariogram(data=Zc, loc=X)
#' as.gstatVariogram(lrvg, V="alr")
as.gstatVariogram <- function(vgemp, ...) UseMethod("as.gstatVariogram", vgemp)
#' @describeIn as.gstatVariogram Represent an empirical variogram in "gstatVariogram" format
#' @method as.gstatVariogram default
as.gstatVariogram.default <- function(vgemp, ...) vgemp
#' @describeIn as.gstatVariogram Represent an empirical variogram in "gstatVariogram" format
#' @method as.gstatVariogram gmEVario
as.gstatVariogram.gmEVario <- function(vgemp,...) stop("not yet available")
#' @describeIn as.gstatVariogram Represent an empirical variogram in "gstatVariogram" format
#' @method as.gstatVariogram logratioVariogram
#' @export
#' @param V eventually, indicator of which logratio should be used (one of: a matrix of logcontrasts, or of the strings "ilr", "alr" or "clr")
#' @param dir.hor eventually, which horizontal direction is captured by the variogram provided (seldom to be touched!)
#' @param dir.ver eventually, which vertical direction is captured by the variogram provided (seldom to be touched!)
#' @param prefix prefix name to use for the variables created (seldom needed)
as.gstatVariogram.logratioVariogram =
function(vgemp, # gstatVariogram object, emprical logratio variogram
V=NULL, # matrix or name of the logratio transformation used
dir.hor=0,
dir.ver=0,
prefix=NULL,
...){
class(vgemp) = NULL
orignames = dimnames(vgemp$vg)[[2]]
D = dim(vgemp$vg)[2]
o = gsi.produceV(V=V, D=D, orignames = orignames, giveInv = F, prefix=prefix )
V = o$V
prefix = o$prefix
newnames = paste(prefix, 1:(D-1), sep="")
Vu = V %*% diag(1/sqrt(colSums(V^2)))
hh = gmApply(vgemp$h, 1, clrvar2ilr, V=Vu^2)
nn = gmApply(vgemp$n, 1, clrvar2ilr, V=Vu^2)
vv = -0.5*apply(vgemp$vg, 1, clrvar2ilr, V=V)
ids = outer(newnames, newnames, paste, sep=".")
diag(ids) = newnames
ordre = NULL
for(i in nrow(ids):1){
ordre = c(ordre, ids[1:i,i])
}
rownames(nn) = ids
rownames(hh) = ids
rownames(vv) = ids
# data.frame: np, dist, gamma, dir.hor, dir.ver=0, id= factor
erg = data.frame(np=c(t(nn[ordre,])), dist=c(t(hh[ordre,])),
gamma=c(t(vv[ordre,])),
dir.hor=rep(dir.hor, each=ncol(nn)),
dir.ver=rep(dir.ver, each=ncol(nn)),
id=factor(rep(1:length(ordre), each=ncol(nn)), labels=ordre)
)
class(erg) = c("gstatVariogram", "data.frame")
return(erg)
}
#' @describeIn as.gstatVariogram Represent an empirical variogram in "gstatVariogram" format
#' @method as.gstatVariogram logratioVariogramAnisotropy
#' @export
as.gstatVariogram.logratioVariogramAnisotropy =
function(vgemp, # gstatVariogram object, emprical logratio variogram
V=NULL, # matrix or name of the logratio transformation used
...){
class(vgemp) = NULL
alphas = gsi.azimuth2angle(colnames(vgemp))
erg = as.gstatVariogram.logratioVariogram(vgemp[,1], V=V, dir.hor=alphas[1],...)
for(i in 2:length(alphas)){
erg = rbind(erg, as.gstatVariogram.logratioVariogram(vgemp[,i], V=V, dir.hor=alphas[i],...))
}
class(erg) = c("gstatVariogram", "data.frame")
return(erg)
}
## as.variogramModel (LMC) -------
#' Convert an LMC variogram model to gstat format
#'
#' Convert a linear model of coregionalisation to the format of package gstat. See [gstat::vgm()] for details.
#'
#' @param m variogram model
#' @param ... further arguments for generic functionality
#'
#' @return The LMC model specified in the format of package gstat, i.e. as the result
#' of using [gstat::vgm()]
#' @export
#' @importFrom gstat gstat
#'
#' @examples
#' data("jura", package = "gstat")
#' X = jura.pred[,1:2]
#' Zc = compositions::acomp(jura.pred[,7:13])
#' lrmd = compositions::CompLinModCoReg(formula=~nugget()+sph(1.5), comp=Zc)
#' as.variogramModel(lrmd, V="alr")
as.variogramModel <- function(m, ...) UseMethod("as.variogramModel", m)
#' @describeIn as.variogramModel Convert an LMC variogram model to gstat format
#' @method as.variogramModel default
#' @export
as.variogramModel.default <- function(m, ...) m
#' @describeIn as.variogramModel Convert an LMC variogram model to gstat format
#' @method as.variogramModel gmCgram
#' @export
as.variogramModel.gmCgram = function(m, ...) stop("not yet available")
#' @describeIn as.variogramModel Convert an LMC variogram model to gstat format
#' @method as.variogramModel LMCAnisCompo
#' @export
#' @param V eventually, specification of the logratio representation to use
#' for compositional data (one of: a matrix of log-contrasts to use, or else one of
#' the strings "alr", "clr" or "ilr")
#' @param prefix optional, name prefix for the generated variables if a transformation is used
#' @param ensurePSD logical, should positive-definiteness be enforced? defaults to TRUE, which may
#' produce several scary looking but mostly danger-free warnings
as.variogramModel.LMCAnisCompo <- function(m, V=NULL, prefix=NULL, ensurePSD=TRUE, ...){
D = ncol(m[,1]$sill)
d=D-1
o = gsi.produceV(V=V, D=D, orignames = dimnames(m["sill",1][[1]])[[2]], giveInv = FALSE, prefix=prefix)
V = o$V
prefix = o$prefix
# which combinations of variables do we have to consider?
noms = paste(prefix, 1:d, sep="")
combs = cbind(rep(1:d, times=d:1), matrix(1:d, ncol=d, nrow=d)[lower.tri(diag(d), diag=TRUE)] )
# equivalence table of correlogram model names
eqtabmodels = factor(c(nugget="Nug", exp="Exp", sph="Sph", gau="Gau"), levels=levels(vgm()[,1]))
models = eqtabmodels[sapply(1:ncol(m), function(j) m[,j]$model)]
# express all sill matrices in the desired logratio
# recode A in azimuth, range and range ratio
for(j in 1:ncol(m)){
aux = -0.5 * t(V) %*% m[,j]$sill %*% V
colnames(aux) <- rownames(aux) <- noms
if(ensurePSD){
e = eigen(aux)
tol = e$values[1]*1e-12
e$values[e$values<tol]=tol
aux = e$vectors %*% diag(e$values) %*% t(e$vectors)
warning("as.variogramModel.LMCAnisCompo: negative eigenvalues corrected")
}
m[,j]$sill = aux
}
anis = matrix(0, nrow=ncol(m), ncol=3)
colnames(anis) = c("range", "ang1", "anis1")
for(j in 1:ncol(m)){
anis[j, "anis1"] <- sqrt(sum((m[,j]$A[,2])^2))
anis[j, "range"] <- m[,j]$range/anis[j, "anis1"]
anis[j, "ang1"] <- atan2(-m[,j]$A[2,1], m[,j]$A[1,1]) * 180/pi
}
anis = data.frame(anis, ang2=0, anis2=1, ang3=0, kappa=0.5*(models!="Nug"))
# if(all(anis$anis1==1)) anis=NULL
# function to build one case
myfun = function(ij){
i = combs[ij,1]
j = combs[ij,2]
sills = sapply(1:ncol(m), function(k) m[,k]$sill[i,j] )
objecte = data.frame(model=models, psill=sills)
if(!is.null(anis)) objecte = cbind(objecte, anis)
# class(objecte) = c("variogramModel","data.frame" )
nugrow = which(objecte$model=="Nug")
nugget = ifelse(length(nugrow)>0, objecte[nugrow,"psill"],0)
first = (1:nrow(objecte))[-nugrow][1]
md = vgm(model=objecte[first,"model"], psill=objecte[first,"psill"], range=objecte[first,"range"],
nugget=nugget, anis=unlist(objecte[first,c("ang1","ang2","ang3", "anis1", "anis2")]),
kappa = objecte[first, "kappa"])
if(length((1:nrow(objecte))[-nugrow])>1)
for(kk in (1:nrow(objecte))[-nugrow][-1])
md = vgm(add.to=md, model=objecte[kk,"model"], psill=objecte[kk,"psill"], range=objecte[kk,"range"],
anis=unlist(objecte[kk,c("ang1","ang2","ang3", "anis1", "anis2")]),
kappa = objecte[kk, "kappa"])
return(md)
}
res = lapply(1:nrow(combs), myfun)
namelist = sapply(1:nrow(combs), function(ij) ifelse(combs[ij,1]==combs[ij,2], noms[combs[ij,1]], paste( noms[combs[ij,1]], noms[combs[ij,2]], sep=".") ) )
names(res) = namelist
#rownames(res) = NULL
class(res) = c("variogramModelList","list")
return(res)
}
#' @describeIn as.variogramModel Convert an LMC variogram model to gstat format
#' @method as.variogramModel CompLinModCoReg
#' @export
#' @param V eventually, specification of the logratio representation to use
#' for compositional data (one of: a matrix of log-contrasts to use, or else one of
#' the strings "alr", "clr" or "ilr")
#' @param prefix optional, name prefix for the generated variables if a transformation is used
#' @param ensurePSD logical, should positive-definiteness be enforced? defaults to TRUE, which may
#' produce several scary looking but mostly danger-free warnings
#' @importFrom compositions vgram.nugget vgram.cardsin vgram.exp vgram.gauss vgram.lin vgram.pow vgram.sph
as.variogramModel.CompLinModCoReg <- function(m, V="alr", prefix=NULL, ensurePSD=TRUE, ...){
strucs = gsi.extractCompLMCstructures(m)
D = ncol(strucs$sills[[1]])
o = gsi.produceV(V=V, D=D, giveInv = FALSE, prefix=prefix)
as.variogramModel(as.LMCAnisCompo(m, varnames=rownames(o$V)), V=o$V, prefix=o$prefix, ensurePSD=ensurePSD, ...)
}
## as.LMCAnisCompo (LMC) -------
#' @describeIn as.LMCAnisCompo Recast compositional variogram model to format LMCAnisCompo
#' @method as.LMCAnisCompo gstat
#' @export
as.LMCAnisCompo.gstat <- function(m,...) as.LMCAnisCompo(m$model, ...)
# @describeIn as.LMCAnisCompo Recast compositional variogram model to format LMCAnisCompo
gstat2LMCAnisCompo <- as.LMCAnisCompo.gstat
#' @describeIn as.LMCAnisCompo Recast compositional variogram model to format LMCAnisCompo
#' @method as.LMCAnisCompo variogramModelList
#' @export
as.LMCAnisCompo.variogramModelList <-
function(m, V=NULL, orignames=NULL, ...){
# prepare constants
DD = length(m)
D = (1+sqrt(1+8*DD))/2
if(is.null(orignames)) orignames = paste("v", 1:D, sep="")
if(length(orignames)!=D) stop("names provided not consistent with number of logratio variables. Did you forget the rest?")
o = gsi.produceV(V=V, D=D, orignames = orignames, giveInv = TRUE)
W = o$W
colnames(W) = orignames
# extract the dimensions and lr-names
lrnames = unique(names(m))
lrnames = sort(lrnames[-grep(".", lrnames, fixed=TRUE)]) # consider only names without "."
# extract the number of structures
K = nrow(m[[1]])
if(length(lrnames)!=(D-1))stop("dimensions of data implied from m and V do not fit!")
# equivalence table of correlogram model names
eqtabmodels = c(Nug="nugget", Exp="exp", Sph="sph", Gau="gau")
# function that extracts one particular sill matrix from the object
darstellung = function(m, i){
sill = matrix(0, ncol=D-1, nrow=D-1)
rownames(sill) <- colnames(sill) <- lrnames
for(jk in 1:DD){
vvns = strsplit(names(m), ".", fixed=TRUE)[[jk]]
if(length(vvns)==1) vvns = rep(vvns, 2)
sill[vvns[1], vvns[2]] <- sill[vvns[2], vvns[1]] <- m[[jk]]$psill[i]
}
return(sill)
}
setvariostructure = function(i){
model = eqtabmodels[ m[[1]]$model[i] ]
sill = clrvar2variation(t(W) %*% darstellung(m, i) %*% W)
range = m[[1]]$range[i] * m[[1]]$anis1[i]
A = anis2D.par2A(ratio=m[[1]]$anis1[i], angle=m[[1]]$ang1[i])
rs = list(model=model, range=range, A=A, sill=sill)
class(rs) = "variostructure"
return(rs)
}
res = sapply(1:K, setvariostructure)
class(res) = "LMCAnisCompo"
return(res)
}
#' @describeIn as.gmCgram Convert theoretical structural functions to gmCgram format
#' @method as.gmCgram variogramModelList
as.gmCgram.variogramModelList = function(m, ...) stop("not yet available")