vdc-showcase.qmd

---
title: "Computational notebook"
bibliography: https://api.citedrive.com/bib/618133a2-4856-443b-b778-a72d316eb599/references.bib?x=eyJpZCI6ICI2MTgxMzNhMi00ODU2LTQ0M2ItYjc3OC1hNzJkMzE2ZWI1OTkiLCAidXNlciI6ICIyNTY4IiwgInNpZ25hdHVyZSI6ICI4MWExNDU3ODE3NWUxMzFiZDI1NDliOTg2NmQzMWJjNzU4OWI4YjI3MDUxOTU5NDg3MDM5MzFlYmQwMDY0OThhIn0=/bibliography.bib
format: 
  gfm: 
    toc: true
    toc-depth: 3
    number-sections: true
    fig-dpi: 300
  html:
    toc: true
    toc-depth: 3
    number-sections: true
    fig-dpi: 300
---

# Context

This is a computational notebook to accompany the short paper manuscript "Vector data cubes for features evolving in space and time" submitted to AGILE 2024. 

In this notebook, the expanded code on how to create vector data cubes (VDC) for shape-shifting features as described in section 3.2 (Vector data cube implementation) is included.

Note that the stand-alone code to recreate the figures of the manuscript can be browsed here ([manuscript-figures.qmd](code/manuscript-figures.qmd)).

These experiments are performed in R v. 4.3.2 and further details on library versions are included at the end of the notebook.

# Libraries

The following chunk installs the needed libraries. Run only if overwriting current installations is not an issue. 
Note that `cubble` and `ggnewscale` come from the GitHub repository.
`tmap` is installed from r-universe to work with version 4.

```{r}
#| eval: false
#| message: false
#| warning: false
pkgs = c(
  "dplyr","ggplot2","here","patchwork",
  "purrr","readr","stringr",
  "tidyr","tsibble","units","zen4R"
)
install.packages(pkgs)
remotes::install_github("huizezhang-sherry/cubble")
remotes::install_github("eliocamp/ggnewscale")
remotes::install_github("loreabad6/sf")
remotes::install_github("r-spatial/stars")
install.packages(
  'tmap',
  repos = c('https://r-tmap.r-universe.dev',
            'https://cloud.r-project.org')
)
```

The libraries needed to reproduce this notebook are listed below.

```{r}
#| message: false
#| warning: false
library(cubble) # tabular VDCs
library(dplyr) # data wrangling
library(ggplot2) # visualisation
library(ggnewscale) # add second color/fill scale in ggplot
library(here) # manage current directory
library(patchwork) # plot composer
library(purrr) # functional programming 
library(readr) # read CSV files
library(sf) # read vector files
library(stars) # array VDCs
library(stringr) # handle strings
library(tidyr) # create tidy data
library(tmap) # spatial visualisation
library(tsibble) # handle time series data
library(units) # set units
# library(zen4R) # download data from zenodo
```

# Data and pre-processing

Both data sets are shared through Zenodo, therefore we will directly fetch the necessary data from this repository using `zen4R`.

## Lava flows in Fagradalsfjall, Iceland

The lava flow outlines for the Fagradalsfjall eruption can be obtained from @pedersen_2023_7866738. 

### Repository description

This is the description attached to the repository.

```{r}
#| eval: false
#| results: "asis"
zenodo = get_zenodo(
  doi = "10.5281/zenodo.7866738"
)
zenodo$metadata$description
```

### Files

As we see, the vector data format is GeoPackage.
These are the files contained in the repository:

```{r}
#| eval: false 
zenodo$files
```

From an inspection on the Zenodo listing we could see that the outline for each date is saved as a different file.

### Download and read

Hence, we will download only the outlines zipped file, unzip and create a list of the GeoPackage files to read into memory.

```{r}
#| label: downloadlf
#| eval: false
dir = tempdir()
download_zenodo(
  doi = "10.5281/zenodo.7866738",
  path = dir,
  files = list(
    "outlines_pedersen_etal2022_v12.zip"
  ),
  overwrite = FALSE,
  timeout = 600
)

# Unzip
files = list.files(here(dir), full.names = TRUE)
lapply(files, unzip, exdir = here(dir, "unzipped"))

# Find geopackage files
fn_gpkg = list.files(
  here(dir, "unzipped"),
  pattern = "Outline.*gpkg$",
  full.names = TRUE,
  recursive = TRUE
)
```

