library(magrittr)
library(knitr)
library(rgl)
library(ggplot2)
knit_hooks$set(webgl = hook_rgl)
view_matrix <- structure(c(0.586383819580078, 0.356217533349991, -0.727502763271332, 
0, -0.810031354427338, 0.257360488176346, -0.526888787746429, 
0, -0.000456457957625389, 0.898260772228241, 0.439460128545761, 
0, 0, 0, 0, 1), .Dim = c(4L, 4L))

bleiglas

bleiglas is an R package that provides some helper functions for 3D tessellation with voro++ and subsequent horizontal cutting of the resulting polygons for 2D plotting. The general workflow is described below.

Installation

You can install bleiglas from github

if(!require('devtools')) install.packages('devtools')
devtools::install_github("nevrome/bleiglas")

For the main function tessellate you also have to install the voro++ software. FOr Linux users: The package is already available in all major software repositories.

Get some data

I decided to use Dirk Seidenstickers Archives des datations radiocarbone d'Afrique centrale dataset for this purpose. It includes radiocarbon datings from Central Africa that combine spatial (x & y) and temporal (z) position with some meta information.

Click here for the data preparation steps

I selected dates from Cameroon between 1000 and 3000 uncalibrated BP and projected them into a worldwide cylindrical reference system (epsg 4088). As Cameroon is close to the equator this projection should represent distances, angles and areas sufficiently correct for this example exercise. I rescaled the temporal data with a factor of 1000 to better show the effect of 3D tessellation. You can imagine the samples to be observations in a 3D geo-time-space where one year equals one kilometre. I had to remove samples with equal position in all three dimensions for the tessellation.

# download raw data
c14_cmr <- c14bazAAR::get_c14data("adrac") %>% 
  # filter data
  dplyr::filter(!is.na(lat) & !is.na(lon), c14age > 1000, c14age < 3000, country == "CMR") 

# remove doubles
c14_cmr_unique <- c14_cmr %>%
  dplyr::mutate(
    rounded_coords_lat = round(lat, 3),
    rounded_coords_lon = round(lon, 3)
  ) %>%
  dplyr::group_by(rounded_coords_lat, rounded_coords_lon, c14age) %>%
  dplyr::filter(dplyr::row_number() == 1) %>%
  dplyr::ungroup()

# transform coordinates
coords <- data.frame(c14_cmr_unique$lon, c14_cmr_unique$lat) %>% 
  sf::st_as_sf(coords = c(1, 2), crs = 4326) %>% 
  sf::st_transform(crs = 4088) %>% 
  sf::st_coordinates()

# create active dataset
c14 <- c14_cmr_unique %>% 
  dplyr::transmute(
    id = 1:nrow(.),
    x = coords[,1], 
    y = coords[,2], 
    z = c14age * 1000, # rescaling of temporal data
    material = material
)

Data: c14

c14 

3D tessellation

Tessellation means filling space with polygons so that neither gaps nor overlaps occur. This is an exciting application for art (e.g. textile art or architecture) and an interesting challenge for mathematics. As a computational archaeologist I was already aware of one particular tessellation algorithm that has quiet some relevance for geostatistical analysis like e.g. spatial interpolation: Voronoi tilings that are created with Delaunay triangulation. These are tessellations where each polygon covers the space closest to one of a set of sample points.

Islamic mosaic with tile tessellations in Marrakech, Morocco. wiki
Delaunay triangulation and its Voronoi diagram. wiki
Output example of voro++ rendered with POV-Ray. math.lbl.gov

I learned that Voronoi tessellation can be calculated not just for 2D surfaces, but also for higher dimensions. The voro++ software library does exactly this for 3D space.

bleiglas::tessellate() is a minimal wrapper function that calls the voro++ command line interface (therefore you have to install voro++ to use it) for datasets like the one introduced above. We can apply it like this:

raw_voro_output <- bleiglas::tessellate(
  c14[, c("id", "x", "y", "z")],
  x_min = min(c14$x) - 150000, x_max = max(c14$x) + 150000, 
  y_min = min(c14$y) - 150000, y_max = max(c14$y) + 150000
)

bleiglas::tessellate(c14[, c("id", "x", "y", "z")]) would be sufficient, but I decided to increase the size of the tessellation box by 150 kilometres to each (spatial) direction to cover the area of Cameroon.

