vector_cubes.Rmd

---
title: "Visualization of vector data cubes with tmap"
author: "Martijn Tennekes"
date: "9/18/2020"
output: rmarkdown::github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Installation of stars and tmap 

We recommend to install the github version of **stars**, and the CRAN version of **tmap**:

```{r, eval=FALSE}
if (!require(remotes))
  install.packages("remotes")
remotes::install_github("r-spatial/stars")
install.packages("starsdata", repos = "https://gis-bigdata.uni-muenster.de/pebesma", type = "source")
install.packages("tmap")
```

Load the packages:

```{r}
library(stars)
library(tmap)
```

## North Carolina example

This example is the North Carolina SIDS (sudden infant death syndrome) data set. See https://edzer.github.io/UseR2020/#Vector_data_cube:_North_Carolina_SIDS

The code to create the vector cube is the following:

```{r}
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))
nc.df = st_set_geometry(nc, NULL) # m is a regular, non-spatial data.frame
head(nc.df)

mat = as.matrix(nc.df[c("BIR74", "SID74", "NWBIR74", "BIR79", "SID79", "NWBIR79")])
dim(mat) = c(county = 100, var = 3, year = 2) # make it a 3-dimensional array
# set dimension values to the array:
dimnames(mat) = list(county = nc$NAME, var = c("BIR", "SID", "NWBIR"), year = c(1974, 1979))
# convert array into a stars object
nc.st = st_as_stars(pop = mat)

# set dimension values to polygons
nc.geom <- st_set_dimensions(nc.st, 1, st_geometry(nc))
nc.sum = sapply(split(nc.geom, 2), sum)

# calculate standardise incidence rates
IR = nc.sum[2]/nc.sum[1]
(nc.SIR = st_apply(nc.geom, c(1,3), function(x) (x[2]/x[1])/IR))
```

## Visualization with tmap

The next version of **tmap** will support those vector cubes natively. The code to create two maps, one for each year, should be as simple as:

```{r,eval=FALSE}
# this code will not work yet, but in the next version of tmap (3.3)
tm_shape(nc.SIR) +
  tm_polygons("pop") +
  tm_facets(by = "year")
```

A workaround using the current version of tmap (3.2) that produces this map is:

```{r}
# extract the years
years = st_get_dimension_values(nc.SIR, 2)

# convert the stars object to an sf object
nc.SIR.sf = st_as_sf(nc.SIR)

# plot it as follows
tm_shape(nc.SIR.sf) + tm_polygons(years, title = "SIR", palette = "Purples", n = 10) + 
  tm_facets(free.scales = FALSE) +
  tm_layout(panel.labels = years)
```