We define the Space dimension

defineDefaultSpace(ESpaceType_RN(), 2)
## NULL

Creating a dummy Model used for simulating a random field

mymodel = Model_createFromParam(ECov_CUBIC(), range=10, sill=1)

Grid representations

Standard Grid

We construct a standard non rotated 2-D grid with two simulated variables

nx = c(70,25)
dx = c(1,2)
x0 = c(-40, 20)
mygrid = DbGrid_create(nx,dx,x0)

err = simtub(NULL,mygrid,mymodel,nbsimu=2)
mygrid$display()
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 5
## Total number of samples      = 1750
## 
## Grid characteristics:
## ---------------------
## Origin :    -40.000    20.000
## Mesh   :      1.000     2.000
## Number :         70        25
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = Simu.1 - Locator = z1
## Column = 4 - Name = Simu.2 - Locator = z2
## NULL

We construct a Selection

mygrid["sel"] = 1. - 
  (mygrid["x1"] > 0) * (mygrid["x1"] < 15) * (mygrid["x2"] > 40) * (mygrid["x2"] < 50)
err = mygrid$setLocator("sel",ELoc_SEL())
mygrid
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 6
## Total number of samples      = 1750
## Number of active samples     = 1694
## 
## Grid characteristics:
## ---------------------
## Origin :    -40.000    20.000
## Mesh   :      1.000     2.000
## Number :         70        25
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = Simu.1 - Locator = z1
## Column = 4 - Name = Simu.2 - Locator = z2
## Column = 5 - Name = sel - Locator = sel

Displaying the grid cells

p = ggDefaultGeographic()
p = p + plot.grid(mygrid, nameRaster="Simu.1", nameContour="Simu.2", 
              flagLegendRaster=TRUE, legendNameRaster="ma Legende")
p = p + plot.decoration(title="Display of Grid Cells")
ggPrint(p)

Using another collor scale. Among the different choices, one makes it possible to create a color scale on the fly, specifying the starting (yellow), middle (red) and ending (blue) colors. The color referred to as ‘naColor’ serves for encoding pixels whose values have been left undefined.

ggDefaultGeographic() + plot.grid(mygrid, palette=c("yellow","red","blue"), naColor="white")

Displaying the grid nodes only

p = ggDefaultGeographic()
p = p + plot.point(mygrid, nameColor="Simu.1",nameSize="Simu.2", 
                   flagLegendColor = TRUE,
                   legendNameColor="myColor")
p = p + plot.decoration(title="Display of Grid Nodes")
ggPrint(p)

Rotated grid

mygrid = DbGrid_create(nx,dx,x0,angles=c(10,0))
err = simtub(NULL,mygrid,mymodel,nbsimu=2)
mygrid
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 5
## Total number of samples      = 1750
## 
## Grid characteristics:
## ---------------------
## Origin :    -40.000    20.000
## Mesh   :      1.000     2.000
## Number :         70        25
## Rotation Angles        =     10.000     0.000
## Direct Rotation Matrix
##                [,  0]    [,  1]
##      [  0,]     0.985    -0.174
##      [  1,]     0.174     0.985
## Inverse Rotation Matrix
##                [,  0]    [,  1]
##      [  0,]     0.985     0.174
##      [  1,]    -0.174     0.985
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = Simu.1 - Locator = z1
## Column = 4 - Name = Simu.2 - Locator = z2

Displaying the cell contents in the rotated grid

p = ggDefaultGeographic()
p = p + plot.grid(mygrid, nameRaster="Simu.1", flagLegendRaster=FALSE)
p = p + plot.decoration(title="Display of Rotated Grid")
ggPrint(p)

As a set of grid nodes

p = ggDefaultGeographic()
p = p + plot.point(mygrid,nameColor="Simu.1", flagLegendColor = TRUE)
p = p + plot.decoration(title="Display of Rotated Grid Nodes")
ggPrint(p)

Points and Polygon

A set of points is sampled from the previous Grid and stored in a new Point Db. The number of samples if fixed to 1% of the number of grid nodes.

mypoint = Db_createSamplingDb(mygrid,0.01)
mypoint$display()
## 
## Data Base Characteristics
## =========================
## 
## Data Base Summary
## -----------------
## File is organized as a set of isolated points
## Space dimension              = 2
## Number of Columns            = 6
## Total number of samples      = 17
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = rank.1 - Locator = NA
## Column = 2 - Name = x1 - Locator = x1
## Column = 3 - Name = x2 - Locator = x2
## Column = 4 - Name = Simu.1 - Locator = z1
## Column = 5 - Name = Simu.2 - Locator = z2
## NULL

