Page 1 of 1

How to create a lithotype rule for Pluri-Gaussian models?

PostPosted: Sat Jul 12, 2014 9:48 pm
by Nicolas Desassis
The rule (or lithotype rule) is an essential object to perform plurigaussian simulation.
Note that in RGeostats, the categories (also called lithotype, facies, ...) are coded by using integers between 1 and the number of categories. The rules involve one or two Gaussian Random Functions (GRFs) and can be created according to two ways:

By using an interface with the function

Code: Select all
rule.input


or with the command line by using

Code: Select all
rule.create


1) A simple rule

We start this tutorial by showing how to create the simplest (non trivial) rule: a rule involving 1 GRF to model 2 facies.

The command

Code: Select all
rule=rule.input()


provides an interface as follows (the answers are indicated in red and explanations are given hereafter):


Choose the Rule Definition mode
0 : Standard definition of GRF
1 : Definition of one GRF with Shift
2 : Definition of one GRF with Shadow
Rule Definition mode [0,2] : 0
Correlation between the two GRFs (Def=0.000000) [0.000000,1.000000] : 0

Enter the Lithotype Rule (recursively)
Use 'F' (Facies), 'S' (Threshold Y1), 'T' (Threshold Y2) or 'stop'
Choice of node 'main' [F,S,T] : S
Use 'F' (Facies), 'S' (Threshold Y1), 'T' (Threshold Y2) or 'stop'
Choice of node 'S1 [Low]' [F,S,T] : F
Enter the Facies value [1,NA] : 1
Use 'F' (Facies), 'S' (Threshold Y1), 'T' (Threshold Y2) or 'stop'
Choice of node 'S1 [Sup]' [F,S,T] : F
Enter the Facies value [1,NA] : 2


The first answer will always be 0 in this tutorial since the other modes should become obsolete in a future version.
The second one indicates the correlation coefficient between the two GRFs (i.e the value of the cross-covariance at distance 0).
This value has to be consistent with the bivariate spatial model.
However, this is currently irrelevant since we are creating a rule for only one GRF (we have answered 0)
After that, we answer S which indicates that we want to divide the current window in two.
Then we have to answer how to fill the left part of the window. We have answered F which indicates that we want to fill it with a facies and then 1 to say that it is the first one.
Then we are asked to say how we want to fill the right window. We have answered F then 2 to fill it with the facies 2.

To obtain the same result as previously without any interface, we can use:

Code: Select all
rule=rule.create(c("S","F1","F2"))
plot(rule)


Try also:

Code: Select all
rule=rule.create(c("S","F2","F1"))
plot(rule)