Newer
Older
#### 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
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#'
#' @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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
#'
#' @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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
#'
#' @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
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
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
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{
Raimon Tolosana-Delgado
committed
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")
}
Raimon Tolosana-Delgado
committed
## 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)
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
#' @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, ...))
}
}
#' @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 , ...)
#' @include gmSpatialMethodParameters.R
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"),
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
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
#' @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
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
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)
# }
#)