```{r}
#| label: temporarychunklfdata
#| eval: true
#| echo: false
data_dir = here("data", "raw", "lava_flows",
                "fagradalsfjall")
fn_gpkg = list.files(
  here(data_dir, "unzipped"),
  pattern = "Outline.*gpkg$",
  full.names = TRUE,
  recursive = TRUE
)
```

With this we can read in the files into R as an `sf` object, where we combine all the files into one. 
The filename is kept to extract the date information from it.

```{r}
#| label: readlf
# Create function to read in each file and to extract 
# the date from the filename
read_fun = function(x){
  read_sf(x) |> 
    mutate(
      fn_gpkg = tail(str_split(x, '/')[[1]], n=1),
      datetime = as.POSIXct(
        paste(str_split(fn_gpkg, "_")[[1]][2:3], collapse = ""),
        format = "%Y%m%d%H%M"
      )
    ) |> 
    # One file does not have a crs set
    st_set_crs(3057)
}

# Call map2 to read in the files
# and bind them into one single sf object
outlines_lf = map(fn_gpkg, read_fun) |> bind_rows()
```

### Explore

We can look at the final object and have an initial idea of the contents. 

```{r}
outlines_lf |> glimpse()
```

As seen, the columns accompanying the data are not complete and also do not provide extra information. 
We will mainly focus on the geometry and the datetime registered.
We can take a short look at the polygons included in the data:

```{r}
ggplot(outlines_lf) +
  geom_sf(aes(color = datetime), fill = NA) 
```

### Prepare for VDC

The final step is to make sure we have a single geometry per timestamp.
For this, we group by datetime and summarise the geometries for each data with a union.
We add a step to make the geometry valid in case the union operation caused any topological issues.

```{r}
outlines_lf = outlines_lf |> 
  # Combine polygons from single date into multipolygons
  group_by(datetime) |> 
  summarise(geom = st_combine(geom)) |> 
  ungroup() |> 
  # Make geometry valid
  st_make_valid()
```

Now our data is ready to be coerced to a VDC.

## Landslide in Butangbunasi, Taiwan

The Butangbunasi landslide outlines can be obtained from @Hoelbling_2024_Butangbunasidata.

### Repository description

This is the description attached to the repository.

```{r}
#| eval: false
#| results: "asis"
zenodo = get_zenodo(
  doi = "10.5281/zenodo.10635102"
)
zenodo$metadata$description
```

### Files

Again, the vector data format is GeoPackage.
These are the files contained in the repository:

```{r}
#| eval: false
zenodo$files
```

### Download and read

From exploring the repository and checking the description, we concluded the "Butangbunasi_OBIA_statistics.csv" also contains valuable information that can be join into the outline data frames.

```{r}
#| label: downloadldsl
#| eval: false
# Fetch also the CSV file with additional info
dir = tempdir()
download_zenodo(
  doi = "10.5281/zenodo.10635102",
  path = dir,
  files = list(
    "outlines.zip",
    "Butangbunasi_OBIA_statistics.csv"
  ),
  overwrite = FALSE,
  timeout = 100
)

# Unzip
files = list.files(here(dir), full.names = TRUE)
lapply(files, unzip, exdir = here(dir, "unzipped"))

# Find geopackage files
mapping_ls = list.files(
  here(dir, "unzipped", "outlines"), 
  pattern = ".gpkg",
  full.names = TRUE
)
```

```{r}
#| label: temporarychunkldsldata
#| eval: true
#| echo: false
dir = here("data", "raw", "landslides",
                "butangbunasi")
mapping_ls = list.files(
  here(dir, "unzipped", "outlines"), 
  pattern = ".gpkg",
  full.names = TRUE
)
```

We can then read in the CSV file and the outline files.
With a left join, we can obtain the date and other columns from the CSV file. 

