In italo-goncalves/geomod3D: Implicit 3D Geological Modeling and Geostatistics

knitr::opts_chunk$set(echo = T, warning = F, message = F)

Welcome to the geomod3D package. This package was designed to allow easy visualization and manipulation of 3D geospatial data, as well as an open source alternative to implicit 3D geological modeling. This vignette contains a simple example of manipulation and modeling of drillhole data from the Ararangua region in southern Brazil.

Importing drillhole data

There are two kinds of data that are needed to generate the 3D object: collar and assay/lithological data. The data looks like this:

library(geomod3D)
library(rgl)

knitr::kable(ara.collar)
knitr::kable(head(ara.lito))

Let's use the example data to generate a lines3DDataFrame object:

ara.dh <- lines3DDataFrame(collar = ara.collar, assay = ara.lito)
ara.dh

If needed, the coordinates and variables can be recovered with getter functions.

knitr::kable(head(GetCoords(ara.dh, as = "data.frame")))
knitr::kable(head(GetData(ara.dh)))

Subsetting

Any object that extends the spatial3DDataFrame class can be subsetted much like a data frame, with the [ operator. The only recommendation is to use column names instead of numbers.

ara.dh[1:5, ]
ara.dh[, c("HOLEID", "Form")]

A single column can be extracted directly with the [[ operator:

unique(ara.dh[["Form"]])

Assignment can also be done the same way as in a data frame:

ara.dh[, "Length.feet"] <- 0.3 * ara.dh[["Length"]]

Visualization

The next step is visualization. The Draw* functions add graphical objects to the current rgl window. If there is none, one will be created.

# Color coding for the geological formations
formcolor <- c("green","cornflowerblue","burlywood4","firebrick",
               "darkorchid","khaki")
formval <- c("Recente", "Rio do Rasto", "Estrada Nova", "Irati",
             "Palermo", "Rio Bonito")

# Opening rgl window to display in html page
invisible(open3d(useNULL = T, windowRect = c(0,0,800,800), zoom = 1))

# Drawing holes
DrawDrillholes(ara.dh, by = "Form", values = formval, col = formcolor,
               size = 100)

# Hole ID
DrawHoleID(ara.dh)

# The aspect3d function can be used to control the vertical exaggeration
aspect3d(1.731788, 1, 0.2)

# Axes
axes3d()

# Showing scene
rglwidget(width = 800, height = 600, elementId = "holes", reuse = NA)

Modeling

In order to make an implicit model, first the data is converted to a suitable format, then we choose a set of parameters and build a Gaussian Process model.

Preparing the data

The drillhole data is processed to extract points inside each geological formation and also in the boundaries between formations. The package contains a series of methods to perform this task.

As we are modeling the geological formation, the other variables can be discarded and adjacent segments of the same formation can be merged in order to reduce the size of the dataset.

ara.dh.proc <- MergeSegments(ara.dh, by = "Form", keep = "Form")
ara.dh.proc

Next we extract points containing the geological information. As a point can belong to a geological contact, each point carries two labels. If the labels differ, then the point lies in a contact.

# Duplication of labels for the line segments
ara.dh.proc[, "Form.down"] <- ara.dh.proc[, "Form.up"] <- ara.dh.proc[, "Form"]

# Conversion to point data
ara.point <- Pointify(ara.dh.proc, locations = seq(0.1, 0.9, 0.2))
ara.point

# Finding the contacts
ara.contacts <- GetContacts(ara.dh.proc, by = "Form")
ara.contacts

# Joining the two objects
ara.point <- Bind(ara.point, ara.contacts)
ara.point

The implicit model

The first step in implciit modeling is the definition of the best covariance model. As the geological layers seem to be stratified, the range in the vertical direction is made smaller than in the horizontal. Due to the sparse data, the same model will be used for all geological classes, and the parameters will be found using an exaustive search. The best fit criterion is the model's log-likelihood.

amplitude = seq(0.1, 1, 0.05)
range = seq(5000, 10000, 500)
nugget = seq(0.025, 0.3, 0.025)
param_df <- expand.grid(amplitude = amplitude, range = range, nugget = nugget)
param_df$logLik <- NA

for (i in seq(nrow(param_df))){
  cov_model <- covarianceStructure3D(type = "cubic", 
                                   contribution = param_df$amplitude[i],
                                   maxrange = param_df$range[i],
                                   midrange = param_df$range[i],
                                   minrange = 0.05 * param_df$range[i])
  gp.model <- GP_geomod(data = ara.point, 
                        value1 = "Form.up", 
                        value2 = "Form.down",
                        model = cov_model, 
                        nugget = param_df$nugget[i])
  param_df$logLik[i] <- logLik(gp.model)
}

# best parameters
param_best <- param_df[which.max(param_df$logLik), ]
param_best

Next we define a point grid over which to make the predictions:

point.grid <- grid3DDataFrame(gridx = seq(647000, 661000, 400),
                              gridy = seq(6791000, 6799000, 400),
                              gridz = seq(-550, 50, 12.5), 
                              fields = "Form")

Then we build a Gaussian Process model and use it to make predictions:

cov_model <- covarianceStructure3D(type = "cubic", 
                                   contribution = param_best$amplitude,
                                   maxrange = param_best$range,
                                   midrange = param_best$range,
                                   minrange = 0.05 * param_best$range)

gp.model <- GP_geomod(data = ara.point, 
                      value1 = "Form.up", 
                      value2 = "Form.down",
                      model = cov_model, 
                      nugget = param_best$nugget)

set.seed(1234) # for reproducibility
point.grid <- Predict(gp.model, target = point.grid, to = "Form")
point.grid

Model visualization

The geological contacts can be drawn by transforming each indicator variable into a 3D array and feeding it to the misc3d::contour3d() function. Note that the model matches the geological contacts in the data.

library(misc3d)

# Opening rgl window to display in html page
invisible(open3d(useNULL = T, windowRect = c(0,0,800,800), zoom = 1))

# Drawing holes
DrawDrillholes(ara.dh, by = "Form", values = formval, col = formcolor,
               size = 100)

# Contouring each indicator variable
form <- Make3DArray(point.grid, "Form..Recente.ind")
display1 <- contour3d(form$value, level = 0, 
                      x = form$x, y = form$y, z = form$z, 
                      color = "green", alpha = 1, add = T)
form <- Make3DArray(point.grid, "Form..Rio.do.Rasto.ind")
display2 <- contour3d(form$value, level = 0, 
                      x = form$x, y = form$y, z = form$z, 
                      color = "cornflowerblue", alpha = 1, add = T)
form <- Make3DArray(point.grid, "Form..Estrada.Nova.ind")
display3 <- contour3d(form$value, level = 0, 
                      x = form$x, y = form$y, z = form$z, 
                      color = "burlywood4", alpha = 1, add = T)
form <- Make3DArray(point.grid, "Form..Irati.ind")
display4 <- contour3d(form$value, level = 0, 
                      x = form$x, y = form$y, z = form$z, 
                      color = "firebrick", alpha = 1, add = T)
form <- Make3DArray(point.grid, "Form..Palermo.ind")
display5 <- contour3d(form$value, level = 0, 
                      x = form$x, y = form$y, z = form$z, 
                      color = "darkorchid", alpha = 1, add = T)
form <- Make3DArray(point.grid, "Form..Rio.Bonito.ind")
display6 <- contour3d(form$value, level = 0, 
                      x = form$x, y = form$y, z = form$z, 
                      color = "khaki", alpha = 1, add = T)

# Aspect ratio and axes
aspect3d(1.731788, 1, 0.2)
axes3d()

# Display scene
rglwidget(width = 800, height = 600, elementId = "model1") %>%
  toggleWidget(ids = display1, label = "Recente") %>%
  toggleWidget(ids = display2, label = "Rio do Rasto") %>%
  toggleWidget(ids = display3, label = "Estrada Nova") %>%
  toggleWidget(ids = display4, label = "Irati") %>%
  toggleWidget(ids = display5, label = "Palermo") %>%
  toggleWidget(ids = display6, label = "Rio Bonito")

The geological layers appear closed away from the data. This happens because the model steadily loses confidence as we move away from the data points. Grid nodes in those regions receive the label "Unknown".

The entropy can be used as another measure of uncertainty:

# Opening rgl window to display in html page
invisible(open3d(useNULL = T, windowRect = c(0,0,800,800), zoom = 1))

# Drawing holes
DrawDrillholes(ara.dh, by = "Form", values = formval, col = formcolor,
               size = 100)

entropy <- Make3DArray(point.grid, "Form..Entropy")
e1 <- contour3d(entropy$value, level = quantile(entropy$value, probs = 0.1), 
                      x = entropy$x, y = entropy$y, z = entropy$z, 
                      color = "green", alpha = 0.5, add = T)
e2 <- contour3d(entropy$value, level = quantile(entropy$value, probs = 0.5), 
                      x = entropy$x, y = entropy$y, z = entropy$z, 
                      color = "yellow", alpha = 0.5, add = T)
e3 <- contour3d(entropy$value, level = quantile(entropy$value, probs = 0.9), 
                      x = entropy$x, y = entropy$y, z = entropy$z, 
                      color = "red", alpha = 0.5, add = T)

# Aspect ratio and axes
aspect3d(1.731788, 1, 0.2)
axes3d()

rglwidget(width = 800, height = 600, elementId = "entropy") %>%
  toggleWidget(ids = e1, label = "Low") %>%
  toggleWidget(ids = e2, label = "Mid") %>%
  toggleWidget(ids = e3, label = "High")