bleiglas is an R package that employs Voro++ for the calculation of three dimensional Voronoi diagrams from input point clouds. This is a special form of tessellation where each polygon is defined as the area closest to one particular seed point. Voronoi diagrams have useful applications in - among others - astronomy, material science or geography and bleiglas provides functions to make 3D tessellation more readily available as a mean for data visualisation and interpolation. It can be used for any 3D point cloud, but the output is optimized for spatiotemporal applications in archaeology.
- This README (see Quickstart guide below) describes a basic workflow with code and explains some of my thought process when writing this package.
- A JOSS paper gives some background, introduces the core functions from a more technical point of view and presents an example application.
- A (rather technical) vignette presents all the code necessary to
reproduce the “real world” example application in said JOSS paper.
When bleiglas is installed you can open the vignette in R with
vignette("bleiglas_case_study").
If you have questions beyond this documentation feel free to open an issue here on Github. Please also see our contributing guide.
You can install bleiglas from github
if(!require('remotes')) install.packages('remotes')
remotes::install_github("nevrome/bleiglas", build_vignettes = TRUE)For the main function tessellate you also have to install the Voro++
software. The package is already
available in all major Linux software repositories (on Debian/Ubuntu you
can simply run sudo apt-get install voro++.). MacOS users should be
able to install it via homebrew (brew install voro++).
For this quickstart, we assume you have packages tidyverse, sf,
rgeos (which in turn requires the Unix package geos) and c14bazAAR
installed.
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. As a minor pre-processing step, I here also remove samples with equal position in all three dimensions for the tessellation.
# download raw data with the data access package c14bazAAR
# c14bazAAR can be installed with
# install.packages("c14bazAAR", repos = c(ropensci = "https://ropensci.r-universe.dev"))
c14_cmr <- c14bazAAR::get_c14data("adrac") %>%
# filter data
dplyr::filter(!is.na(lat) & !is.na(lon), c14age > 1000, c14age < 3000, country == "CMR") ## | | | 0% | |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% | |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
# 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 = seq_len(nrow(.)),
x = coords[,1],
y = coords[,2],
z = c14age,
period = period
)Data: c14
c14 ## # A tibble: 393 × 5
## id x y z period
## <int> <dbl> <dbl> <int> <chr>
## 1 1 1284303. 450340. 1920 EIA
## 2 2 1101276. 321798. 2340 EIA
## 3 3 1101276. 321798. 2520 LSA
## 4 4 1093159. 264311. 2000 <NA>
## 5 5 1132077. 340034. 1670 <NA>
## 6 6 1101276. 321798. 2200 <NA>
## 7 7 1101276. 321798. 2030 <NA>
## 8 8 1101276. 321798. 1760 EIA
## 9 9 1093159. 264311. 1710 <NA>
## 10 10 1093159. 264311. 1940 <NA>
## # … with 383 more rows
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 quite some relevance for geostatistical analysis like 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 |
|---|---|---|
It turns out that Voronoi tessellation can be calculated not just for 2D surfaces, but also for higher dimensions. The Voro++ software library does exactly this for 3 dimensions. This makes it useful for spatio-temporal applications.
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,
unit_scaling = c(0.001, 0.001, 1)
)A critical step when using tessellation for spatio-temporal data is a
suitable conversion scale between time- and spatial units. Since 3D
tessellation crucially depends on the concept of a 3D-distance, we need
to make a decision how to combine length- and time-units. Here, for the
purpose of this example, we have 1 kilometre correspond to 1 year. Since
after the coordinate conversion our spatial units are given in meters,
we divide all spatial distances by a factor 1000 to achieve this
correspondence: unit_scaling = c(0.001, 0.001, 1).
I decided to increase the size of the tessellation box by 150 kilometres
to each (spatial) direction to cover the area of Cameroon. Mind that the
scaling factors in unit_scaling are also applied to the box size
parameters x_min, x_max, ….
The output of Voro++ is highly customizable, and structurally complex.
With the -v flag, the voro++ CLI interface prints some config info,
which is also the output of bleiglas::tesselate:
Container geometry : [937.154:1936.57] [63.1609:1506.58] [1010:2990]
Computational grid size : 3 by 5 by 6 (estimated from file)
Filename : /tmp/RtmpVZjBW3/file3aeb5f400f38
Output string : %i*%P*%t
Total imported particles : 392 (4.4 per grid block)
Total V. cells computed : 392
Total container volume : 2.8563e+09
Total V. cell volume : 2.8563e+09
It then produces an output file (*.vol) that contains all sorts of
geometry information for the calculated 3D polygons. tesselate returns
the content of this file as a character vector with the additionally
attached attribute unit_scaling
(attributes(raw_voro_output)$unit_scaling), which is just the scaling
vector we put in above.
I focussed on the edges of the polygons and wrote a parser function
bleiglas::read_polygon_edges() that can transform the complex Voro++
output for this specific output case to a tidy data.table with six
columns: the coordinates (x, y, z) of the start (a) and end point (b) of
each polygon edge. A data.table is a tabular R data structure very
similar to the standard data.frame. Read more about it
here.
polygon_edges <- bleiglas::read_polygon_edges(raw_voro_output)read_polygon_edges automatically reverses the rescaling introduced in
tesselate with the unit_scaling attribute.
Data: polygon\_edges
## x.a y.a z.a x.b y.b z.b polygon_id
## 1: 937154 374130 1307.99 1201480 392161 1299.80 25
## 2: 1289460 241706 1324.42 1201480 392161 1299.80 25
## 3: 1212280 387619 1290.18 1201480 392161 1299.80 25
## 4: 1190480 335990 1202.59 1233970 377206 1268.57 25
## 5: 1352310 233958 1240.81 1233970 377206 1268.57 25
## ---
## 24916: 1683200 887252 2655.00 1645270 892489 2655.00 290
## 24917: 1622180 900165 2682.50 1645270 892489 2655.00 290
## 24918: 1393030 1012170 2682.50 1622180 900165 2682.50 290
## 24919: 1596200 911750 2731.50 1622180 900165 2682.50 290
## 24920: 1645270 892489 2655.00 1622180 900165 2682.50 290
We can plot these polygon edges (black) together with the input sample points (red) in 3D.
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)
This 3D plot, even if rotatable using mouse input, is of rather limited
value since it’s 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 on the z axis.
bleiglas::cut_polygons() yields a rather raw format for specifying
polygons. Another function, bleiglas::cut_polygons_to_sf(), transforms
it to sf. Here 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)
) %>%
bleiglas::cut_polygons_to_sf(crs = 4088)Data: cut\_surfaces
## Simple feature collection with 76 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 937154 ymin: 63160.9 xmax: 1936570 ymax: 1506580
## Projected CRS: World Equidistant Cylindrical (Sphere)
## First 10 features:
## x z id
## 1 POLYGON ((1195386 319810.5,... 2500 3
## 2 POLYGON ((1936570 809055.4,... 2500 31
## 3 POLYGON ((1146675 374628.2,... 2500 38
## 4 POLYGON ((1215947 365177.1,... 2500 40
## 5 POLYGON ((1416056 455852, 1... 2500 69
## 6 POLYGON ((1083149 968036.7,... 2500 103
## 7 POLYGON ((1936570 315020.3,... 2500 105
## 8 POLYGON ((1386575 333838.1,... 2500 135
## 9 POLYGON ((1116416 63160.9, ... 2500 144
## 10 POLYGON ((1377347 63160.9, ... 2500 185
With this data we can plot a matrix of maps that show the cut surfaces.
cut_surfaces %>%
ggplot() +
geom_sf(
aes(fill = z),
color = "white",
lwd = 0.2
) +
geom_sf_text(aes(label = id)) +
facet_wrap(~z) +
theme(
axis.text = element_blank(),
axis.ticks = element_blank()
)
As all input dates come from Cameroon it makes sense 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 = z),
color = "white",
lwd = 0.2
) +
facet_wrap(~z) +
theme(
axis.text = element_blank(),
axis.ticks = element_blank()
)
Finally, we can also visualise any point-wise information in our input data as a feature of the tessellation polygons.
cut_surfaces_material <- cut_surfaces_cropped %>%
dplyr::left_join(
c14, by = "id"
)cut_surfaces_material %>%
ggplot() +
geom_sf(
aes(fill = period),
color = "white",
lwd = 0.2
) +
facet_wrap(~z.x) +
theme(
axis.text = element_blank(),
axis.ticks = element_blank()
)
This quickstart was a simple primer on how to use this package. If you think the final use case wasn’t too impressive, take a look at this analysis of Bronze Age burial types through time, as performed in our JOSS paper and the vignette.
To cite bleiglas in publications use:
Schmid and Schiffels (2021). bleiglas: An R package for interpolation
and visualisation of spatiotemporal data with 3D tessellation.
Journal of Open Source Software, 6(60), 3092,
https://doi.org/10.21105/joss.03092
A BibTeX entry for LaTeX users is
@Article{,
title = {{bleiglas}: An {R} package for interpolation and visualisation of spatiotemporal data with 3D tessellation},
author = {Clemens Schmid and Stephan Schiffels},
journal = {Journal of Open Source Software},
volume = {6},
number = {60},
pages = {3092},
year = {2021},
doi = {10.21105/joss.03092},
url = {https://doi.org/10.21105/joss.03092},
}