```{r}
#| label: readldsl
#| warning: false
#| message: false
# Read in CSV file
stats = read_csv(
  here(dir, "Butangbunasi_OBIA_statistics.csv")
)

# Set-up a read function that fetches filename
read_fun = function(x){
  read_sf(x) |> 
    mutate(
      fn_gpkg = tail(str_split(x, '/')[[1]], n=1),
    )
}

# Read files and combine into single sf object
# Join with CSV file
outlines_ldsl = lapply(mapping_ls, read_fun) |> 
  bind_rows() |> 
  left_join(stats, by = c("fn_gpkg" = "file"))
```

### Explore

We can look at the final object and have an initial idea of the contents. 

```{r}
outlines_ldsl |> glimpse()
```

For this case, the CSV file added some extra information to the data, such as the satellite sensor and the landslide area.
We see also the column class and lake area which suggests there is more than one class for this dataset. 
Let's look at the polygons included in the data:

```{r}
ggplot(arrange(outlines_ldsl, desc(date))) +
  geom_sf(aes(color = date, fill = Class), alpha = 0.2) +
  scale_color_viridis_c(trans = "date") +
  scale_fill_manual(values = c("deepskyblue3", "peru"))
```

We see here that the data does not only contain the landslide outlines but also the delineation of a dammed lake. 

### Prepare for VDC

```{r}
outlines_ldsl = outlines_ldsl |> 
  # Coerce date column to Date class
  mutate(date = as.Date(date)) |> 
  # Remove outline for 2018
  # this is a reference outline and does not follow a
  # Typhoon event according to Hoelbling et al., 2020
  filter(date != as.Date("2018-11-08")) |>
  # Group by date which in combination with summarise
  # will union all single polygons into a multipolygon
  group_by(date, Class) |> 
  # Summarise sensor and area information for landslide
  summarise(
    sensor = first(satellite_sensor),
    area_ldsl = set_units(first(landslide_area_ha), 'ha'),
    area_lake = set_units(first(lake_area_ha), 'ha')
  ) |> 
  ungroup() |> 
  # collapse the area column into a single one depending on Class
  mutate(
    area = case_when(
      Class == "landslide" ~ area_ldsl,
      Class == "lake" ~ area_lake
    )
  ) |> 
  select(-starts_with("area_")) |> 
  # Complete data for missing cases of date-Class combination
  complete(date, Class) |> 
  # rename the class column (lowercase)
  rename(class = Class) |> 
  # convert back to sf
  st_as_sf()
```

# VDC creation

## Array format

We will start creating an array VDC for the shape-evolving features.
For this, we create an array object containing the data, i.e., the measure, that in this case are the changing geometries, and the dimensions of the array that correspond to the length of unique summary geometries (in this case 1), and of unique datetime strings.

Finally, we create a dimensions object that includes the *summary geometry* for our shape-evolving features.
This can be the union of all the changing geometries or the centroid of those or the point from which the lava erupts or the landslide scar, or any other point the analyst sees adequate.

### Lava flow VDC

For the lava flow, we do not use a unique identifier (`id`) because this is the only feature of this class under analysis.
That means we are not including a second lava flow, e.g. the one in Grindavík, in the same cube, although we could when data is available.

```{r}
# Create array
a = array(
  data = outlines_lf$geom, 
  dim = c(
    1,
    length(unique(outlines_lf$datetime))
  ),
  dimnames = list(
    geom_sum = 1,
    datetime = unique(outlines_lf$datetime)
  )
)

# Create dimensions object
# Summary geometry is the centroid of the union of all geometries
# The point parameter indicates if the value refers to a point (location)
# or to a pixel (area) value
dim_cent = st_dimensions(
  geom_sum = st_centroid(st_make_valid(st_union(outlines_lf$geom))), 
  datetime = unique(outlines_lf$datetime),
  point = c(TRUE, FALSE)
)

# Coerce to cube
(cube_arr_lf = st_as_stars(
  list(geometry = a), 
  dimensions = dim_cent)
)
```

If we would use the union we will see under the values column of the `geom_sum` dimension not a `POINT` geometry as above but a `POLYGON` geometry.
We illustrate this with an unsaved object, so we will keep working with the centroid as the `geom_sum`.

