recexcavAAR: Vignette >>Semiautomatic spit attribution<<

library(devtools)
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)
)

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),
  back = "lines"
)

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
)

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)

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
)

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
)

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
)

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

recexcavAAR: Vignette >>Semiautomatic spit attribution<<

Clemens Schmid

July 2016