---
title: "recexcavAAR: Vignette >>Semiautomatic spit attribution<<"
author: "Clemens Schmid"
date: "July 2016"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Vignette >>Semiautomatic spit attribution<<}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, echo=FALSE}
# check if pandoc is available
if (requireNamespace("rmarkdown") && !rmarkdown::pandoc_available("1.13.1"))
stop("These vignettes assume pandoc version 1.13.1; older versions will not work.")
# see https://r-forge.r-project.org/forum/message.php?msg_id=43797&group_id=234
```

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

```{r, message=FALSE}
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. 

```{r}
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`:

```{r, echo=FALSE, results="hide"}
# avoid plotting in X11 window
open3d(useNULL = TRUE)
```

```{r}
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"
)
```

```{r, echo=FALSE, fig.width=7, fig.height=5}
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.  

```{r}
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`.

```{r}
points3d(
  df1$x, df1$y, df1$z,
  col = "darkgreen",
  add = TRUE
)

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

```{r, echo=FALSE, fig.width=7, fig.height=5}
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. 

```{r }
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`.

```{r}
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.

```{r}
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
)
```

```{r, echo=FALSE, fig.width=7, fig.height=5}
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`:

```{r fig.width=7, fig.height=5}
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.

```{r}
cx = fillhexa(hexatestdf, 0.1)
```

```{r}
completeraster <- points3d(
  cx$x, cx$y, cx$z,
  col = "red",
  add = TRUE
)

```

```{r, echo=FALSE, fig.width=7, fig.height=5}
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`.

```{r}
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
)
```

```{r, echo=FALSE, fig.width=7, fig.height=5}
rglwidget()
```

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

```{r}
sapply(
  crlist, 
  function(x){
    x$pos %>%
    table() %>%
    prop.table() %>%
    multiply_by(100) %>%
    round(2)
  }
  ) %>% t
```