```{r}
dim_union = st_dimensions(
  geom_sum = st_make_valid(st_union(outlines_lf$geom)), 
  datetime = unique(outlines_lf$datetime),
  point = c(TRUE, FALSE)
)

# Coerce to cube
st_as_stars(
  list(geometry = a), 
  dimensions = dim_union)
```

### Landslide VDC

For the landslide example we have added also the landslide-dammed lake outlines. 
The lake has been present in seven of the Landsat scenes mapped and therefore can serve as an example on how to handle absence in the data, as well as a second feature set in a VDC.

In this example, we add `class` as an extra dimension, which will mark whether the shape-evolving geometry is a landslide or a lake.
For this particular case this dimension becomes redundant since there are two classes and two summary geometries, however, when we deal with more landslides and associated lakes, this could become useful.
We will notice data redundancy by checking the amount of geometries in the summary of the attributes, we had initially 40 geometries and now we have 80, with 53 empty geometries.

We will perform some extra data wrangling to create this cube. 
Namely, we compute the *summary geometry* on the original `sf` object and perform again a complete case check to fill values with `NA`s for every combination of `date`, `class` and `geom_sum`.

Further, this dataset contains other information that can be mapped to the attributes such as the sensor and area.
We will construct separate arrays for these and add them all to the cube.

```{r}
# Add geom_sum columns to original sf data
# will be used also to create tabular VDC
outlines_ldsl_m = outlines_ldsl |> 
  # Summary geometry is the centroid of the union 
  # of all geometries per type
  group_by(class) |> 
  mutate(
    geom_sum = st_centroid(st_make_valid(st_union(geom)))
  ) |> 
  relocate(geom_sum, .after = everything()) |> 
  ungroup() |> 
  st_as_sf(sf_column_name = 'geom_sum') 

# Complete cases, saved as separate object, see tabular format
outlines_ldsl_c = outlines_ldsl_m |> 
  complete(date, class, geom_sum) |> 
  st_as_sf(sf_column_name = 'geom_sum')

# Create arrays
dims = c(
    length(unique(outlines_ldsl_c$geom_sum)),
    length(unique(outlines_ldsl_c$class)),
    length(unique(outlines_ldsl_c$date))
  )
dnms = list(
    geom_sum = unique(outlines_ldsl_c$geom_sum),
    class = unique(outlines_ldsl_c$class),
    datetime = unique(outlines_ldsl_c$date)
  )
geom = array(
  data = outlines_ldsl_c$geom, 
  dim = dims,
  dimnames = dnms
)
area = array(
  data = outlines_ldsl_c$area, 
  dim = dims,
  dimnames = dnms
)
sensor = array(
  data = outlines_ldsl_c$sensor, 
  dim = dims,
  dimnames = dnms
)

# Summarise geometries to add in dimensions object
sumgeom = outlines_ldsl_c |> 
  group_by(class) |> 
  summarise(geom_sum = first(geom_sum))

# Create dimensions object
# The point parameter indicates if the value refers to a point (location)
# or to a pixel (area) value
dim_cent = st_dimensions(
  geom_sum = sumgeom$geom_sum, 
  class = unique(outlines_ldsl_c$class),
  date = unique(outlines_ldsl_c$date)
  # point = c(TRUE, TRUE, FALSE)
)

# Coerce to cube
(cube_arr_ldsl = st_as_stars(
  list(geometry = geom, area = set_units(area, "ha"), sensor = sensor), 
  dimensions = dim_cent)
)
```

## Tabular format

For the tabular format we use `cubble`. 
We can either coerce the `stars` VDC into a `cubble` (as shown for the lava flow example), or we can create a `cubble` from the original `sf` object (as shown for the landslide example).

For this format, we need to define a `key` and `index.` The key is the identifier of the spatial face of the cube while the index is the identifier of the temporal face.
We then use `index = time` and we should use `key = geom_sum`. 
This works when creating the `cubble` object but the algorithm to display the temporal face gets stuck as it tries to add a WKT geometry to every entry. 
There is one "hack" to handle this when coercing from a `stars` object.
Internally, `cubble` adds an `id` column to the data, so we can call `key = id` and the cube creation in its spatial and temporal forms will work.
This is the approach we use for the lava flow VDC.
For the landslide VDC we added the column class as the key, which becomes then the identifier.
Better ways to handle this will be suggested to `cubble` in future work. 

