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"))
# 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)
The simulated variables are the thickness \(Z_i^{(k)}(s)\).
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)