Introduction

This R Markdown document tests the utilities developed to compute the layer proportions for each selected cells of a three dimensional grid. The $N$ layers are defined by the elevation of their top limit.

knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
library(gstlearn)

## 
## Attaching package: 'gstlearn'

## The following objects are masked from 'package:base':
## 
##     message, toString

OptCst_defineByKey("ASP",0)

## NULL

library(ggplot2)
library(ggpubr)
library(ggnewscale)
set.seed(43243)
opers = EStatOption_fromKeys(c("NUM","MINI","MAXI","MEAN","STDV"))

Parameters

# The 3D grid is defined on  [0,1]^3 
# 2D limits of the 3D grid 
limits <- list(XP = 1.0*c(0, 1, 1, 0, 0),
               YP = 1.0*c(0, 0, 1, 1, 0))
p_lim   = PolyElem(limits$XP, limits$YP)
pol_lim = Polygons()
err     = pol_lim$addPolyElem(p_lim)

flag.rotation = TRUE
# ----------------------------------------------
# ----------------------------------------------
if(flag.rotation) {
# 3D grid
nx = c(10, 13, 7)
dx = 1/nx
x0 = dx/2
# 2D grid
nx_2d = c(256, 256)
dx_2d = sqrt(2.0) / nx_2d
x0_2d = c(1/2, -1/2)
angles_2d = c(45, 0) # in degrees
} else { # no rotation
# 3D grid
nx = c(10, 10, 10)
dx = 1/nx
x0 = dx/2
# 2D grid
nx_2d = nx[1:2]
dx_2d = dx[1:2]
x0_2d = dx_2d/2
angles_2d = c(0, 0) # in degrees
}

grid = DbGrid_create(nx = nx, dx = dx, x0 = x0)
surfaces = DbGrid_create(nx = nx_2d, dx = dx_2d, x0 = x0_2d, angles = angles_2d)

# plots
p = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "rank", useSel = TRUE, flagLegendRaster = FALSE) +
    plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
    plot.decoration(xlab = "Easting", ylab = "Northing", title = paste0("3D grid and limits"))
ggPrint(p)

p = ggDefaultGeographic() + 
    plot.grid(surfaces, nameRaster = "rank", useSel = TRUE, flagLegendRaster = FALSE) +
    plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
    plot.decoration(xlab = "Easting", ylab = "Northing", title = paste0("2D grid and limits"))
ggPrint(p)

# ----------------------
# selection of some grid cells
# ----------------------
sel  = (grid["x1"] >= 0.20)&(grid["x1"] <= 0.70) & 
       (grid["x3"] >= 0.10)&(grid["x3"] <= 0.70)
sel  = as.numeric(sel)
err  = grid$setColumn(tab = sel, name = "sel")
err  = grid$setLocator(name = "sel", ELoc_SEL(), cleanSameLocator = TRUE)

# ----------------------
# plot of a slice
# ----------------------
var_nm = "rank"
posx = 1; xlab = "Ox"
posy = 3; ylab = "Oz"
k = 3
corner = rep(k - 1, 3); 
corner[posx] = 0; 
corner[posy] = 0
slice = DbGrid_create(nx = c(nx[posx], nx[posy]), dx = c(dx[posx], dx[posy]), x0 = c(x0[posx], x0[posy]))
val = grid$getOneSlice(name = var_nm, posx = posx - 1, posy = posy - 1 , corner = c(0,0,k-1), 
                       useSel = TRUE)
err   = slice$setColumn(tab = val, name = paste0("K = ", k))
p = ggDefaultGeographic() + 
    plot.grid(slice, nameRaster = paste0("K = ", k), useSel = TRUE, legendNameRaster = var_nm,
            flagLegendRaster = TRUE) +
    plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
    plot.decoration(xlab = xlab, ylab = ylab, 
                  title = paste0("Slice K = ", k))
ggPrint(p)

# -----------------------------------
# simulation of the surfaces
# -----------------------------------
N = 3
h_mean = c(0.5, 0.25, 0.25)
nbsim = c(3, 4, 3)

# N = 2
# h_mean = c(0.75, 0.25)
# nbsim = c(3, 4)

model = Model_createFromParam(ECov_GAUSSIAN(), range=0.5, sill=1.0)
for (i in 1:N) {
  err = simtub(dbin = NULL, dbout = surfaces, model = model, nbsimu = nbsim[i], namconv = NamingConvention(paste0("Y_", i)))
}

# conversion into uniform thickness
for (i in 1:N) {
  val = matrix(surfaces$getColumns(names = paste0("Y_", i , ".", 1:nbsim[i]), useSel = FALSE),
               nrow = surfaces$getSampleNumber(useSel = FALSE), ncol = nbsim[i])
  val = h_mean[i] * pnorm(val)
  for (k in 1:nbsim[i]) {
    err = surfaces$setColumn(tab = val[, k], name = paste0("Z_", i , ".", k), useSel = FALSE)
  }
}

