--- title: "PluriGaussian Simulations" author: "D. Renard" date: "26 july 2017" output: pdf_document: toc: yes html_document: fig_caption: yes highlight: tango toc: yes vignette: > %\VignetteIndexEntry{PluriGaussian Simulations} %\VignetteEngine{knitr::rmarkdown} \usepackage[utf8]{inputenc} --- ```{r Loading_library,include=FALSE} library(RGeostats) rm(list=ls()) ``` # Introduction The ingredients for a plurigaussian simulation are: • The number of lithotypes to be simulated • Their proportions - constant in the domain for a stationary case or variable in space for a non stationary case • Their relationships given by the lithotype rule which involves one or two underlying Gaussian Random Function (GRF) • The variogram models of the underlying GRF(s) When all these parameters are defined, the non-conditional simulations can be performed using the classical parameters such as the seed, the number of bands (when the GRF simulation is performed using the Turning Bands method) and the number of simulations. Finally all these simulations can match the (Gaussian) data when available: they are then called conditional simulations. This feature is not addressed in this case study. ## Define the simulation grid The non-conditional simulations will be stored in a Db (called outgrid1) organized as a 2-D grid with 200 nodes along X and 100 nodes along Y (square mesh of 1m side) and which is created in the next command: ```{r Creating_Grid} outgrid1<-db.create(flag.grid=T,x0=c(0,0),nx=c(200,100),dx=c(1,1)) ``` ## Number of lithotypes to be simulated For this illustration we have chosen four lithotypes. Their colors are given in the following colormap: ```{r Defining_color_map} pal4fac <- c("red","orange","green","blue") ``` ## Lithotype Rule The relationships between lithotypes are summarized in a lithotype rule. It consists of a 2-D diagram (represented as a square) where the horizontal axis gives the range of values for the first GRF and the vertical axis for the second GRF. Both GRF vary between[). This square is subdivided into as many rectangles as there are lithotypes. The procedure rule.input allows the user to define the lithotype rule interactively. The first question concerns the simulation method: • The Plurigaussian simulation: this corresponds to the current illustration, where the two underlying GRFs are not correlated • The Shifted simulation • The Shadow simulation The rule construction is performed sequentially, describing the rectangles starting from the lower left corner. The rule is read from left to right and from bottom to top. Each rectangle is defined by the ordered thresholds. Therefore S the user must choose which threshold is encountered: (if the next limit is T defined along the first GRF) or (if the next limit is defined along the second F GRF). When the rectangle is defined, the user chooses to define the facies number. ```{r Defining_Lithotype_Rule} rule4fac = rule.create(c("S","F1","T","F2","S","F3","F4")) ``` The visualization of the lithotype rule is obtained using the following command: ```{r Visualizing_Lithotype_Rule} plot (rule4fac,col=pal4fac) ``` ## Variogram Models In the scope of the Plurigaussian model, the lithotypes are obtained by thresholding the underlying GRFs. Consequently, the variograms of the lithotype indicators are related to the models of these GRFs, taking the proportions and the lithotype rule into account. The fitting procedure is a trial and error method: we postulate the model of the underlying GRFs and check that the quality of the t for the variograms of the corresponding lithotype indicators. Here, the model is assumed to be known and is defined interactively. The model of each underlying GRF is defined using the command model.input. The model for the first GRF is an anisotropic Cubic structure and an isotropic Exponential structure for the second GRF. The models for both underlying GRFs are stored in the file modY1 and modY2. ```{r} modY1 = model.create(vartype="Cubic",sill=1,aniso.angles=-60,range=c(35,20)) modY2 = model.create(vartype="Exponential",sill=1,range=15) ``` We check the contents of these models: ```{r Checking_First_Model} modY1 ``` ```{r Checking_Second_Model} modY2 ``` # Stationary PluriGaussian Simulation The simulation is performed and stored in the output grid. In the following example, a single non-conditional simulation is performed with: • the lithotype rule (rule4fac) • the models of the underlying GRFs (modY1 and modY2) • the constant proportions: here, they are considered as all equal to 0.25 (as they must sum up to 1) • 600 turning bands and a seed equal to 126543. ```{r Performing_PGS} outgrid1 <- simpgs(dbout=outgrid1,rule=rule4fac, model1=modY1,model2=modY2,props=c(0.25,0.25,0.25,0.25), nbsimu=1,nbtuba=600,seed=126543) ``` The display of the simulated outcome is obtained with the following command: ```{r Representing_Stationary_PGS_Simulation} plot(outgrid1,col=pal4fac,flag.scale=TRUE, title="Stationary PluriGaussian Simulation") ``` # Non-stationary PluriGaussian Simulation For illustrating the non-stationary case, a second output grid (called outgrid2) is created. ```{r Creating_second_Grid} outgrid2<-db.create(flag.grid=TRUE,x0=c(0,0),nx=c(200,100),dx=c(1,1)) ``` The proportions of each lithotype vary in the space: they are stored from the grid. In our example, four variables are defined (called Prop1, Prop2, Prop3 and Prop4) which contain the proportions. For instance: ```{r Creating_proportions} outgrid2 <- db.add(outgrid2,Prop1=x1/200) outgrid2 <- db.add(outgrid2,Prop2=0.1*(1-Prop1)) outgrid2 <- db.add(outgrid2,Prop3=0.5*(1-Prop1-Prop2)) outgrid2 <- db.add(outgrid2,Prop4=1-Prop1-Prop2-Prop3) ``` The type (locators) of these variables (located in fields 4 to 7) must be defined as «p» : ```{r Defining_Proportion_Locators} outgrid2 <- db.locate(outgrid2,seq(4,7),loctype="p") outgrid2 ``` The different proportions can be displayed using: ```{r Displaying_Proportions} plot(outgrid2,name="Prop1",pos.legend=1, title="Proportion of First Lithotype") plot(outgrid2,name="Prop4",pos.legend=1, title="Proportion of fourth Lithotype") ``` The rule and the models for the two underlying GRF are kept the same as in the stationary case. The command for the non stationary simulation is: ```{r Non_stationary_PluriGaussian_Simulation} outgrid2 <- simpgs(dbout=outgrid2, rule=rule4fac, dbprop=outgrid2, model1=modY1,model2=modY2, nbsimu=1, nbtuba=600,seed=126543,verbose=TRUE) ``` The corresponding result is displayed in the next figure: ```{r Representing_Non_Stationary_PGS_Simulation} plot(outgrid2,col=pal4fac,flag.scale=TRUE, title="Non Stationary PluriGaussian Simulation") ```