The output of voro++ is highly customizable, but structurally complex. With -v it first of all prints some config info on the command line.

Container geometry        : [937143:1.90688e+06] [63124.2:1.50658e+06] [1.01e+06:2.99e+06]
Computational grid size   : 3 by 5 by 6 (estimated from file)
Filename                  : /tmp/RtmpL9VlIm/file23304de43e55
Output string             : %i*%P*%t
Total imported particles  : 379 (4.2 per grid block)
Total V. cells computed   : 379
Total container volume    : 2.77155e+18
Total V. cell volume      : 2.77168e+18

It then produces an output file (*.vol) that can contain all sorts of geometry information for the calculated 3D polygons. I focussed on the edges of the polygons and wrote a parser function bleiglas::read_polygon_edges() that can transform the complex voro++ output to a tidy data.frame with six columns: the coordinates (x, y, z) of the start (a) and end point (b) of each polygon edge.

polygon_edges <- bleiglas::read_polygon_edges(raw_voro_output)

Data

polygon_edges

We can plot these polygon edges (black) together with the input sample points (red) in 3D.

Before plotting I wanted to changed the scaling of the temporal information back again to increase the readability of the plot.

polygon_edges %<>% dplyr::mutate(
  z.a = z.a / 1000,
  z.b = z.b / 1000
)

c14 %<>% dplyr::mutate(
  z = z / 1000
)
rgl::axes3d()
rgl::points3d(c14$x, c14$y, c14$z, color = "red")
rgl::aspect3d(1, 1, 1)
rgl::segments3d(
  x = as.vector(t(polygon_edges[,c(1,4)])),
  y = as.vector(t(polygon_edges[,c(2,5)])),
  z = as.vector(t(polygon_edges[,c(3,6)]))
)
rgl::view3d(userMatrix = view_matrix, zoom = 0.9)

rgl::axes3d()
rgl::points3d(c14$x, c14$y, c14$z, color = "red")
rgl::aspect3d(1, 1, 1)
rgl::segments3d(
  x = as.vector(t(polygon_edges[,c(1,4)])),
  y = as.vector(t(polygon_edges[,c(2,5)])),
  z = as.vector(t(polygon_edges[,c(3,6)]))
)
rgl::view3d(userMatrix = view_matrix, zoom = 0.9)

Cutting the polygons

The static 3D plot is of rather dubious value: Very hard to read. I therefore wrote bleiglas::cut_polygons() that can cut the 3D polygons at different levels of the z-axis. As the function assumes that x and y represent geographic coordinates, the cuts produce sets of spatial 2D polygons for different values of z -- in our example different points in time. The parameter cuts takes a numeric vector of cutting points and crs defines the spatial coordinate reference system of x and y to project the resulting 2D polygons correctly.

cut_surfaces <- bleiglas::cut_polygons(
  polygon_edges, 
  cuts = c(2500, 2000, 1500), 
  crs = 4088
)

Data

cut_surfaces

With this data we can plot a matrix of maps that show the cut surfaces.

cut_surfaces %>%
  ggplot() +
  geom_sf(
    aes(fill = time), 
    color = "white",
    lwd = 0.2
  ) +
  geom_sf_text(aes(label = id)) +
  facet_wrap(~time) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

cut_surfaces %>%
  ggplot() +
  geom_sf(
    aes(fill = time), 
    color = "white",
    lwd = 0.2
  ) +
  facet_wrap(~time) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

As all input dates come from Cameroon it might be a sensible decision to cut the polygon surfaces to the outline of this administrative unit.

cameroon_border <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>% 
  dplyr::filter(name == "Cameroon") %>% 
  sf::st_transform(4088)

cut_surfaces_cropped <- cut_surfaces %>% sf::st_intersection(cameroon_border)
cut_surfaces_cropped %>%
  ggplot() +
  geom_sf(
    aes(fill = time), 
    color = "white",
    lwd = 0.2
  ) +
  facet_wrap(~time) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

wzxhzdk:17
Finally we can also visualise any point-wise information in our input data as a feature of the tessellation polygons.

wzxhzdk:18 wzxhzdk:19

wzxhzdk:20

nevrome/bleiglas documentation built on Jan. 6, 2020, 12:45 a.m.