Preamble

rm(list=ls())
library(gstlearn)

Creating the Environment

Defining some essential parameters:

ndim = 2
defineDefaultSpace(ESpaceType_RN(), ndim);

## NULL

nbsimu = 50
nbcc   = 4
cmap   = c('red', 'blue', 'yellow')

Downloading the data base (from the distribution Exdemo_PGS.db) and creating the output Grid, the Model (Cubic) and the Neighborhood (Unique):

fileNF = loadData("PluriGaussian", "Data.NF")
data = Db_createFromNF(fileNF)
grid = DbGrid_create(nx=c(110,110))
model = Model_createFromParam(type=ECov_CUBIC(), ranges=c(50,30))
neigh = NeighUnique()

Defining Internal Display Function

Defining the internal function plot.dat used to visualize the data set with convenient color representation:

first facies in red
second facies in blue
third facies in yellow
samples which must belong to the same connected component in black

Display Data

p = plot.init(dims=c(6,6)) +
  plot(data, nameColor="facies", nameSize="connect", palette=cmap) +
  plot.decoration(title="Conditioning Information")
plot.end(p)

Proportions and Lithotype Rule

props = dbStatisticsFacies(data)
rule = Rule_createFromNames(c("S","S","F1","F2","F3"))

p = plot.init(dims=c(6,6)) +
  plot(rule, proportions = props, maxG=3, cols=cmap, legendName="Facies") +
  plot.decoration(title="Lithotype Rule")

## Warning in geom_rect(data = df, mapping = aes(xmin = xmin, xmax = xmax, :
## Ignoring unknown parameters: `cols`

plot.end(p)

Variography

Model of Underlying GRF

Calculate the Experimental Variogram of the Underlying Gaussian Random Function and fit the Model (used in PGS).

varioparam = VarioParam_createOmniDirection(dlag=5, nlag=20)
ruleprop = RuleProp_createFromRule(rule, props)
vario = variogram_pgs(data, varioparam, ruleprop)

model_gaus = Model()
err = model_gaus$fit(vario, types=ECov_fromKeys(c("CUBIC")),
                     constraints=Constraints(1.))

plot.init(dims = c(6, 6)) + plot.varmod(vario, model_gaus, asCov = TRUE)

Simulations

PluriGaussian Simulation

err = simpgs(dbin=data, dbout=grid, ruleprop=ruleprop, model1=model_gaus,
             neigh=neigh, nbsimu=nbsimu,
             namconv = NamingConvention("SimuPGS"))

p1 = plot.init(dims=c(6,6)) +
  plot(grid, name=NC_getNameEncoded("SimuPGS", NULL, 1, 1, 1, nbsimu), palette=cmap) +
  plot(data, nameColor="facies", nameSize="connect", palette=cmap)
p2 = plot.init(dims=c(6,6)) +
  plot(grid, name=NC_getNameEncoded("SimuPGS", NULL, 1, 1, 2, nbsimu), palette=cmap) +
  plot(data, nameColor="facies", nameSize="connect", palette=cmap)
p3 = plot.init(dims=c(6,6)) +
  plot(grid, name=NC_getNameEncoded("SimuPGS", NULL, 1, 1, 3, nbsimu), palette=cmap,) +
  plot(data, nameColor="facies", nameSize="connect", palette=cmap)
p4 = plot.init(dims=c(6,6)) +
  plot(grid, name=NC_getNameEncoded("SimuPGS", NULL, 1, 1, 4, nbsimu), palette=cmap) +
  plot(data, nameColor="facies", nameSize="connect", palette=cmap)
ggarrange(p1,p2,p3,p4, nrow=2, ncol=2)

Simulations under constraints

Acceptation Function

Acceptation internal function: Select a Target Facies and build its Connected Components. For each simulation outcome, check the ranks of the connected component(s) at constraining wells and accept the simulation if all ranks are similar.

accept <- function(data, grid, name, verbose=FALSE, transBinary=TRUE, faccc=2)
{
  # Get the indices of samples which should be connected (starting from 0)
  rankData = which(data["connect"] == 1) - 1
  rankGrid = grid$locateDataInGrid(data, rankData)
    if (verbose)
    {
      cat("Number of conditioning data (connected) =",length(rankData),"\n")
    cat("Their ranks in the input Data Base =",rankData,"\n")
    cat("Their ranks in the output Data Base =",rankGrid,"\n")
    }

  # Perform the labelling into connected components
  err = grid$setLocator(name, ELoc_Z(), cleanSameLocator=TRUE)
  err = dbMorpho(grid, EMorpho_CC(), vmin=faccc-0.5, vmax=faccc+0.5)
  cc_list = grid[rankGrid,"Morpho*"]
  if (verbose)
    cat("List of their connected components indices =",cc_list,"\n")

  # Check that the data points belong to the same connected component
  number = length(unique(cc_list))
  retval = (number == 1)
  if (verbose)
    cat("Acceptation score =",retval,"\n")

  # Convert the valid Simulation outcome into a binary image
  if (retval && transBinary)
    grid[name] = (grid["Morpho*"] == cc_list[1])

  grid$deleteColumn("Morpho*")
  retval
}

