recexcavAAR: Vignette >>Trench visualisation<<

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.

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.

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

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

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.

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.

# 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
  )
}

# 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.

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

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_…”).

ve <- KT_vessel
vesselsingle <- ve[grep("KTF", ve$inv), ]

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

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

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.

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.

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)

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

# 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.

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
)

recexcavAAR: Vignette >>Trench visualisation<<

Clemens Schmid

August 2016