---
title: "recexcavAAR: Vignette >>Trench visualisation<<"
author: "Clemens Schmid"
date: "August 2016"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Vignette >>Trench visualisation<<}
  %\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:

- **kriging:** simple surface modelling tool
- **rgl:** nice, interactive 3D plots

```{r, message=FALSE}
library(devtools)
library(recexcavAAR)
library(kriging)
library(rgl)
```

Like in the first vignette `Ifri el Baroud` we start by setting up an artificial and pretty simple excavation trench with a depth of 1.4 meters, a length of 6 meters and a width of 8 meters. This one is not parallel to the main axis of our coordinate system. 

```{r}
edges <- data.frame(
  x = c(6.899, 10.658, 4.428, 0.669, 6.899, 10.658, 4.428, 0.669),
  y = c(19.292, 14.616, 9.597, 14.273, 19.292, 14.616, 9.597, 14.273),
  z = c(9.7, 9.7, 9.7, 9.7, 8.3, 8.3, 8.3, 8.3)
)
```

We can plot the edges of this trench with `rgl::plot3d`, but first of all we have to calculate a reasonable aspectratio. Remember: The trench is tilted in relation to the main axis. 

```{r}
rangex <- abs(max(edges$x) - min(edges$x))
rangey <- abs(max(edges$y) - min(edges$y))

edgesordered = rbind(
  edges[1:4, ], 
  edges[1, ], 
  edges[5:8, ], 
  edges[5, ],
  edges[c(6,2), ],
  edges[c(3,7), ],
  edges[c(8,4), ]
)
```

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

```{r}
plot3d(
  edgesordered$x, edgesordered$y, edgesordered$z,
  type="l",
  aspect = c(rangex, rangey, 5),
  xlab = "x", ylab = "y", zlab = "z",
  sub = "Grab me and rotate me!",
  col = "darkgreen"
)

bbox3d(
  xat = c(2, 4, 6, 8, 10),
  yat = c(10, 12, 14, 16, 18),
  zat = c(8.5, 9, 9.5),
  back = "lines"
)
```

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

The approach for the excavation of this trench is to go deeper in squares of 1 meter by 1 meter. The spit depth is about 30 centimeters, but as we also try to follow natural layers the individual spit depth varies a lot. Fortunately we have niveau measurements of the resulting surfaces in the dataset `recexcavAAR::KT_spits`. 

```{r}
sp <- KT_spits

splist <- list()
spitnames <- c("^surface", "^spit1", "^spit2", "^spit3", "^bottom")

for (i in 1:length(spitnames)){
  splist[[i]] <- sp[grep(spitnames[i], sp$id), ]
}
```

Let's apply kriging with `recexcavAAR::kriglist`, transform the result with `recexcavAAR::spatialwide` and add the surfaces to the plot. 

```{r}
# I had to choose a very low pixel value to keep the vignette small enough 
maps <- kriglist(splist, x = 2, y = 3, z = 4, lags = 3, model = "spherical", pixels = 30)

surf <- list()
for (i in 1:length(maps)) {
  surf[[i]] <- spatialwide(maps[[i]]$x, maps[[i]]$y, maps[[i]]$pred, digits = 3)
}

idvec <- c()
for (i in 1:length(surf)) {
  idvec[i] <- surface3d(
    surf[[i]]$x, surf[[i]]$y, t(surf[[i]]$z),
    color = c("black", "white"),
    alpha = 0.5,
    add = TRUE
  )
}
```

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

```{r}
# remove surfaces from plot
for (i in 1:length(idvec)) {
  rgl.pop(id = idvec[i])
}
```

Hm... ok. But the surfaces extend beyond the surfaces of the trench... We can use `recexcavAAR::pnpmulti` to decide which points of the kriging result are within the trench polygon. Then we can remove the others and plot the result again.

```{r}
for (i in 1:length(maps)) {
  rem <- recexcavAAR::pnpmulti(edges$x[1:4], edges$y[1:4], maps[[i]]$x, maps[[i]]$y)
  maps[[i]] <- maps[[i]][rem, ]
}

surf2 <- list()
for (i in 1:length(maps)) {
  surf2[[i]] <- recexcavAAR::spatialwide(maps[[i]]$x, maps[[i]]$y, maps[[i]]$pred, 3)
}

for (i in 1:length(surf)) {
  surface3d(
    surf2[[i]]$x, surf2[[i]]$y, t(surf2[[i]]$z),
    color = c("black", "white"),
    alpha = 0.5,
    add = TRUE
  )
}
```

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

Now to the actual task: During the excavation we found a lot of pottery from all over the trench. Many of the individual sherds fit together -- it's possible to reconstruct complete vessels. We would like to visualize how the sherds of one particular vessel are distributed. The spatial information for this vessel is stored in the dataset `recexcavAAR::vessel`.

Some of the sherds were considered special during the excavation and were therefore measured individualy (Inventar Number *"KTF_..."*).

```{r}
ve <- KT_vessel
vesselsingle <- ve[grep("KTF", ve$inv), ]

points3d(
  vesselsingle$x, vesselsingle$y, vesselsingle$z,
  col = "red",
  size = 8, 
  add = TRUE
)
```

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

Unfortunately many were not. For these we just have the information from which spit and square they are coming from (Inv. Nr. *"KTM_..."*). 

```{r}
vesselmass <- ve[grep("KTM", ve$inv), ]
```

To visualize their position in the trench we could maybe set a point to the center of the respective square. The horizontal center of the square is easy to determine - but due to the irregular way of excavation which respects natural layers it's much more complex to get the vertical center of the square. 

Let's first of all load the dataset `recexcavAAR::KT_squarecorners` with the square point raster and create a list with the 4 individual corner points of each square. 

```{r}
sqc <- KT_squarecorners

squares <- list()
sqnum <- 1
for (i in 1:(nrow(sqc) - 9)) {
  if (i %% 9 == 0) {
    next
  } else {
    a <- sqc[i, ]
    b <- sqc[i + 1, ]
    c <- sqc[i + 9, ]
    d <- sqc[i + 10, ]
  }
  squares[[sqnum]] <- data.frame(
    x = c(a[, 1], b[, 1], c[, 1], d[, 1]), 
    y = c(a[, 2], b[, 2], c[, 2], d[, 2])
  )
  sqnum <- sqnum + 1
}
```

Now we can use `recexcavAAR::spitcenternatlist` to determine the list of spit center points in relation to the defined documentation surfaces.  

```{r}
sqcenters <- recexcavAAR::spitcenternatlist(squares, maps)

for (i in 1:length(sqcenters)) {
  sqcenters[[i]] <- data.frame(sqcenters[[i]], square = i, spit = c("spit1", "spit2", "spit3", "bottom"))
}

sqcdf <- do.call(rbind.data.frame, sqcenters)
```

Let's plot all of them for once. 

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

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

```{r}
# remove point raster from plot
rgl.pop(id = completeraster)
```

Finally we can merge the info about the vessel sherds in `vesselmass` with the new list of square center points `sqcdf`. Now we can plot single find sherds and mass find sherds together.

```{r warning=FALSE}
vmsq <- merge(vesselmass[, 1:4], sqcdf, by = c("square", "spit"), all.x = TRUE)

vesselm <- vmsq[complete.cases(vmsq), ]

points3d(
  vesselm$x, vesselm$y, vesselm$z,
  col = "orange",
  size = 8, 
  add = TRUE
)
```

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