Preamble

rm(list=ls())
library(gstlearn)
library(ggplot2)
library(ggpubr)
library(ggnewscale)

flagInternetAvailable=TRUE

Creating the Environment

Defining some essential parameters:

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

## NULL

nbsimu = 20
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 = "Data.NF"
if(flagInternetAvailable){
  download.file(paste0("https://soft.minesparis.psl.eu/gstlearn/data/PluriGaussian/",fileNF), fileNF, quiet=TRUE)
}
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

# TODO: Utiliser CMAP
plot.setDefaultGeographic(dims=c(6,6))
p = ggDefaultGeographic()
p = p + plot(data, nameColor="facies", nameSize="connect")
#         edgecolors='black', sizmin=40, cmap=cmap)
p = p + plot.decoration(title="Conditioning Information")
ggPrint(p)

Proportions and Lithotype Rule

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

# TODO: utiliser cmap
p = ggplot()
p = p + plot(rule, proportions = props, maxG=3)
p = p + plot.decoration(title="Lithotype Rule")
ggPrint(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(dpas=5, npas=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.))

ggplot() + 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 = ggDefaultGeographic()
p1 = p1 + plot(grid, nameRaster="SimuPGS.1")
p1 = p1 + plot(data, nameColor="facies", nameSize="connect")
p2 = ggDefaultGeographic()
p2 = p2 + plot(grid, nameRaster="SimuPGS.2")
p2 = p2 + plot(data, nameColor="facies", nameSize="connect")
p3 = ggDefaultGeographic()
p3 = p3 + plot(grid, nameRaster="SimuPGS.3")
p3 = p3 + plot(data, nameColor="facies", nameSize="connect")
p4 = ggDefaultGeographic()
p4 = p4 + plot(grid, nameRaster="SimuPGS.4")
p4 = p4 + plot(data, nameColor="facies", nameSize="connect")
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, "SimuPGS.1", 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 NA NA 1 
## 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 = paste("SimuPGS.",i, sep="")
  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.1 is invalid
## Simulation SimuPGS.2 is valid
## Simulation SimuPGS.3 is invalid
## Simulation SimuPGS.4 is valid
## Simulation SimuPGS.5 is invalid
## Simulation SimuPGS.6 is invalid
## Simulation SimuPGS.7 is invalid
## Simulation SimuPGS.8 is invalid
## Simulation SimuPGS.9 is invalid
## Simulation SimuPGS.10 is valid
## Simulation SimuPGS.11 is invalid
## Simulation SimuPGS.12 is invalid
## Simulation SimuPGS.13 is invalid
## Simulation SimuPGS.14 is valid
## Simulation SimuPGS.15 is invalid
## Simulation SimuPGS.16 is invalid
## Simulation SimuPGS.17 is invalid
## Simulation SimuPGS.18 is invalid
## Simulation SimuPGS.19 is invalid
## Simulation SimuPGS.20 is invalid

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

## Number of valid simulations = 4 out of 20

Derive the Probability Map

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

## NULL

Display

p = ggDefaultGeographic()
p = p + plot(grid,"Stats.MEAN")
p = p + plot(data,nameColor="facies", nameSize="connect")
p = p + plot.decoration(title="Probability of Connecting Wells")
ggPrint(p)

Tutorial on PluriGaussian Simulations

gstlearn Team

Preamble

Creating the Environment

Variography

Simulations

Simulations under constraints