We create a polygon as the convex hull of the samples

mypoly = Polygons_createFromDb(mypoint)

We now display the points and the polygon on top of the grid: the overlay is ensured by using the argument ‘p’.

p = ggDefaultGeographic()
p = p + plot(mygrid)
p = p + plot(mypoly, fill=NA, color='yellow', linewidth=1)
p = p + plot(mypoint,color="black")
p = p + plot.decoration(title="mon titre", xlab="mon axe des X", ylab="Mon axe des Y")
ggPrint(p)

Several plots on the same figure

We create two layers containing the Point and the Grid representation. We then use ggarrange to display them side-by-side.

p1 = ggDefaultGeographic() + plot(mygrid)
p2 = ggDefaultGeographic() + plot(mypoint)
ggarrange(p1, p2, labels = c("Plot #1", "Plot #2"), ncol = 2, nrow = 1)

Variograms and Models

We calculate the variogram along the two main directions of the 2-D grid. We compute it for the 2 variables currently defined in the Grid file.

varioparam = VarioParam_createMultipleFromGrid(mygrid,npas=10)
vario = Vario(varioparam)
err = vario$compute(mygrid,ECalcVario_VARIOGRAM())

We fit a Model (automatic procedure)

model = Model()
err = model$fit(vario,type=ECov_fromKeys(c("SPHERICAL", "CUBIC")))

We display the experimental variogram for the first variable in the first direction.

ggplot() + plot(model)

p = ggDefault()
p = p + plot(vario, ivar=0, jvar=0, idir=0, drawPsize=-1, color="red")
p = p + plot.decoration(title="First Variable - First Direction")
ggPrint(p)

We display the experimental variogram for the first variable in the second direction.

p = ggDefault()
p = p + plot(vario, ivar=0, jvar=0, idir=1)
p = p + plot.decoration(title="First Variable - Second Direction")
ggPrint(p)

We simply overlay the contents of the two previous plots.

p = ggDefault()
p = p + plot(vario, ivar=0, jvar=0, idir=0, color="black")
p = p + plot(vario, ivar=0, jvar=0, idir=1, color='red')
p = p + plot.decoration(title="First Variable - Two directions")
ggPrint(p)

We display the cross-variogram between both variables along both directions.

p = ggDefault()
p = p + plot(vario, ivar=1, jvar=0, idir=-1)
p = p + plot.decoration(title="Cross-variogram - All directions")
ggPrint(p)

We display the simple and cross variograms in both calculation directions.

p = ggDefault()
p = p + plot(vario, ivar=-1, jvar=-1, idir=-1)
p = p + plot.decoration(title="Simple and Cross Variograms in all directions")
ggPrint(p)

We display the Model as calculated in the first direction of the variogram

p = ggDefault()
p = p + plot(model, ivar=1, jvar=1, vario=vario, idir=0)
p = p + plot.decoration(title="Model for the Second Variable in First Direction")
ggPrint(p)

We now represent all the simple and cross variograms together with the fitted models.

p = ggDefault()
p = p + plot.varmod(vario, model)
p = p + plot.decoration(title="All variograms in First Direction")
ggPrint(p)

multi.varmod(vario=vario, model=model)

Slice of a 3-D Grid

g3D = DbGrid_create(nx=c(10,15,20))
g3D
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 3
## Number of Columns            = 4
## Total number of samples      = 3000
## 
## Grid characteristics:
## ---------------------
## Origin :      0.000     0.000     0.000
## Mesh   :      1.000     1.000     1.000
## Number :         10        15        20
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = x3 - Locator = x3

Plot a horizontal slice:by default, this corresponds to the first level IZ=1

p = ggDefault()
p = p + plot.grid(g3D, "rank", flagLegendRaster=TRUE, legendNameRaster="Rank")
p = p + plot.decoration(title="Horizontal Slice (IZ=1)")
ggPrint(p)

Plot another horizontal slice (IZ=10): the color scale is adapter to the values (ranks) of the new set of cells to be represented

p = ggDefault()
p = p + plot.grid(g3D, nameRaster="rank", flagLegendRaster=TRUE, legendNameRaster="Rank", corner=c(0,0,10))
p = p + plot.decoration(title="Horizontal Slice (IZ=10)")
ggPrint(p)

Display a vertical slice (XOZ) for the index IY=3

p = ggDefault()
p = p + plot.grid(g3D, nameRaster="rank", flagLegendRaster=TRUE, legendNameRaster="Rank", 
                  posY=2, corner=c(0,3,0))
p = p + plot.decoration(title="Vertical Slice XOZ (IY=3)")
ggPrint(p)