recexcavAAR: Vignette >>Ifri el Baroud<<

Clemens Schmid

July 2016

First we need to load recexcavAAR and some external packages:

library(devtools)
devtools::load_all() # in your script: library(recexcavAAR)
library(dplyr)
library(kriging)
library(magrittr)
library(rgl)

Now let’s imagine an artificial and pretty simple excavation trench with a depth of 2 meters, a length of 3 meters and a width of 1 meter.

edges <- data.frame(
  x = c(0, 3, 0, 3, 0, 3, 0, 3),
  y = c(0, 0, 0, 0, 1, 1, 1, 1),
  z = c(0, 0, 2, 2, 0, 0, 2, 2)
)

We can plot the corner points of this trench with rgl::plot3d:

# avoid plotting in X11 window
open3d(useNULL = TRUE)
## null 
##    1
plot3d(
  edges$x, edges$y, edges$z,
  type="s",
  aspect = c(3, 1, 2),
  xlab = "x", ylab = "y", zlab = "z",
  sub = "Grab me and rotate me!"
)

bbox3d(
  xat = c(0, 1, 2, 3),
  yat = c(0, 0.5, 1),
  zat = c(0, 0.5, 1, 1.5, 2)
)

rglwidget()

When we look at the profiles of our fictional trench we can trace three clearly separated horizons following the natural slope. Let’s figuratively take some measurements of the two horizon borders by creating two data.frames df1 and df2 with points along the profiles. The z axis values are randomly computed.

df1 <- data.frame(
  x = c(rep(0, 6), seq(0.2, 2.8, 0.2), seq(0.2, 2.8, 0.2), rep(3,6)),
  y = c(seq(0, 1, 0.2), rep(0, 14), rep(1, 14), seq(0, 1, 0.2)), 
  z = c(seq(0.95, 1.2, 0.05), 0.9+0.05*rnorm(14), 1.3+0.05*rnorm(14), seq(0.95, 1.2, 0.05))
)

df2 <- data.frame(
  x = c(rep(0, 6), seq(0.2, 2.8, 0.2), seq(0.2, 2.8, 0.2), rep(3,6)),
  y = c(seq(0, 1, 0.2), rep(0, 14), rep(1, 14), seq(0, 1, 0.2)),
  z = c(seq(0.65, 0.9, 0.05), 0.6+0.05*rnorm(14), 1.0+0.05*rnorm(14), seq(0.65, 0.9, 0.05))
)

Looks complicated? Becomes pretty simple when we look at an other plot. For this one we add the points to the previous plot object by calling rgl::points3d.

points3d(
  df1$x, df1$y, df1$z,
  col = "darkgreen",
  add = TRUE
)

points3d(
  df2$x, df2$y, df2$z,
  col = "blue",
  add = TRUE
)

rglwidget()

We can put this two or even more data.frames df1 and df2 into a list lpoints and feed it to recexcavAAR::kriglist. This function serves as an interface to kriging::kriging. We’ll get a list maps of data.frames with surface estimations for the two input layers.

lpoints <- list(df1, df2)

maps <- kriglist(lpoints, lags = 3, model = "spherical", pixels = 30)

The result of recexcavAAR::kriging is in a tall format – we have to transform it. For this purpose we use recexcavAAR::spatialwide.

surf1 <- spatialwide(maps[[1]]$x, maps[[1]]$y, maps[[1]]$pred, 3)
surf2 <- spatialwide(maps[[2]]$x, maps[[2]]$y, maps[[2]]$pred, 3)

After the transformation rgl can visualize the generated surfaces.

surface3d(
  surf1$x, surf1$y, t(surf1$z),
  color = c("black", "white"),
  alpha = 0.5,
  add = TRUE
)

surface3d(
  surf2$x, surf2$y, t(surf2$z),
  color = c("black", "white"),
  alpha = 0.5,
  add = TRUE
)

rglwidget()

During the excavation we created an artificial surface every 20 centimeters. Also we seperated the material of 1m*1m squares. Like this we get bodies of 1m*1m*0.2m. Let’s set up one by writing the corner coordinates for one into the data.frame hexatest:

hexatestdf <- data.frame(
  x = c(1, 1, 1, 1, 2, 2, 2, 2),
  y = c(0, 1, 0, 1, 0, 1, 0, 1),
  z = c(0.8, 0.8, 1, 1, 0.8, 0.8, 1, 1)
)

Now we can fill the shape with an equidistant three dimensional point raster using recexcavAAR::fillhexa and take a look at it. recexcavAAR::fillhexa can deal with completly amorphous hexahedrons.

cx = fillhexa(hexatestdf, 0.1)
completeraster <- points3d(
  cx$x, cx$y, cx$z,
  col = "red",
  add = TRUE
)

rglwidget()
# remove point raster from plot
rgl.pop(id = completeraster)

Damn! This spit penetrates both reconstructed surfaces. We should try to determine how his volume is distributed among the three major horizons of our trench. For this purpose we apply recexcavAAR::posdeclist (there’s also recexcavAAR::posdec to apply this to just one data.frame). It makes a position decision for every point of the artificial point raster we created with recexcavAAR::fillhexa.

cuberasterlist <- list(cx)

crlist <- posdeclist(cuberasterlist, maps)

hexa <- crlist[[1]]

a <- filter(
  hexa,
  pos == 0
)

b <- filter(
  hexa,
  pos == 1
)

c <- filter(
  hexa,
  pos == 2
)

points3d(
  a$x, a$y, a$z,
  col = "red",
  add = TRUE
)

points3d(
  b$x, b$y, b$z,
  col = "blue",
  add = TRUE
)

points3d(
  c$x, c$y, c$z,
  col = "green",
  add = TRUE
)

rglwidget()

Finally we can find out the percentual distribution. Could be an interesting information to determine the possible origin of finds from this spit.

sapply(
  crlist, 
  function(x){
    x$pos %>%
    table() %>%
    prop.table() %>%
    multiply_by(100) %>%
    round(2)
  }
  ) %>% t
##          0     1     2
## [1,] 24.49 64.24 11.27