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)