-
Raimon Tolosana-Delgado authoredclass determination now with is(...) or inherits(...); calls to FORTRAN within C with string arguments platform independent
Raimon Tolosana-Delgado authoredclass determination now with is(...) or inherits(...); calls to FORTRAN within C with string arguments platform independent
gmSpatialModel.R 31.98 KiB
#### gmSpatialModel ------
#' Conditional spatial model data container
#'
#' This class is devised to contain a conditional spatial model, with: some conditioning data
#' (a [sp::SpatialPointsDataFrame()]), an unconditional geospatial model (a structure with e.g.
#' a training image; or the information defining a Gaussian random field); and eventually some
#' extra method parameters. The class extends [sp::SpatialPointsDataFrame()] and has therefore its slots,
#' plus `model` (for the unconditional model) and `parameters` (for the extra method information)
#'
#' @slot data a data.frame (or class extending it) containing the conditional data
#' @slot coords a matrix or dataframe of 2-3 columns containing the sampling locations of the conditional data
#' @slot coords.nrs see [sp::SpatialPointsDataFrame()]
#' @slot bbox see [sp::SpatialPointsDataFrame()]
#' @slot proj4string see [sp::SpatialPointsDataFrame()]
#' @slot model gmUnconditionalSpatialModel. Some unconditional geospatial model. It can be NULL.
#' @slot parameters gmSpatialMethodParameters. Some method parameters. It can be NULL
#' @family gmSpatialModel
#'
#' @return You will seldom create the spatial model directly. Use instead the creators `make.gm*` linked below
#' @export
#' @include abstractClasses.R
#' @importClassesFrom sp SpatialPointsDataFrame
#'
#' @examples
#' data("jura", package="gstat")
#' library(sp)
#' X = jura.pred[,1:2]
#' Zc = jura.pred[,7:13]
#' spdf = sp::SpatialPointsDataFrame(coords=X, data=Zc)
#' new("gmSpatialModel", spdf)
#' make.gmCompositionalGaussianSpatialModel(data=Zc, coords=X, V="alr")
setClass("gmSpatialModel", contains="SpatialPointsDataFrame",
slots = list(model="gmUnconditionalSpatialModel",
parameters="gmSpatialMethodParameters")
)
### creators --------
# make.gmMultilayerSpatialModel
# make.gmMultivariateSpatialModel
# make.gmAmountsSpatialModel
# make.gmUnivariateSpatialModel
# ...
#' Construct a Gaussian gmSpatialModel for regionalized multivariate data
#'
#' Construct a regionalized multivariate data container to be used for Gaussian-based geostatistics: variogram modelling, cokriging and simulation.
#' @param data either a data set of any data.frame similar class, or else a [sp::SpatialPointsDataFrame()] containing it
#' @param coords the coordinates of the sampling locations, if no SpatialPointsDataFrame was provided
#' @param model a variogram model, of any relevant class
#' @param beta (see `formula`) the coefficients of dependence of the mean of the random field, if these are known; e.g. if `formula=~1` constant mean,
#' and the mean is indeed known, `beta` would be a compositional mean; seldom used directly
#' @param formula a formula without left-hand-side term, e.g. ~1 or ~Easting+Northing, specifying what do we know of the
#' dependence of the mean of the random field; this follows the same ideas than in [gstat::gstat()]
#' @param ng optional neighborhood information, typically created with [KrigingNeighbourhood()]
#' @param nmax optional, neighborhood description: maximum number of data points per cokriging system
#' @param nmin optional, neighborhood description: minimum number of data points per cokriging system
#' @param omax optional, neighborhood description: maximum number of data points per cokriging system per quadrant/octant
#' @param maxdist optional, neighborhood description: maximum radius of the search neighborhood
#' @param force optional logical, neighborhood description: if not `nmin` points are found inside `maxdist` radius,
#' keep searching. This and all preceding arguments for neighborhood definition are borrowed from [gstat::gstat()]
#'
#' @return A "gmSpatialModel" object with all information provided appropriately structured. See [gmSpatialModel-class].
#' @export
#' @family gmSpatialModel
#' @seealso [SequentialSimulation()], [TurningBands()] or [CholeskyDecomposition()] for specifying the exact
#' simulation method and its parameters, [predict_gmSpatialModel] for running predictions or simulations
#'
#' @examples
#' data("jura", package="gstat")#' X = jura.pred[,1:2]
#' Zc = jura.pred[,7:13]
#' make.gmMultivariateGaussianSpatialModel(data=Zc, coords=X)
make.gmMultivariateGaussianSpatialModel <- function(
data, coords = attr(data, "coords"), ## data trench
model=NULL, beta=model$beta, formula=model$formula, ## model trench
ng=NULL, nmax=ng$nmax, nmin=ng$nmin, omax=ng$omax, maxdist=ng$maxdist, force=ng$force){ # method trench
# consider the class of 'data'
if(is(data,"Spatial")){
dt = data@data
}else{
dt = data
}
spdf = sp::SpatialPointsDataFrame(sp::SpatialPoints(coords), data=dt)
ng = KrigingNeighbourhood(nmax=c(nmax,Inf)[1], nmin=c(nmin,0)[1], omax=c(omax,0)[1], maxdist=c(maxdist,Inf)[1], force=c(force,FALSE)[1])
if(is.null(formula)) formula=~1
if(is.null(beta)) beta=Inf
if(is.null(model)){
vgm = new("gmGaussianModel", structure=NULL, beta=beta, formula=formula)
}else{
vgm = new("gmGaussianModel", structure=model, beta=beta, formula=formula)
}
res = new("gmSpatialModel", spdf, model=vgm, parameters=ng)
return(res)
}
#' Construct a Gaussian gmSpatialModel for regionalized compositions
#'
#' Construct a regionalized compositional data container to be used for Gaussian-based geostatistics: variogram modelling, cokriging and simulation.
#' @param data either a [compositions::acomp()] compositional data set, or else a [sp::SpatialPointsDataFrame()] containing it
#' @param coords the coordinates of the sampling locations, if no SpatialPointsDataFrame was provided
#' @param V optionally, a matrix of logcontrasts, or else one of the following strings: "alr", "ilr" or "clr";
#' to produce a plot of the empirical variogram in the corresponding representation; default to variation-variograms
#' @param prefix the desired prefix name for the logratio variables, if this is wished to be forced; otherwise derived from `V`
#' @param model a variogram model, of any relevant class
#' @param beta (see `formula`) the coefficients of dependence of the mean of the random field, if these are known; e.g. if `formula=~1` constant mean,
#' and the mean is indeed known, `beta` would be a compositional mean; seldom used directly
#' @param formula a formula without left-hand-side term, e.g. `~1` or `~Easting+Northing`, specifying what do we know of the
#' dependence of the mean of the random field; this follows the same ideas than in [gstat::gstat()]
#' @param ng optional neighborhood information, typically created with [KrigingNeighbourhood()]
#' @param nmax optional, neighborhood description: maximum number of data points per cokriging system
#' @param nmin optional, neighborhood description: minimum number of data points per cokriging system
#' @param omax optional, neighborhood description: maximum number of data points per cokriging system per quadrant/octant
#' @param maxdist optional, neighborhood description: maximum radius of the search neighborhood
#' @param force optional logical, neighborhood description: if not `nmin` points are found inside `maxdist` radius,
#' keep searching. This and all preceding arguments for neighborhood definition are borrowed from [gstat::gstat()]
#'
#' @return A "gmSpatialModel" object with all information provided appropriately structured. See [gmSpatialModel-class].
#' @export
#' @family gmSpatialModel
#' @seealso [SequentialSimulation()], [TurningBands()] or [CholeskyDecomposition()] for specifying the exact
#' simulation method and its parameters, [predict_gmSpatialModel] for running predictions or simulations
#'
#' @examples
#' data("jura", package="gstat")
#' X = jura.pred[1:20,1:2]
#' Zc = compositions::acomp(jura.pred[1:20,7:13])
#' make.gmCompositionalGaussianSpatialModel(data=Zc, coords=X, V="alr")
make.gmCompositionalGaussianSpatialModel <- function(
data, coords = attr(data, "coords"), V="ilr", prefix=NULL, ## data trench
model=NULL, beta=model$beta, formula=model$formula, ## model trench
ng=NULL, nmax=ng$nmax, nmin=ng$nmin, omax=ng$omax, maxdist=ng$maxdist, force=ng$force) # method trench
{
if(is(data,"Spatial")){
dt = data@data
}else{
dt = data
}
if(is.rmult(dt)){ W = gsi.getV(dt)
V = t(gsiInv(W))
warning("make.gmCompositionalSpatialModel: data of class 'rmult'! provided V will be ignored")
spdf = sp::SpatialPointsDataFrame(sp::SpatialPoints(coords), data=dt)
}else{
o = gsi.produceV(V=V,D=ncol(dt),orignames=colnames(dt),giveInv=FALSE, prefix=prefix)
prefix = o$prefix
V = o$V
if(is.null(rownames(V))) tryCatch(rownames(V) <- colnames(data))
if(is.null(rownames(V))) rownames(V) = paste("v", 1:nrow(V), sep="")
if(is.null(colnames(V))) colnames(V) = paste(prefix, 1:ncol(V), sep="")
spdf = sp::SpatialPointsDataFrame(sp::SpatialPoints(coords), data=ilr(dt,V))
}
if(is.null(ng)) ng = KrigingNeighbourhood(nmax=c(nmax,Inf)[1], nmin=c(nmin,0)[1], omax=c(omax,0)[1], maxdist=c(maxdist,Inf)[1], force=c(force,FALSE)[1])
if(is.null(formula)) formula=~1
if(is.null(beta)) beta=Inf
if(is.null(model)){
vgm = new("gmGaussianModel", structure=NULL, beta=beta, formula=formula)
}else{
vgm = new("gmGaussianModel", structure=model, beta=beta, formula=formula)
}
res = new("gmSpatialModel", spdf, model=vgm, parameters=ng)
return(res)
}
#' Construct a Multi-Point gmSpatialModel for regionalized compositions
#'
#' Construct a regionalized compositional data container to be used for multipoint geostatistics.
#' @param data either a [compositions::acomp()] compositional data set, or else a [sp::SpatialPointsDataFrame()] containing it
#' @param coords the coordinates of the sampling locations, if no SpatialPointsDataFrame was provided
#' @param V optionally, a matrix of logcontrasts, or else one of the following strings: "alr", "ilr" or "clr";
#' to produce a plot of the empirical variogram in the corresponding representation; default to variation-variograms
#' @param prefix the desired prefix name for the logratio variables, if this is wished to be forced; otherwise derived from `V`
#' @param model a training image, of any appropriate class (typically a [sp::SpatialGridDataFrame()] or [sp::SpatialPixelsDataFrame()])
#'
#' @return A "gmSpatialModel" object with all information provided appropriately structured. See [gmSpatialModel-class].
#' @export
#' @family gmSpatialModel
#' @seealso [DirectSamplingParameters()] for specifying a direct simulation method parameters,
#' [predict_gmSpatialModel] for running the simulation
make.gmCompositionalMPSSpatialModel = function(
data, coords = attr(data, "coords"), V="ilr", prefix=NULL, ## data trench
model=NULL) ## model trench
{
# extract data
if(is(data,"Spatial")){
dt.dt = data@data
}else{
dt.dt = data
}
# extract grid topology
if("grid" %in% slotNames(data)){
grid.dt = data@grid
}else if("grid" %in% slotNames(coords)){
grid.dt = coords@grid
}else{
grid.dt = NULL
warning("make.gmCompositionalMPSSpatialModel: grid topology is missing, it will be inferred from the data")
}
# check the adequacy of the model provided
if(is(model,"Spatial")){
dt.ti = model@data
grid.ti = if("grid" %in% slotNames(model)) model@grid else NULL
coords.ti = if(is(model, "SpatialGridDataFrame")) NULL else sp::coordinates(model)
}else if(is.data.frame(model) | is.array(model)){
warning("make.gmCompositionalMPSSpatialModel: model provided is not a Spatial object; conversion will be attempted") dt.ti = model
grid.ti = NULL
coords.ti = attr(model, "coords")
}else if(is.null(model)){
dt.ti = dt.dt
grid.ti = NULL
coords.ti = NULL
}else stop("make.gmCompositionalMPSSpatialModel: model provided cannot be interpreted")
# create grid for the TI or else check consistency
if(!is.null(grid.dt)){
if(is.null(grid.ti)){
x0 = grid.dt@grid@cellcentre.offset
Dx = grid.dt@grid@cellsize
n = dim(dt.ti)[1:2]
grid.ti = sp::GridTopology(cellcentre.offset = x0, cellsize = Dx, cells.dim = n)
}else{
stopifnot(all(grid.dt@cellsize==grid.ti@cellsize))
}
}else{
if(!is.null(grid.ti)){
x0 = grid.ti@grid@cellcentre.offset
Dx = grid.ti@grid@cellsize
n = grid.ti@grid@cellsize
grid.dt = sp::GridTopology(cellcentre.offset = x0, cellsize = Dx, cells.dim = n)
warning("make.gmCompositionalMPSSpatialModel: grid topology found on the TI but not on the data; TI from the grid will be used for the data!")
}
}
# check consistency between data parts of the data and model
if(ncol(dt.dt)!=ncol(dt.ti)) stop("make.gmCompositionalMPSSpatialModel: number of variables for the data and model provided do not coincide!")
# compute object data
if(is.rmult(dt.dt)){
W = gsi.getV(dt.dt)
V = t(gsiInv(W))
warning("make.gmCompositionalMPSSpatialModel: data of class 'rmult'! provided V will be ignored")
spdf = sp::SpatialPixelsDataFrame(sp::SpatialPoints(coords), data=dt.dt, grid=grid.dt)
}else{
o = gsi.produceV(V=V,D=ncol(dt.dt),orignames=colnames(dt.dt),giveInv=FALSE, prefix=prefix)
prefix = o$prefix
V = o$V
if(is.null(rownames(V))) tryCatch(rownames(V) <- colnames(dt.dt))
if(is.null(rownames(V))) rownames(V) = paste("v", 1:nrow(V), sep="")
if(is.null(colnames(V))) colnames(V) = paste(prefix, 1:ncol(V), sep="")
spdf = sp::SpatialPixelsDataFrame(sp::SpatialPoints(coords),
data=compositions::ilr(dt.dt,V), grid=grid.dt)
}
myilr = function(x,V){
res = log(as.matrix(x)) %*% V
tk = gmApply(is.na(x),1, any)
res[tk,] =NA
return(data.frame(res))
}
# compute object model
if(is(model, "SpatialGridDataFrame")){
model = sp::SpatialGridDataFrame(grid = sp::getGridTopology(model),
data = myilr(model@data,V))
}else if(is(model, "SpatialPixelsDataFrame")){
model = sp::SpatialPixelsDataFrame(points =
sp::SpatialPoints(sp::coordinates(model)),
grid = sp::getGridTopology(model),
data = myilr(model@data,V))
}else if(is.null(model)){
}else if(is.null(coords.ti)){
model = sp::SpatialGridDataFrame(grid = grid.ti, data = myilr(dt.ti,V))
}else{
model = sp::SpatialPixelsDataFrame(points = sp::SpatialPoints(coords.ti),
grid = grid.ti, data = myilr(dt.ti,V))
}
res = new("gmSpatialModel", spdf, model=model)
return(res)}
### exporter to gstat -----
as.gstat.gmSpatialModel <- function(object, ...){
# extra arguments
lldots = list(...)
# data elements
coords = sp::coordinates(object)
X = compositions::rmult(object@data, V= gsi.getV(object@data), orig=gsi.orig(object@data))
V = gsi.getV(X)
if(!is.null(V)){
# compositional case
compo = backtransform(X)
Vinv = t(gsiInv(V))
prefix = sub("1","",colnames(V)[1])
if(is.null(prefix) | length(prefix)==0) prefix = sub("1","",colnames(object@data)[1])
if(is.null(prefix) | length(prefix)==0) prefix = "ilr"
# model elements
if(!is(object@model, "gmGaussianModel")) stop("as.gstat: object@model must be of class 'gmGaussianModel'!")
if(is.null(aux <- object@model@structure)){
lrvgLMC = NULL
}else{
lrvgLMC = as.LMCAnisCompo(aux, V=Vinv)
}
formulaterm = paste(as.character(object@model@formula), collapse="")
beta = object@model@beta
if(any(is.infinite(beta))) beta = NULL
# neighbourhood
ng = object@parameters
# manage manual changes of parameters given in dots...
for(nm in names(ng)){
tk = grep(nm, names(lldots))
if(length(tk)>1) ng[[tk]] = lldots[[tk]]
}
# convert
if(!is(ng, "gmKrigingNeighbourhood")) stop("as.gstat: object@parameters must be of class 'gmGaussianMethodParameters'!")
res = compo2gstatLR(coords=coords, compo=compo, V=Vinv, lrvgLMC=lrvgLMC,
nscore=FALSE, formulaterm = formulaterm, prefix=prefix, beta=beta,
nmax=ng$nmax, nmin=ng$nmin, omax=ng$omax, maxdist=ng$maxdist, force=ng$force,
...)
return(res)
}else{
# non-compositional case
prefix = sub("1","",colnames(V)[1])
vgLMC <- object@model@structure
formulaterm = paste(as.character(object@model@formula), collapse="")
beta = object@model@beta
if(any(is.infinite(beta))) beta = NULL
# neighbourhood
ng = object@parameters
# manage manual changes of parameters given in dots...
for(nm in names(ng)){
tk = grep(nm, names(lldots))
if(length(tk)>1) ng[[tk]] = lldots[[tk]]
}
res = rmult2gstat(coords=coords, data=X, V=V, vgLMC=vgLMC,
nscore=FALSE, formulaterm = formulaterm, prefix=prefix, beta=beta,
nmax=ng$nmax, nmin=ng$nmin, omax=ng$omax, maxdist=ng$maxdist, force=ng$force,
...)
}
}
#' Recast spatial object to gmSpatialModel format
#'
#' Recast a spatial data object model to format gmSpatialModel
#'
#' @param object object to recast#' @param ... extra parameters for generic functionality
#'
#' @return The same spatial object re-structured as a "gmSpatialModel", see [gmSpatialModel-class]
#' @export
#' @family gmSpatialModel
as.gmSpatialModel <- function(object, ...) UseMethod("as.gmSpatialModel", object)
#' @describeIn as.gmSpatialModel Recast spatial object to gmSpatialModel format
#' @method as.gmSpatialModel default
as.gmSpatialModel.default = function(object, ...) object
#' @describeIn as.gmSpatialModel Recast spatial object to gmSpatialModel format
#' @param V optional, if the original data in the sptail object was compositional, which logcontrasts
#' were used to express it? Thsi can be either one string of "alr", "ilr" or "clr", or else a
#' (Dx(D-1))-element matrix of logcontrasts to pass from compositions to logratios
#' @method as.gmSpatialModel gstat
as.gmSpatialModel.gstat = function(object, V=NULL, ...){
stop("as.gmSpatialModel.gstat: not yet implemented")
}
## predict and simulate methods -------------
#' Predict method for objects of class 'gmSpatialModel'
#'
#' This is a one-entry function for several spatial prediction and simulation methods, for model objects
#' of class [gmSpatialModel-class]. The several methods are chosen by means of `pars` objects of the
#' appropriate class.
#'
#' @param object a complete "gmSpatialModel", containing conditioning data and unconditional model
#' @param newdata a collection of locations where a prediction/simulation is desired; this is typically
#' a [sp::SpatialPoints()], a data.frame or similar of X-Y(-Z) coordinates; or perhaps for gridded data
#' an object of class [sp::GridTopology()], [sp::SpatialGrid()] or [sp::SpatialPixels()]
#' @param pars parameters describing the method to use, *enclosed in an object of appropriate class*
#' (according to the method; see below)
#' @param ... further parameters for generic functionality, currently ignored
#'
#' @return Depending on the nature of `newdata`, the result will be a data container of the same kind,
#' extended with the predictions or simulations. For instance, if we want to obtain predictions on the
#' locations of a "SpatialPoints", the result will be a [sp::SpatialPointsDataFrame()]; if we want to obtain
#' simulations on the coordinates provided by a "data.frame", the result will be a [DataFrameStack()] with
#' the spatial coordinates stored as an extra attribute; or if the input for a simulation is a masked grid of class
#' [sp::SpatialPixels()], the result will be of class [sp::SpatialPixelsDataFrame()] which `data` slot will be
#' a [DataFrameStack].
#'
#' @details Package "gmGeostats" aims at providing a broad series of algorithms for geostatistical prediction
#' and simulation. All can be accesses through this interface, provided that arguments `object` and `pars` are of the
#' appropriate kind. In `object`, the most important criterion is the nature of its slot `model`. In `pars`
#' its class counts: for the creation of informative parameters in the appropriate format and class, a series
#' of accessory functions are provided as well.
#'
#' Classical (gaussian-based two-point) geostatistics are obtained if `object@model` contains a covariance function,
#' or a variogram model. Argument `pars` can be created with functions such as [KrigingNeighbourhood()],
#' [SequentialSimulation()], [TurningBands()] or [CholeskyDecomposition()] to respectively trigger a cokriging, as
#' sequential Gaussian simulation, a turning bands simulation, or a simulation via Cholesky decomposition.
#' The kriging neighbourhood can as well be incorporated in the "gmSpatialModel" `object` directly, or even be
#' nested in a "SequentialSimulation" parameter object.
#'
#' Conversely, to run a multipoint geostatistics algorithm, the first condition is that `object@model` contains a
#' training image. Additionally, `pars` must describe the characteristics of the algorithm to use. Currently, only
#' direct sampling is available: it can be obtained by providing some parameter object created with a call to
#' [DirectSamplingParameters()]. This method requires `newdata` to be on a gridded set of locations (acceptable
#' object classes are `sp::gridTopology`, `sp::SpatialGrid`, `sp::SpatialPixels`, `sp::SpatialPoints` or `data.frame`,
#' for the last two a forced conversion to a grid will be attempted).
#' @family gmSpatialModel
#' @name predict_gmSpatialModel
NULL
#' @rdname predict_gmSpatialModel
#' @export
#' @method predict gmSpatialModel
predict.gmSpatialModel <- function(object, newdata, pars=NULL, ...){
if(is.null(pars)){
return(Predict(object, newdata, ...))
}else{
return(Predict(object, newdata, pars, ...))
}
}
#' @rdname predict_gmSpatialModel
#' @export
setMethod("predict", signature(object="gmSpatialModel"), definition = predict.gmSpatialModel)
#' @rdname predict_gmSpatialModel
#' @include gmAnisotropy.R
#' @include preparations.R
#' @export
setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY"),
function(object, newdata, ...){
if(is.null(object@parameters)) object@parameters = KrigingNeighbourhood()
Predict(object, newdata, pars = object@parameters , ...)
}
)
#' @rdname predict_gmSpatialModel
#' @include gmSpatialMethodParameters.R
#' @export
setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmNeighbourhoodSpecification"),
function(object, newdata, pars, ...){
cat("starting cokriging \n")
object@parameters = pars
out = predict(as.gstat(object), newdata=newdata, ...)
return(out)
}
)
#' @rdname predict_gmSpatialModel
#' @export
setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmTurningBands"),
function(object, newdata, pars, ...){
stop("Turning Bands method not yet interfaced here; use")
cat("starting turning bands \n")
}
)
#' @rdname predict_gmSpatialModel
#' @export
setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmCholeskyDecomposition"),
function(object, newdata, pars, ...){
stop("Choleski decomposition method not yet implemented")
cat("starting Choleski decomposition \n")
}
)
#' @rdname predict_gmSpatialModel
#' @export
setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmSequentialSimulation"),
function(object, newdata, pars, ...){ cat("starting SGs \n")
object@parameters = pars$ng
erg = predict(as.gstat(object), newdata=newdata, nsim=pars$nsim, debug.level=pars$debug.level, ...)
Dg = ncol(object@coords)
erg = DataFrameStack(erg[,-(1:Dg)],
dimnames=list(
loc=1:nrow(erg), sim=1:pars$nsim,
var=colnames(object@data)
),
stackDim="sim")
attr(erg, "coords") = newdata
return(erg)
}
)
#' @rdname predict_gmSpatialModel
#' @export
setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmDirectSamplingParameters"),
function(object, newdata, pars, ...){
cat("starting direct sampling \n")
# extract training image
gt.ti = sp::getGridTopology(object@model)
dt.ti = as(object@model,"SpatialGridDataFrame")@data
# extract newdata mask
mask = getMask(newdata)
# fuse newdata, conditioning data and mask
if(is.null(newdata)){
gt.nw = NULL # object@grid
newdata = sp::SpatialPixelsDataFrame(points = as(object, "SpatialPoints"), data=object@data)
# tolerance = sqrt(sum(gt.nw@cellsize^2))/2 , grid=gt.nw)
}else if(is(newdata, "GridTopology")){
gt.nw = newdata
newdata = sp::SpatialPixelsDataFrame(points = as(object, "SpatialPoints"), data=object@data,
tolerance = sqrt(sum(gt.nw@cellsize^2))/2 , grid=gt.nw)
}else if(is(newdata, "SpatialGrid")){
gt.nw = sp::getGridTopology(newdata)
newdata = sp::SpatialPixelsDataFrame(points = as(object, "SpatialPoints"), data=object@data,
tolerance = sqrt(sum(gt.nw@cellsize^2))/2 , grid=gt.nw)
}else if(is(newdata, "SpatialPixels")){
gt.nw = sp::getGridTopology(newdata)
cc = rbind(sp::coordinates(newdata), sp::coordinates(object))
### WARNING: This is how it works!! COrrect from SpatialPoints and data.frame
aux = matrix(NA, nrow=nrow(sp::coordinates(newdata)), ncol=ncol(object@data))
colnames(aux) = colnames(object@data)
dt = rbind( aux , object@data)
newdata = sp::SpatialPixelsDataFrame(points = sp::SpatialPoints(cc), data=dt,
tolerance = sqrt(sum(gt.nw@cellsize^2))/2 , grid=gt.nw)
}else if(is(newdata, "SpatialPoints")){
gt.nw = object@model@grid
cc = rbind(sp::coordinates(object), sp::coordinates(newdata))
### WARNING: NOT YET TESTED WHAT HAPPENS WITH DUPLICATE LOCATIONS!!!
aux = matrix(NA, nrow=nrow(sp::coordinates(newdata)), ncol=ncol(object@data))
colnames(aux) = colnames(object@data)
dt = rbind( object@data, aux )
newdata = sp::SpatialPixelsDataFrame(points = sp::SpatialPoints(cc), data=dt,
tolerance = sqrt(sum(gt.nw@cellsize^2))/2) # , grid=gt.nw)
}else if(is.data.frame(newdata)){
gt.nw = NULL # object@grid
cc = rbind(sp::coordinates(object), newdata[,1:2])
if( all( colnames(object@data) %in% colnames(newdata) ) ){
aux = as.matrix(newdata[,colnames(object@data)])
}else{
aux = matrix(NA, nrow=nrow(sp::coordinates(newdata)), ncol=ncol(object@data))
}
colnames(aux) = colnames(object@data)
### WARNING: NOT YET TESTED WHAT HAPPENS WITH DUPLICATE LOCATIONS!!!
dt = rbind( object@data, aux ) newdata = sp::SpatialPixelsDataFrame(points = sp::SpatialPoints(cc), data=dt)
# tolerance = sqrt(sum(gt.nw@cellsize^2))/2 , grid=gt.nw)
}
if(is.null(gt.nw)) gt.nw = sp::getGridTopology(newdata)
dt.nw = as(newdata,"SpatialGridDataFrame")@data
dt.ti = as(object@model,"SpatialGridDataFrame")@data
if(all(gt.nw@cellsize!=gt.ti@cellsize))
stop("predict for gmSpatialModel with gmDirectSamplingParameters: inferred grid topologies for newdata, conditioning data and model do not coincide")
if(is.null(mask)) mask = rep(TRUE, nrow(dt.nw))
#erg = gsi.DS4CoDa(n=pars$patternSize, f=pars$scanFraction, t=pars$gof, n_realiz=pars$nsim,
# nx_TI=gt.ti@cells.dim[1], ny_TI=gt.ti@cells.dim[2],
# nx_SimGrid= gt.nw@cells.dim[1], ny_SimGrid=gt.nw@cells.dim[2],
# TI_input=as.matrix(dt.ti),
# SimGrid_input=as.matrix(dt.nw),
# V = "I", W=gsi.getV(object@data),
# ivars_TI = colnames(dt.ti),
# SimGrid_mask = mask,
# invertMask = FALSE)
erg = gsi.DS(n=pars$patternSize, f=pars$scanFraction, t=pars$gof, n_realiz=pars$nsim,
dim_TI=gt.ti@cells.dim, dim_SimGrid=gt.nw@cells.dim,
TI_input=as.matrix(dt.ti), SimGrid_input=as.matrix(dt.nw),
ivars_TI = colnames(dt.ti),
SimGrid_mask = mask,
invertMask = FALSE)
erg = gmApply(erg, FUN=ilrInv, V=gsi.getV(object@data), orig=gsi.orig(object@data))
erg = sp::SpatialGridDataFrame(grid = gt.nw, data=erg)
return(erg)
}
)
#' @describeIn gmSpatialModel Compute a variogram, see [variogram_gmSpatialModel()] for details
#' @inheritParams variogram_gmSpatialModel
#' @include variograms.R
#' @export
setMethod("variogram", signature=c(object="gmSpatialModel"),
def=variogram_gmSpatialModel)
#' @describeIn gmSpatialModel S4 wrapper method around [logratioVariogram()] for `gmSpatialModel`
#' objects
#' @inheritParams logratioVariogram
#' @inheritParams logratioVariogram_gmSpatialModel
#' @include compositionsCompatibility.R
#' @export
setMethod("logratioVariogram", "gmSpatialModel", logratioVariogram_gmSpatialModel)
#' @describeIn gmSpatialModel convert from gmSpatialModel to gstat; see [as.gstat()]
#' for details
#' @inheritParams as.gstat
#' @export
#' @include gstatCompatibility.R
setMethod("as.gstat", signature="gmSpatialModel", def=as.gstat.gmSpatialModel)
### tests -----------
# if(!exists("do.test")) do.test=FALSE
# if(do.test){
# library("gstat")
# library("magrittr")
# data("jura", package = "gstat")
# dt = jura.pred %>% dplyr::select(Cd:Zn)
# X = jura.pred[,1:2]
# a1 = make.gmCompositionalGaussianSpatialModel(acomp(dt), X)
# spc = SpatialPointsDataFrame(SpatialPoints(X),acomp(dt))
# a2 = make.gmCompositionalGaussianSpatialModel(spc)
# a3 = make.gmCompositionalGaussianSpatialModel(spc, V="alr")# a3gs = as.gstat(a3)
# a3vg = variogram(a3gs)
# plot(a3vg)
# a3gs = fit_lmc(a3vg, a3gs, vgm("Sph", psill=1, nugget=0, range=1))
# plot(a3vg, a3gs$model)
# a4 = make.gmCompositionalGaussianSpatialModel(spc, V="alr", model=a3gs$model)
# a4gs = as.gstat(a4)
# plot(a3vg, a4gs$model)
# }
### converters -----
# gmSpatialModel2SpatialGridDataFrame = function(from, to){
# tol = ifelse(is(grid, "GridTopology"), 0.5*sqrt(sum(from@grid@cellsize^2)), sqrt(.Machine$double.eps))
# to = SpatialPixelsDataFrame(
# grid = from@grid, data = from@data, points = SpatialPoints(from@coords),
# proj4string = proj4string(from), tolerance = tol
# )
# return(as(to, "SpatialGridDataFrame"))
# }
# SpatialGridDataFrame4gmSpatialModel = function(from, value){
# gt.mdl = from@model@grid
# gt.new = value@grid
# if(!is.null(gt.mdl))
# if(!all(gt.new@cellsize==gt.mdl@cellsize)) stop("replace(SpatialGridDataFrame->gmSpatialModel): provided sgdf grid topology inconsistent with existing model topology")
# from@grid = gt.new
# from@data = value@data
# from@coords = sp::coordinates(value)
# from@bbox = bbox(value)
# from@proj4string = proj4string(value)
# return(from)
# }
#setIs(class1 = "gmSpatialModel", class2 ="SpatialGridDataFrame",
# coerce = gmSpatialModel2SpatialGridDataFrame,
# replace = SpatialGridDataFrame4gmSpatialModel)
#setIs(class1 = "gmSpatialModel", class2 ="SpatialPointsDataFrame",
# coerce = function(from, to){
# to = SpatialPointsDataFrame(
# coords = SpatialPoints(from@coords), data = from@data,
# proj4string = proj4string(from), bbox = bbox(from)
# )
# return(to)
# }, replace = function(from, value){
# from@grid = NULL
# from@data = value@data
# from@coords = sp::coordinates(value)
# from@bbox = bbox(value)
# from@proj4string = proj4string(value)
# return(from)
# }
#)