gsi.DS: Workhorse function for direct sampling

In gmGeostats: Geostatistics for Compositional Analysis

Workhorse function for direct sampling

Description

This function implements in R the direct sampling algorithm

Usage

gsi.DS(
  n,
  f,
  t,
  n_realiz,
  dim_TI,
  dim_SimGrid,
  TI_input,
  SimGrid_input,
  ivars_TI = 3:ncol(TI_input),
  SimGrid_mask = ncol(SimGrid_input),
  invertMask = TRUE
)

Arguments

n	size of the conditioning data event (integer)
f	fraction of the training image to scan (numeric between 0 and 1)
t	maximal acceptable discrepance between conditioning data event and TI event (numeric between 0 and 1)
n_realiz	number of simulations desired
dim_TI	dimensions of the grid of the training image (ie. either (n_x, n_y) for dimension k=2 or (n_x, n_y, n_z) for dimension k=3)
dim_SimGrid	dimensions of the simulation grid (ie. either (m_x, m_y) or (m_x, m_y, m_z))
TI_input	training image, as a matrix of (n_x\cdot n_y\cdot n_z, k+D) elements; WITH NAMED COLUMNS and including spatial coordinates
SimGrid_input	simulation grid with conditioning data, as a matrix of (m_x\cdot m_y\cdot m_z, k+D) elements; with same columns as TI_input
ivars_TI	which colnames of TI_input and SimGrid_input identify variables to consider in the data event
SimGrid_mask	either a logical vector of length m_x\cdot m_y\cdot m_z, or else a column name of SimGrid_input giving a logical column
invertMask	logical, does SimGrid_mask identify with TRUE the data OUTSIDE the simulation area?

Value

A sp::SpatialPixelsDataFrame() or sp::SpatialGridDataFrame(), depending on whether the whole grid is simulated. The '@data' slot of these objects contains a DataFrameStack() with the stacking dimension running through the realisations. It is safer to use this functionality through the interface make.gmCompositionalMPSSpatialModel(), then request a direct simulation with DSpars() and finally run it with predict_gmSpatialModel.

Author(s)

Hassan Talebi (copyright holder), Raimon Tolosana-Delgado

Examples

## training image:
x = 1:10
y = 1:7
xy_TI = expand.grid(x=x, y=y)
TI_input = cbind(xy_TI, t(apply(xy_TI, 1, function(x) c(sum(x), abs(x[2]-x[1]))+rnorm(2, sd=0.01))))
colnames(TI_input) = c("x", "y", "V1", "V2")
o1 = image_cokriged(TI_input, ivar="V1")
o2 = image_cokriged(TI_input, ivar="V2")
## simulation grid:
SimGrid = TI_input
SimGrid$mask = with(SimGrid, x==1 | x==10 | y==1 | y==7)
tk = SimGrid$mask
tk[sample(70, 50)] = TRUE 
SimGrid[tk,3:4]=NA
image_cokriged(SimGrid, ivar="V1", breaks=o1$breaks, col=o1$col)
image_cokriged(SimGrid, ivar="V2", breaks=o2$breaks, col=o2$col)
image_cokriged(SimGrid, ivar="mask", breaks=c(-0.0001, 0.5, 1.001))
## Not run: 
res = gsi.DS(n=5, f=0.75, t=0.05, n_realiz=2, dim_TI=c(10,7),  dim_SimGrid=c(10,7), 
       TI_input=as.matrix(TI_input), SimGrid_input=as.matrix(SimGrid), 
       ivars_TI = c("V1", "V2"), SimGrid_mask="mask", invertMask=TRUE)
image_cokriged(cbind(xy_TI, getStackElement(res,1)), ivar="V1", breaks=o1$breaks, col=o1$col)
image_cokriged(cbind(xy_TI, getStackElement(res,2)), ivar="V1", breaks=o1$breaks, col=o1$col)
image_cokriged(cbind(xy_TI, getStackElement(res,1)), ivar="V2", breaks=o2$breaks, col=o2$col)
image_cokriged(cbind(xy_TI, getStackElement(res,2)), ivar="V2", breaks=o2$breaks, col=o2$col)

## End(Not run)

gmGeostats documentation built on April 18, 2023, 5:08 p.m.