# recexcavAAR: Vignette >>Ifri el Baroud<<

#### July 2016

First we need to load recexcavAAR and some external packages:

• dplyr: filter function
• kriging: simple surface modelling tool
• magrittr: introduces piping via %>% operator
• rgl: nice, interactive 3D plots
library(devtools)
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",
)

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,
)

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",
)

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",
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 %>%
) %>% t
##          0     1     2
## [1,] 24.49 64.24 11.27