gsi.DS {gmGeostats}

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)