manuscript-figures.qmd

---
format: 
  gfm: default
  html:
    toc: true
    toc-depth: 3
    number-sections: true
    fig-dpi: 300
---

# Install packages
Run this chunk to get the necessary packages. 
Note that, to reproduce, some packages require the development version fetched from GitHub.

```{r}
#| eval: false
#| message: false
#| warning: false
pkgs = c(
  "dplyr","ggplot2","here","patchwork",
  "purrr","readr","stringr",
  "tidyr","units","zen4R"
)
install.packages(pkgs)
remotes::install_github("huizezhang-sherry/cubble")
remotes::install_github("eliocamp/ggnewscale")
remotes::install_github("r-spatial/sf")
remotes::install_github("r-spatial/stars")
```

# Libraries 
```{r}
#| message: false
#| warning: false
library(cubble)
library(dplyr)
library(ggnewscale)
library(ggplot2)
library(here)
library(patchwork)
library(purrr)
library(readr)
library(sf)
library(stars)
library(stringr)
library(tidyr)
library(units)
library(zen4R)
```

# Figure 1 
Created with sketchpad.io and draw.io
See source file in code/vdc-schema.drawio

# Figure 2
Code snippet with the result of the following script
The script creates a VDC for Fagradalsfjall
It is in array format, using the `stars` package
 
Data comes from Zenodo so the script shows demo code with steps to
1. download to temporary directory
2. pre-process
3. coerce to a vdc

## Download
```{r}
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
)
```

## Pre-process
```{r}
# 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"
      )
    ) |> 
    st_set_crs(3057)
}

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

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

## VDC coercion
```{r}
# Create array
a = array(
  data = outlines$geom, 
  dim = c(
    1,
    length(unique(outlines$datetime))
  ),
  dimnames = list(
    geom_sum = 1,
    datetime = unique(outlines$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$geom))), # approach centroid
  datetime = unique(outlines$datetime),
  point = c(TRUE, FALSE)
)

# Coerce to cube
# The output of this code constitutes Figure 2
(cube_arr = st_as_stars(
  list(geometry = a), 
  dimensions = dim_cent)
)
```

# Figure 3
Base R plot of lava flow outlines
a) shows all the lava flows for the whole time period
b) shows the result of filtering the VDC between two dates
Plots are combine with the patchwork library,
hence the use of `wrap_elements()`

```{r}
#| fig-width: 4.75
#| fig-height: 8
#| fig-dpi: 300
oldpar = par(no.readonly = TRUE)
par(mar = c(1.8,1,1.1,0.3), bg = "transparent")
wrap_elements(
  ~cube_arr |>
    # use lambda-function to extract geometries from vdc attributes
    (\(x) plot(x$geometry, axes = TRUE))(),
  clip = FALSE
) /
wrap_elements(
  ~cube_arr |>
    filter(datetime > "2021-03-18", datetime < "2021-03-25") |> 
    # use lambda-function to extract geometries from vdc attributes
    (\(x) plot(x$geometry, axes = TRUE))(),
  clip = FALSE
) +
  plot_annotation(tag_levels = "a", tag_suffix = ")")
par(oldpar)
```

# Figure 4
Code snippet with the result of the following script
The script creates a VDC for Butangbunasi
It is in tabular format, using the `cubble` package

Data comes from Zenodo so the script shows demo code with steps to
1. download to temporary directory
2. pre-process
3. coerce to a vdc

## Download
```{r}
# Fetch also the CSV file with addiontional 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
)

# 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
mapping = lapply(mapping_ls, read_fun) |> 
  bind_rows() |> 
  left_join(stats, by = c("fn_gpkg" = "file"))
```

## Pre-process
```{r}
# All pre-processing steps are chained but have 
# documentation in each line
landslides = mapping |> 
  # 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 and class 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
  ungroup() |> 
  # set area to respective class
  mutate(
    area = case_when(
      Class == "landslide" ~ area_ldsl,
      Class == "lake" ~ area_lake
    )
  ) |> 
  # remove wide area columns
  select(-starts_with("area_")) |> 
  # complete cases for every combination of date and class
  # for geometries this is filled with an empty geometry
  complete(date, Class) |> 
  # rename the class column (lowercase)
  rename(class = Class) |> 
  # convert back to sf
  st_as_sf()

# Compute geom_sum as the centroid of the union of 
# the geometries per feature class
landslidescent = landslides |> 
  group_by(class) |> 
  mutate(geom_sum = st_centroid(st_union(geom))) |> 
  relocate(geom_sum, .after = everything()) |> 
  ungroup() |> 
  st_as_sf(sf_column_name = 'geom_sum')
```

## VDC coercion
```{r}
# Call as cubble to coerce from sf to cubble
# key is the gid column while index is date
cube_tab = as_cubble(
  landslidescent, key = class, index = date
)

# The output of this code corresponds to Figure 4a
cube_tab |> 
  face_spatial()
# The output of this code corresponds to Figure 4b
cube_tab |> 
  face_temporal() |> 
  arrange(date) |> 
  print(n = 4)
```

# Figure 5
This figure shows in 
a) a multi-dimensional representation
of the Butangbunasi landslide and 
b) a time series of the landslide area
To create a) the geometry of the original landslide data
(outside of the cube) was distorted by multiplying it
by a shear matrix and adding a shift in the y axis to 
stack them on top of each other

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

# Apply shear matrix
ldsl_shear = landslides |>
  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(landslides)) 

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

# Plot a)
spaceplot = 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(2576000, 2646000), clip = "off") +
  theme_void()+
  theme(
    text = element_text(size = 18)
  )

# Wrangle data to add in rug of area plot
rug_data = cube_tab |> 
  face_temporal() |> 
  filter(class == 'lake', !st_is_empty(geom)) |> 
  mutate(
    class = "Landslide-dammed lake",
    Area = 0.3
  )

# Plot b), directly from cubble object in temporal face
area = cube_tab |> 
  face_temporal() |> 
  rename(Area = area) |> 
  filter(class == 'landslide') |> 
  ggplot() +
  aes(
    x = date, y = Area,
    shape = sensor, 
    group = 1
  ) +
  geom_point(size = 3, aes(color = date)) + 
  geom_line(aes(color = date)) +
  scale_y_units(unit = "km2") +
  scale_color_viridis_c(
    "Date", direction = 1,
    trans = "date", option = "D",
    guide = guide_colorbar(barwidth = 20, barheight = 0.8)
  ) +
  scale_shape("Sensor") +
 ggnewscale::new_scale_color() +
  geom_rug(data = rug_data, aes(color = class), sides = "b") +
  scale_color_manual("", values = "red") +
  theme_bw() +
  theme(
    text = element_text(size = 18),
    axis.title.x = element_blank(),
    legend.box = "horizontal",
    legend.position = "bottom",
    legend.title.position = "top"
  )

# Set-up layout
layout = "
AAAABBBBBB
AAAABBBBBB
AAAABBBBBB
AAAABBBBBB
AAAABBBBBB
AAAABBBBBB
AAAABBBBBB
AAAABBBBBB
AAAABBBBBB
AAAACCCCCC
AAAACCCCCC
AAAACCCCCC
AAAACCCCCC
AAAACCCCCC
"

# Combine plots with patchwork
spaceplot + area + guide_area() + 
  plot_layout(guides = "collect", design = layout) + 
  plot_annotation(tag_levels = "a", tag_suffix = ")")
```

# R Session Info {#sec-session}

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