Preamble

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

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 = 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:

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 = NA 1 1 1 
## Acceptation score = FALSE
cat("Connectivity for Simulation #1 :",isValid,"\n")
## Connectivity for Simulation #1 : FALSE

Probability Map

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 valid
## Simulation SimuPGS.4 is invalid
## Simulation SimuPGS.5 is valid
## Simulation SimuPGS.6 is valid
## Simulation SimuPGS.7 is valid
## Simulation SimuPGS.8 is valid
## Simulation SimuPGS.9 is valid
## Simulation SimuPGS.10 is invalid
## Simulation SimuPGS.11 is valid
## Simulation SimuPGS.12 is valid
## Simulation SimuPGS.13 is valid
## Simulation SimuPGS.14 is invalid
## Simulation SimuPGS.15 is valid
## Simulation SimuPGS.16 is valid
## Simulation SimuPGS.17 is valid
## Simulation SimuPGS.18 is invalid
## Simulation SimuPGS.19 is valid
## Simulation SimuPGS.20 is valid
cat("Number of valid simulations =", nb.valid, "out of", nbsimu, "\n")
## Number of valid simulations = 15 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)