Experiment the Acceptation Function on one Simulation outcome

isValid = accept(data, grid, NC_getNameEncoded("SimuPGS", NULL, 1, 1, 1, nbsimu), TRUE)

## Number of conditioning data (connected) = 4 
## Their ranks in the input Data Base = 4 12 15 16 
## Their ranks in the output Data Base = 6000 1270 6648 9712 
## List of their connected components indices = 1 1 1 NA 
## Acceptation score = FALSE

cat("Connectivity for Simulation #1 :",isValid,"\n")

## Connectivity for Simulation #1 : FALSE

Probability Map

For each simulation, convert a pixel into 1 if it matches the Target Facies and 0 otherwise
Calculate the mean per pixel over all simulation outcomes

This operation provides the probability that each pixel belongs to the Target Facies, calculated over all simulations that fulfill the Connectivity Constraint.

nb.valid = 0
for (i in 1:nbsimu)
{
  name = NC_getNameEncoded("SimuPGS", NULL, 1, 1, i, nbsimu)
  cat("Simulation",name)
  isValid = accept(data, grid, name)
  if (isValid)
  {
    cat(" is valid\n")
    nb.valid = nb.valid + 1
  }
  else
  {
    cat(" is invalid\n")
    grid$deleteColumn(name)
  }
}

## Simulation SimuPGS.S1 is invalid
## Simulation SimuPGS.S2 is valid
## Simulation SimuPGS.S3 is valid
## Simulation SimuPGS.S4 is invalid
## Simulation SimuPGS.S5 is invalid
## Simulation SimuPGS.S6 is valid
## Simulation SimuPGS.S7 is invalid
## Simulation SimuPGS.S8 is invalid
## Simulation SimuPGS.S9 is invalid
## Simulation SimuPGS.S10 is valid
## Simulation SimuPGS.S11 is invalid
## Simulation SimuPGS.S12 is valid
## Simulation SimuPGS.S13 is invalid
## Simulation SimuPGS.S14 is valid
## Simulation SimuPGS.S15 is valid
## Simulation SimuPGS.S16 is invalid
## Simulation SimuPGS.S17 is invalid
## Simulation SimuPGS.S18 is invalid
## Simulation SimuPGS.S19 is invalid
## Simulation SimuPGS.S20 is valid
## Simulation SimuPGS.S21 is invalid
## Simulation SimuPGS.S22 is invalid
## Simulation SimuPGS.S23 is invalid
## Simulation SimuPGS.S24 is invalid
## Simulation SimuPGS.S25 is invalid
## Simulation SimuPGS.S26 is valid
## Simulation SimuPGS.S27 is invalid
## Simulation SimuPGS.S28 is invalid
## Simulation SimuPGS.S29 is invalid
## Simulation SimuPGS.S30 is invalid
## Simulation SimuPGS.S31 is invalid
## Simulation SimuPGS.S32 is valid
## Simulation SimuPGS.S33 is invalid
## Simulation SimuPGS.S34 is invalid
## Simulation SimuPGS.S35 is valid
## Simulation SimuPGS.S36 is invalid
## Simulation SimuPGS.S37 is valid
## Simulation SimuPGS.S38 is valid
## Simulation SimuPGS.S39 is invalid
## Simulation SimuPGS.S40 is invalid
## Simulation SimuPGS.S41 is invalid
## Simulation SimuPGS.S42 is invalid
## Simulation SimuPGS.S43 is valid
## Simulation SimuPGS.S44 is valid
## Simulation SimuPGS.S45 is invalid
## Simulation SimuPGS.S46 is invalid
## Simulation SimuPGS.S47 is valid
## Simulation SimuPGS.S48 is invalid
## Simulation SimuPGS.S49 is invalid
## Simulation SimuPGS.S50 is invalid

cat("Number of valid simulations =", nb.valid, "out of", nbsimu, "\n")

## Number of valid simulations = 16 out of 50

Derive the Probability Map

grid$statisticsBySample(c("SimuPGS*"),EStatOption_fromKeys(c("MEAN")))

## NULL

Display

p = plot.init(dims=c(6,6)) +
  plot(grid,"Stats.MEAN") +
  plot(data,nameColor="facies", nameSize="connect", palette=cmap) +
  plot.decoration(title="Probability of Connecting Wells")
plot.end(p)

Tutorial on PluriGaussian Simulations

gstlearn Team

Preamble

Creating the Environment

Variography

Simulations

Simulations under constraints