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(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`

:

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

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

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
```

```
## 0 1 2
## [1,] 22.54 62.21 15.25
```