# plot of a variable
var_nm = "Z_1.3"
p = ggDefaultGeographic() + 
    plot.grid(surfaces, nameRaster = var_nm, useSel = TRUE, legendNameRaster = var_nm,
            flagLegendRaster = TRUE) +
    plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
    plot.decoration(xlab = "Easting", ylab = "Northing", 
                  title = paste0("2D Surface: variable = ", var_nm))
ggPrint(p)

Computing the proportions

The simulated variables are the thickness $Z_i^{(k)}(s)$.

Sample selection

For a block selected by $(i, j, k)$ the samples $(x,y,z)$ of the surfaces are selected by:

$|x - (x_0[1]+ dx[1] (i - 1)| dx[1] /2 $
$|y - (x_0[2]+ dx[2] (j - 1)| dx[2] /2 $

# compute the rank from ijk of the values of an array defined by its dimensions nn
# This is a vector version of ijk2rank and rank2ijk for any number of dimensions
#  rank = i + nn[1] (j - 1) + nn[1]*nn[2]*(k - 1) + ...
rank2ijkVec <- function(rank, nn) {
    stopifnot((range(rank)[1] > 0)&(range(rank)[2] <= prod(nn)))
    N = length(nn)
    IJK = matrix(NaN, nrow = length(rank), ncol = N)
    rr = rank - 1
    for (i in 1:N) {
      IJK[,i] = 1 + rr %% nn[i]
      if(i < N) {rr = (rr - (IJK[,i] - 1)) %/% nn[i]}
    }
    if(length(rank) == 1) {IJK = as.numeric(IJK)}
    IJK
  }
  
# convert the ijk corresponding to the rank of an array defined by its dimensions nn
ijk2rankVec <- function(ijk, nn) {
    if(!is.matrix(ijk)){ijk = matrix(ijk, nrow = 1, ncol = length(ijk))}
    stopifnot(dim(ijk)[2] == length(nn))
    as.numeric(1 + (ijk-1) %*% cumprod(c(1, nn[-length(nn)])))
  }

# test on a 4d array  
nn = c(5, 13, 7, 31) 
rank = 1:prod(nn)
ijk  = rank2ijkVec(rank = rank, nn = nn)
stopifnot(max(abs(rank - ijk2rankVec(rank2ijkVec(rank,nn),nn))) == 0)
stopifnot(max(abs(ijk - rank2ijkVec(ijk2rankVec(ijk,nn),nn))) == 0)

# test on a cell
cell = c(3, 4, 5)
rank = ijk2rankVec(ijk = cell, nn = grid$getNXs())
stopifnot(rank == ijk2rankVec(ijk = rank2ijkVec(rank = rank, nn = grid$getNXs()), nn = grid$getNXs()))
stopifnot(max(abs(cell - rank2ijkVec(rank = ijk2rankVec(ijk = cell, nn = grid$getNXs()), nn = grid$getNXs()))) == 0)

# select the samples in dbin belonging to the grid cell and return the number of selected samples
# the selection works only with "not rotated" grid 
select_samples <- function(dbin, rank, grid, nameSel = "inBlock") {
  stopifnot(max(abs(grid$getAngles())) == 0) # no rotation
  cell = rank2ijkVec(rank = rank, nn = grid$getNXs())
  x0  = grid$getX0s()
  dx  = grid$getDXs()
  xCell = x0 + dx * (cell - 1)
  u1  = dbin$getColumn(name = "x1", useSel = FALSE)
  u2  = dbin$getColumn(name = "x2", useSel = FALSE)
  sel = as.numeric( (abs(u1 - xCell[1]) <= dx[1]/2)&(abs(u2 - xCell[2]) <= dx[2]/2))
  err = dbin$deleteColumn(name = nameSel)
  err = dbin$setColumn(tab = sel, name = nameSel, useSel=FALSE)
  err = dbin$setLocator(name = nameSel, locatorType = ELoc_SEL(), cleanSameLocator = TRUE)
  dbin$getSampleNumber(useSel = TRUE)
}

err = surfaces$deleteColumn("inBlock")
err = surfaces$deleteColumn("sel_*")
print(paste0(">>> number of selected samples = ", select_samples(surfaces, rank = ijk2rankVec(ijk = cell, nn = grid$getNXs()), grid)))

## [1] ">>> number of selected samples = 238"

# plot of selected samples of the surface grid

err = surfaces$clearSelection()
err = surfaces$setLocator(name = "inBlock", locatorType = ELoc_SEL(), cleanSameLocator = TRUE)
var_nm = "Z_1.3"
p = ggDefaultGeographic() + 
    plot.grid(surfaces, nameRaster = var_nm, useSel = TRUE, legendNameRaster = var_nm,
            flagLegendRaster = TRUE) +
    plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
    plot.decoration(xlab = "Easting", ylab = "Northing", 
                  title = paste0("2D surface: variable = ", var_nm))
ggPrint(p)

Simulation Post-Processing for calculating Layer Proportions

Xavier Freulon

2023-07-24

Introduction

Parameters

Computing the proportions

Sample selection