### Lava flow VDC

As mentioned above, this VDC is created from the `stars` VDC. 

```{r}
cube_tab_lf = as_cubble(
  cube_arr_lf, key = id, index = datetime
)
```

We can see the temporal and spatial faces of the cube.
As we will see `ts` in the spatial table is a list column where the time series is stored as displayed in the temporal table.

```{r}
cube_tab_lf |> 
  face_spatial()
cube_tab_lf |> 
  face_temporal()
```

What is interesting from this tabular format is that the `tsibble` package is used under the hood to handle the temporal aspect.
Therefore, information on time intervals is inferred from the given timestamps and it is therefore determined if the time series has any gaps.

### Landslide VDC

For this VDC we will use the modified `sf` object for the landslide dataset.
We have included here the version without the complete cases for the three dimensions because `cubble` throws an error for this (there is a duplication of indices).
Further investigation into why this is the case will follow for future work.

```{r}
## Create cubble
cube_tab_ldsl = as_cubble(
  outlines_ldsl_m, key = class, index = date
)
```

The `cubble` creation is now based on the object with complete cases for every combination of `date` and `class`.
This means `NA` values and empty geometries are included.
However, we could also build the `cubble` without this `NA` values, and hence avoid data redundancy (note the dimensions of the tibble in the list-column `ts`.

```{r}
as_cubble(
  drop_na(outlines_ldsl_m), key = class, index = date
)
```

For compatibility with the array format, we will keep the first approach as the tabular VDC.
We can see the temporal and spatial faces of the cube.
Here we note how the "lake" and "landslide" become the identifiers for the feature sets.

```{r}
cube_tab_ldsl |> 
  face_spatial()
cube_tab_ldsl |> 
  face_temporal()
```

# Showcase

In this section, we will demonstrate some of the operations that are supported now for VDCs, and some of the operations that we would like to implement in future work.

We will use the created cubes (lava flow, landslide) interchangeably, considering that most of the operations (given the data structure, e.g., extra dimensions, extra attributes) can be applied to both examples.

## Cube operations

Most of `dplyr` verbs and `sf` functions work for the array and tabular cubes.
For the tabular cubes, it is important to remember which face we want to apply the operation to.

### Geometric measurements

```{r}
cube_arr_lf |> 
  mutate(area = st_area(geometry))

cube_tab_lf |> 
  face_temporal() |> 
  mutate(area = st_area(geometry))
```

### Subsetting dimensions

Subsetting refers to lookup operations via an index.
A way to subset a cube is to use the slice verb. 
Similar results can be achieved with `[`, but the syntax is not so intuitive and we need to know the dimension order.
<!-- When a dimension after slicing has a single index it is dropped by default. -->
<!-- Note how the number of geometries in the attributes changes. -->

```{r}
cube_arr_lf |> 
  slice(index = 3:12, along = "datetime")

cube_arr_lf[,,3:12]
```

```{r}
#| echo: false
#| eval: false
cube_arr_ldsl |> 
  slice(index = 1, along = "class")
```

### Selecting attributes

We can select single attributes or columns. 
Note that for the tabular format the `key` and `index` behave as sticky columns and always stay.

```{r}
cube_arr_ldsl |> 
  select(sensor)

cube_tab_ldsl |> 
  face_temporal() |> 
  select(sensor)
```

### Filtering time

```{r}
cube_arr_lf |> 
  filter(datetime > "2021-03-18", datetime < "2021-03-25")

cube_tab_lf |> 
  face_temporal() |> 
  filter(datetime > "2021-03-18", datetime < "2021-03-25")
```

### CRS transformation

Transforming the CRS of the VDC works for both the tabular and array formats.

For the tabular format, the changing geometries preserve their original CRS, even though the *summary_geometry* is modified for both formats.

For the array format, this issue is fixed thanks to [@edzer](https://github.com/edzer), see [#1](https://github.com/loreabad6/vdc-space-time-feats/issues/1).

In future implementations, the tabular format will have a better integration of the *summary geometry* and *shape-evolving geometries*, where both geometry-columns get updated with such operations.

```{r}
test_arr = cube_arr_lf |> 
  st_transform(4326)
test_arr
test_arr$geometry

test_tab = cube_tab_lf |>  
  st_transform(4326)
test_tab
test_tab |> 
  face_temporal() 
```

### Differences in time series

For the tabular format, we can use `tsibble` functions when coercing our cube into a temporal tsibble.
One useful function is to calculate differences between timestamps.
We use this for computing the difference in area (per feature class for the landslide example) below.

```{r}
cube_tab_lf |> 
  face_temporal() |> 
  mutate(area = set_units(st_area(geometry), "ha")) |> 
  make_temporal_tsibble() |> 
  mutate(diff = difference(area, order_by = datetime))

cube_tab_ldsl |> 
  face_temporal() |> 
  make_temporal_tsibble() |> 
  group_by(class) |> 
  mutate(diff = difference(area, order_by = date)) |> 
  ungroup() |> 
  arrange(rev(class))
```

Functions like spatial filtering and data aggregations would become more relevant when we deal with a larger number of feature sets, and will be explored for future implementations of the VDC concepts.

## Visualisation

### Maps with array VDCs

`stars` plot methods in base R are a fast way to inspect information in data cubes.
For the VDC array implementation, however, the `plot.stars()` method does not work by default,
as the inclusion of the geometry as an attribute is not supported.

```{r}
#| error: true
plot(cube_arr_lf)
```

A workaround is to call the attribute directly, while if we want to perform operations in a pipeline and then plot we would need to use a lambda function.

```{r}
oldpar = par(no.readonly = TRUE)
par(mar = c(1,1,1,1), bg = "white")
plot(cube_arr_lf$geometry)
```

```{r}
par(mar = c(1,1,1,1), bg = "white")
cube_arr_lf |>
  filter(datetime > "2021-03-18", datetime < "2021-03-25") |> 
  (\(x) plot(x$geometry))()
```


Tidyverse plotting methods, i.e., `ggplot`s implementation for stars `geom_stars()` also does not know how to handle this geometry attribute.

```{r}
#| error: true
ggplot(cube_arr_lf) +
  geom_stars()
```

Here, going to the tabular approach is a better solution. 
A pre-step is to assign the geometry to the changing geometry instead of `geom_sum`.

```{r}
cube_tab_lf |> 
  face_temporal() |> 
  # rearrange df for plot order
  arrange(desc(datetime)) |> 
  # assign sf column to changing geometry
  (\(x) st_as_sf(x, sf_column_name = "geometry"))() |> 
  ggplot() +
  geom_sf()
```

The advantage of using `ggplot2` is for example the use of facets.

```{r}
cube_tab_lf |> 
  face_temporal() |> 
  # rearrange df for plot order
  arrange(desc(datetime)) |> 
  # assign sf column to changing geometry
  (\(x) st_as_sf(x, sf_column_name = "geometry"))() |> 
  ggplot() +
  geom_sf() +
  facet_wrap(~datetime) +
  theme_void()
```

Other plotting packages like tmap (v4) can also be used in this manner.

```{r}
#| fig-height: 5
#| fig-width: 8
cube_tab_ldsl |> 
  face_temporal() |> 
  # rearrange df for plot order
  arrange(desc(date)) |> 
  # assign sf column to changing geometry
  (\(x) st_as_sf(x, sf_column_name = "geom"))() |> 
  tm_shape() +
  tm_polygons(
    fill = "date", fill_alpha = 0.9, col = "white", lwd = 0.3,
    fill.scale = tm_scale(values = "-viridis")
  ) +
  tm_facets(columns = "class")
```

### Multi-dimensional plot

In the paper we included an example of a multi-dimensional visualisation of the VDC data.
This was achieved by adding a distortion to the geometries to make them appear stacked over each other.

```{r}
#| fig-height: 10
#| fig-width: 7
# Create shear matrix
sm = matrix(c(2.5, 1.2, 0, 1), 2, 2)

# Apply shear matrix
ldsl_shear = outlines_ldsl |>
  mutate(
    geom = geom * sm,
    # sequence along date, i.e. 1 to 20 per class
    shift_y = rep(
      1:length(unique(date)),
      each = length(unique(class))
    )
  ) |> 
  # Add lost crs
  st_set_crs(st_crs(outlines_ldsl)) 

ldsl_shift = ldsl_shear |> 
  rowwise() |> 
  # add a shift to stack outlines on top of each other
  mutate(
    geom = geom + c(0, shift_y * 4000),
    # Marker to place date labels
    y_label = st_coordinates(st_centroid(geom))[,'Y']
  ) |> 
  ungroup() |> 
  st_as_sf()

ggplot(ldsl_shift) +
  # plot landslide geometries
  geom_sf(
    data = filter(ldsl_shift, class == "landslide"),
    aes(fill = date), 
    color = "black",
    show.legend = FALSE
  ) +
  # plot lake geometries
  geom_sf(
    data = filter(ldsl_shift, class == "lake"),
    color = "red", fill = "pink",
    show.legend = FALSE
  ) +
  # Add date label, the value of x is added after visual inspection
  geom_text(
    data = filter(ldsl_shift, class == "landslide"),
    aes(label = date, y = y_label),
    x = 3751500, size = 3.5
  ) +
  scale_fill_viridis_c("Date", direction = 1, trans = "date", option = "D") +
  # Limits are expanded after visual inspection
  coord_sf(xlim = c(3745000, 3767000), ylim = c(2575500, 2647200), clip = "off") +
  theme_void()+
  theme(
    text = element_text(size = 18)
  )
```

We hope that future work allows us to address the unsupported features of VDC maps for standard plotting methods like base R, `ggplot2` and `tmap`
and that we can add more visualisation options like the multi-dimensional plot.


### Time series plots

As we are working with spatio-temporal data, to visualise the time component comes as natural.
We can extract information both from the array and tabular format to do so.

```{r}
par(mar = c(4,5,1,1))
plot(
  cube_arr_ldsl$area,
  col = rep(c('deepskyblue3','peru'), 10),
  pch = 19
)
```

If we want to use tidyverse methods, the tabular VDC presents a more straightforward syntax.

```{r}
cube_tab_ldsl |> 
  face_temporal() |> 
  ggplot() +
  aes(
    x = date, y = area,
    color = date,
    shape = sensor, 
    group = class
  ) +
  geom_point() + geom_line() +
  scale_color_viridis_c(
    "Date", direction = 1,
    trans = "date", option = "D"
  ) 
```

### Glyph maps

`cubble` includes glyph maps as one of their main features. 
The idea is to summarise time series data in a map format.
The concept works as an innovative way to communicate spatio-temporal variation.

Here we have used it to show the landslides and lake area time series.
The glyphs are located in the `geom_sum` coordinates.
This type of visualisation becomes handy when dealing with several feature sets in an area and to have a fast overview of large changes per location.
We included the landslide and lake delineations as background for reference.

```{r}
cube_tab_ldsl |> 
  face_temporal() |> 
  # make the x and y coordinates explicit on the temporal table 
  # to pass onto ggplot aesthetics
  unfold(x, y) |> 
  ggplot() +
  # add delineations as background for reference
  geom_sf(
    data = st_as_sf(face_temporal(cube_tab_ldsl), sf_column_name = "geom"),
    inherit.aes = F, fill = NA,
    aes(color = class), 
    show.legend = T
  ) +
  scale_color_manual("Landform type", values = c("deepskyblue3", "peru")) +
  new_scale_color() +
  # map aesthetics for glyph maps
  aes(x_major = x, y_major = y,
      x_minor = date, 
      y_minor = unclass(area)) +
  # plot glyphs
  geom_glyph_box(
    width = 510, height = 400,
    color = "black", linewidth = 0.1
  ) +
  geom_glyph(
    width = 500,  height = 390,
    aes(color = drop_units(set_units(area, "km^2")))
  ) +
  scale_color_viridis_c(
    expression("Area [ha]"),
    direction = 1, option = "D"
  ) +
  # bring CRS back
  coord_sf(
    crs = st_crs(cube_tab_ldsl)
  ) +
  theme_void() 
```

```{r}
par(oldpar)
```

# References
::: {#refs}
:::

# R Session Info {#sec-session}

```{r}
sessioninfo::session_info()
```