simpgs {RGeostats}R Documentation

Plurigaussian Conditional Simulations


Plurigaussian conditional simulations


simpgs(dbin = NA, dbout, dbprop = NA,
       rule = rule.input(), model1 = NA, model2 = NA, 
       neigh=neigh.input(), props = NA, nostat=NA, domain=0,
       flag.gaus = FALSE, flag.prop=FALSE, flag.check=FALSE, = FALSE, seed = 232131, nbsimu = 1, nbtuba = 100,
       ngburn = 10, ngiter = 100, ngint=5, percent=0.1, toleps=1, 
       delta = NA, flag.spde=FALSE, triswitch="nqQ", gext=NA, numparts=NA,
       npiter=10, verbose=FALSE, = NA, accept.nloop.max = 10,
       accept.verbose=TRUE, radix = "SimPGS", = db.locmod())



The db-class structure containing the data file. If not provided, non-conditional simulations are performed


The db-class structure which will contain the simulation outputs


The db-class where the proportion variables must be located in the non-stationary case. See details in rule.check.


The rule-class structure which defines the Lithotype rule


The model-class structure containing the Model information for the first underlying gaussian random function. This model must be provided if the first gaussian random function is involved in the Lithotype rule.


The model-class structure containing the Model information for the second underlying gaussian random function. This model must be provided if the second gaussian random function is involved in the Lithotype rule


The neigh-class structure containing the Neighborhood information. Only the Unique or Bench neighborhoods are authorized.


Array giving the (constant) proportions of the different facies involved in the Lithotype rule. See details in rule.check.


List of non-stationary parameters. For details see model.param.define.


When (strictly) positive, it corresponds to the reference domain value. If a Db attribute is defined as a Domain variable in a Db file, a sample will be selected only if the corresponding Domain variable is equal to the reference domain value. This Domain selection can serve as a secondary selection for the sample in the Input Db or as a mask in the Output Db.


When FALSE, the simulation outcome is coded into facies. When TRUE, the (conditional) simulation is kept in the gaussian scale


When TRUE, the procedure returns the proportion of facies, rather than the simulations


When TRUE (and if dbin is defined), the facies of conditioning data is compared to the facies of the closest grid node.

When TRUE (and if dbin is defined), the grid node which coincides with a data is represented with the data facies (with a negative sign).


Seed used for the generation of random numbers. When 0, the seed is not initialized.


Number of simulations


Number of turning bands


Number of iterations performed for bootstrap


Maximum number of iterations for calculating the conditional expectation


Number of iterations inside the Gibbs sampler iterative algorithm (only for SPDE technology)


Amount of nugget effect added to the model if this initial model only contains Gaussian structures. If set to zero, the initial model is kept unchanged. This amount is defined as a percentage of the global variance.


Relative difference between the previous and the new value at a sample location below which the sample is considered as immobile


Spatial distance for generating the replicates in case of conditional simulations for shadow option.


When TRUE, the simulation is carried out using SPDE technique. Otherwise the Turning Bands method is used.


Command line for the internal triangulation step. For more information see meshing.


When 'dbout' is organized as a grid, it may be dilated by gext. This argument designates an array, with its dimension equal to the dimension of the space and which contains the extension defined in number of grid nodes. If not defined, the grid is not dilated and the simulated results may suffer some edge effect problems.


Array for subdividing the field into parts. Its dimension must be equal to the space dimension. If not defined (or equal to 1 in each space direction), the space is not subdivided. Subdividing in parts reduces the dimension of the matrices and can be used in the case of large files (input or output). In order to reduce the artifacts that may be induced by the subdivision, two steps of subdivisions are actually processed. For illustration purpose, let us assume that the field (SX by SY) with origin (X0,Y0) is subdivided into NX ny NY parts.

  • Step 1: The field is subdivided into NX by NY parts, with dimensions equal to SX/NX and SY/NY

  • Step 2: The field (with origin X0-SX/2, Y0-SY/2) is subdivided into NX+1 by NY+1 parts, with dimensions equal to SX/NX and SY/NY.


When the field is subdivided into several parts, several iterations are necessary to glue the parts. This parameter defined the number of iterations to be processed.


Verbose flag

An acceptation function. See details for more information.


Maximum number of iterations before 'nbsimu' acceptable simulation outcomes is reached.


When TRUE, the activity of the acceptation function is echoed.


Radix of the name given to the simulation outcomes stored in the output Db.

Decides whether or not the newly created variables will have their locator defined or not. For more information, see db.locmod.


The data Db where the simulations outcomes have been added.


# Acceptation function #
accept <- function(db,iatt) 
                                        # Locations of wells in the grid #
  ranks = c(5021, 5081)

                                        # Perform the Connected components #
  a = db
  a = db.locate(a,iatt,"z")
  a = morpho(a,0.5,1.5,oper="cc")

                                        # Do all wells belong to same CC #
  test = TRUE
  if (length(unique(a[ranks,a$natt])) > 1) test = FALSE

# Perform a conditional simulation of facies #
nbsimu = 25
grid   = db.create(nx=c(100,100))
model  = model.create(vartype="Cubic",range=30)
rule   = rule.create(c("S","F1","F2"))
x1 = c(20,80,50)
x2 = c(50,50,50)
z1 = c(1,1,2)
neigh = neigh.create(type=0)

grid = simpgs(data,grid,rule=rule,model1=model,neigh=neigh,,nbsimu=nbsimu,nbtuba=1000)
grid[,"SimPGS*"] = ifelse(grid[,"SimPGS*"] == 1,1,0)
grid =,names="SimPGS*",fun="mean")
plot(grid,title="Connectivity probability",pos.legend=1)


[Package RGeostats version 11.1.2 Index]