palaeomap_presentation.Rmd

---
title: "Plotting Palaeogeographical Maps in R: an Example"
author: Ben M
date: 2021-03-11
fig_caption: true
output: beamer_presentation
---

# Plotting Palaeogeographical Maps in R

- Run in the cloud on Binder: <https://bit.ly/2N86wYf>
- _OR_
- Download code from GitHub: <https://bit.ly/38tBH7U>


# Palaeogeographical Maps

>- Hand made maps are a drag
>- Deep Time Maps (<https://deeptimemaps.com>)
>- Paleomap Project (<http://scotese.com>)

***

**GPlates** (<https://www.gplates.org/>) reconstructs palaeogeography at arbitrary
time.

# Automatic map plotting

>- NonaR/paleoMap (<https://github.com/NonaR/paleoMap>)
>- LunaSare/gplatesr (<https://github.com/LunaSare/gplatesr>)

```{r setup, echo = FALSE}
library(broom)
library(ggthemes)
library(mapproj)
library(rgdal)
library(tidyverse)

list.files("functions/", full.names = TRUE, pattern = "\\.R") %>%
  walk(source)
```

# GPlates Web Service

- <https://gws.gplates.org/>

```{r gws_data, eval = FALSE}
coastline_gws_url <-
  "http://gws.gplates.org/reconstruct/coastlines/?time=155&model=GOLONKA"
polygons_gws_url <-
  "http://gws.gplates.org/reconstruct/static_polygons/?time=155&model=GOLONKA"

kimmeridgian_coastlines <-
  rgdal::readOGR(coastline_gws_url) %>%
    broom::tidy()
kimmeridgian_polygons <-
  rgdal::readOGR(polygons_gws_url) %>%
    broom::tidy()
```

***

```{r gws_plot, echo = FALSE, fig.width = 7, fig.height = 4, fig.cap = "**Outlines of continental plates in the Kimmeridgian (155 Ma).** Data downloaded from the GPlates web service.", warning = FALSE, message = FALSE, results = "hide"}
kimmeridgian_coastlines <-
  readOGR(
    "data/GWS/Matthews_etal_GPG_2016_Coastlines_reconstructed_155.00Ma.gmt"
  ) %>%
  tidy()
kimmeridgian_polygons <-
  readOGR(
    "data/GWS/Matthews_etal_GPG_2016_Polygons_reconstructed_155.00Ma.gmt"
  ) %>%
  tidy()
ggplot() +
  geom_map(
    data = kimmeridgian_polygons,
    map = kimmeridgian_polygons,
    aes(x = long, y = lat, map_id = id),
    fill = "#D8D8D8"
  ) +
  geom_map(
    data = kimmeridgian_coastlines,
    map = kimmeridgian_coastlines,
    aes(x = long, y = lat, map_id = id),
    colour = "#222222", fill = NA, size = 0.3
  ) +
  coord_map("mollweide") +
  map_border() +
  theme_map()
```


# 'True' Palaeogeographical Outlines

GPlates web service gives:

>- outlines of plates from the model
>- modern coastline outlines
>- I want _ancient coastal outlines_ instead

***

- GPlates software has more models
- Including Cao et al. (2017) palaeogeographical reconstructions

```{r map_layers}
map_layers <-
  c(
    "Land"              = "#FFD23A",
    "Mountain"          = "#FF8D51",
    "Shallow marine"    = "#45D8FF"
  )
```

# Method

1. Load data into GPlates
2. Export for desired
3. Load into R and plot

```{r read_polygon_data, include = FALSE}
polygon_data <-
  list.files("data/palaeogeography/", pattern = ".gmt", full.names = TRUE) %>%
    map(readOGR) %>%
    map(tidy) %>%
    map2(
      names(map_layers),
      ~ add_column(.x, geog_layer = .y)
    ) %>%
    bind_rows() %>%
    mutate(
      polygon_id = str_c(id, group, geog_layer) %>% as_factor(),
      geog_layer = factor(geog_layer, labels = names(map_layers))
    )
```

***

But showing the separate layers (land, mountain, shallow marine, ice) doesn't
work easily with `geom_map`.

- have to add a new geom for every layer:

```{r multiple_geom_map, eval = FALSE}
ggplot() +
  geom_map() +
  geom_map() +
  geom_map() + …
```

***

Instead I assign names to the layers (land, mountain…) and use `geom_polygon` to
plot.

NB:

- The data has 'groups' and 'subgroups' to close off the polygons, otherwise
  shapes may cross the whole globe.
- Layers must be plot in order be marine > land > mountain > ice – use factors.

```{r rearrange_polygon_labels, include = FALSE}
id_level_order <-
  map_layers[c(3, 1, 2)] %>%
    map(
      function(layr) polygon_data$polygon_id %>% levels() %>% str_which(layr)
    ) %>%
    unlist()

polygon_data <-
  polygon_data %>%
    mutate(
      polygon_id = fct_relevel(polygon_id, levels(polygon_id)[id_level_order])
    )
```

***

```{r plot_base_map, fig.width = 7, fig.height = 4, fig.cap = "**Palaeogeography in the Kimmeridgian (155 Ma).**", echo = FALSE}
# dev.new(width = 7, height = 4)
map_plot <-
  ggplot() +
    geom_polygon(
      data = polygon_data,
      aes(
        x      = long,
        y      = lat,
        fill   = geog_layer,
        colour = geog_layer,
        group  = polygon_id
      ),
    ) +
    coord_map("mollweide") +
    theme_map() +
    palaeogeog_map_niceties()
map_plot
```

***

I also wanted to add coastlines to show where modern countries are, but not the
full plates from the web service.

>- use coastlines reconstructed polylines

```{r read_coastline_polygons, echo = FALSE, results = "hide"}
modern_coastlines <-
    readOGR("data/coastlines/Matthews_etal_GPC_2016_Coastlines_Polyline_reconstructed_155.00Ma.gmt") %>%
    tidy() %>%
    add_column(geog_layer = "Modern coastlines") %>%
    mutate(polygon_id = str_c(id, group, geog_layer) %>% as_factor())
```

***

```{r plot_coast_outlines, fig.width = 7, fig.height = 4, fig.cap = "**Palaeogeography in the Kimmeridgian (155 Ma).** Outlines of modern coastlines (where known) are included in grey.", echo = FALSE}
map_plot +
  geom_path(
    data = modern_coastlines,
    aes(
      x = long,
      y = lat,
      group = polygon_id,
    ),
    colour = "#888888", size = 0.3
  )
```


# Add Fossil Occurrences

Now to add some fossil occurrences.

Of course it's ichthyosaurs, because I have no imagination:

- Callovian–Tithonian (166-145 Ma)
- all levels, taxonomy – no filtering

```{r download_pbdb_data, eval = FALSE}
pbdb_url <-
  "https://paleobiodb.org/data1.2/occs/list.csv?base_name=Ichthyosauromorpha&interval=Callovian,Oxfordian,Kimmeridgian,Tithonian&show=paleoloc"

occ_ichthyosaurs <-
  read_csv(pbdb_url)
```

Use `geom_point` to overlay this on the base map.

```{r read_ichthyosaur_occurrences, echo = FALSE, results = "hide", warning = FALSE, message = FALSE}
occ_ichthyosaurs <-
  read_csv("data/occurrences/2021-02-18-ichthyosaur-occurrences.csv")
```

***

```{r plot_occurrences, fig.height = 4, fig.width = 7, fig.cap = "**Occurrences of ichthyosaurs from the Callovian–Tithonian.** Palaeogeographical map shows the distribution of land in the Kimmeridgian (155 Ma).", echo = FALSE, warning = FALSE}
occ_plot <-
  map_plot +
    geom_point(
      data = occ_ichthyosaurs,
      aes(
        x = paleolng,
        y = paleolat
      )
    )
occ_plot
```

***

These plots can also be split automagically using _facets._

```{r facet_plot, fig.width = 7, fig.height = 4, fig.cap = "**Occurrences of ichthyosaurs in the Callovian–Tithonia separated by identified rank.** Palaeogeographical maps shows distribution of land in the Kimmeridgian (155 Ma).", echo = FALSE, warning = FALSE}
occ_plot +
  facet_wrap(vars(accepted_rank))
```