Mapping the Arctic

Robert Schlegel

Problem

How do we plot maps near the poles more accurately?
What effect does this have on other map layers?
This is time intensive, is there a faster way?

Solution

We must learn more about map projections
We will look at packages that help us project raster layers etc.
ggOceanMaps provides the full suite of solutions

Setup

We will use many libraries and datasets

# Load libraries
library(tidyverse)
library(ggOceanMaps)
library(ggOceanMapsData)

# NB: ggOceanMapsData is not on CRAN
  # Uncomment and run this line of code:
# remotes::install_github("MikkoVihtakari/ggOceanMapsData")
  # If this causes an error, think about why

# Fixed base map
map_global_fix <- map_data('world') %>% 
  rename(lon = long) %>% 
  mutate(group = ifelse(lon > 180, group+2000, group),
         lon = ifelse(lon > 180, lon-360, lon))

# Load SST from 2022-12-25 to 2022-12-31
load("../data/OISST_2022.RData")

The Earth is round

No matter what you’ve heard, its not flat. Proof of this is that if it were flat, it would be much easier to plot parts of it on a map! Instead, because the surface of the Earth is curved, we cannot accurately plot it on a computer screen.

The way we deal with this is by changing the measurements between longitude/latitude coordinates. The technical word is ‘projecting’. The noun is ‘projection’ or ‘projections’.

For much of the Earth this is not that noticeable. But for the Arctic it is a big issue.

Map projections

There are many different sorts of map projections in use. Over tea time I recommend googling them. There are some interesting ones.

The Spilhaus projection is oriented around the global ocean.

Map projections

R has many built-in projections, and ggplot2 gives us easy access to many of them.

# Type the following in your console,
# see what comes up as auto-complete suggestions
# Use the up and down arrows to scroll through
coord_

Cartesian projections

This is the default coordinate projection
Basically, it doesn’t alter the lon/lat decimal degree coordinates in any way

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  # Numeric sizing for lon/lat 
  coord_cartesian()

Equal projections

This projects numeric x/y (e.g. lon/lat) values the same

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  # Equal sizing for lon/lat 
  coord_equal()

Fixed projections

Similar to Fixed projections, but let’s the user decide on the ratio.

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  # Ratio (Y divided by X) sizing for lon/lat 
  coord_fixed(ratio = 2)

Map projections

There is also a quick mapping option
This forces straight lines, but runs quickly

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  # Behind the scenes this adapts the "mercator" projection
  coord_quickmap()

Map projections

We access crs map projections with the _sf() range of functions
crs = Coordinate Reference System
These are supported with the packages sf and stars
Which require additional software outside of R

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  coord_sf() # sf = simple feature

Polar projections

coord_polar() does not quite give us what we want
More designed for bar plots etc.

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  scale_y_reverse() +
  # A very different projection
  coord_polar()

Polar projections

We can use the ‘ortho’ projections with coord_map()

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  # Look up the help file for moer info
  coord_map(projection = "ortho", orientation = c(90, 0, 0))

Projecting layers

Don’t try running this, the background calculations are enormous

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  geom_tile(data = OISST_2022, aes(fill = temp)) +
  coord_map(projection = "ortho", orientation = c(90, 0, 0))

Projecting layers

Making the plot smaller can help

# Filter map data and plot it in one code chunk
map_global_fix %>% 
  filter(lon > 9, lon < 28, lat > 76, lat < 81) %>% 
  ggplot(aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  # Filtering the OISST_2022 data directly in geom_tile()
  geom_tile(data = filter(OISST_2022,
                          lon > 9, lon < 28, lat > 76, lat < 81), 
            aes(fill = temp)) +
  coord_map(projection = "ortho", orientation = c(90, 0, 0))

Projecting layers

Or plot it on a flat surface

# Filter map data and plot it in one code chunk
map_global_fix %>% 
  filter(lon > 9, lon < 28, lat > 76, lat < 81) %>% 
  ggplot(aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group)) +
  # Filtering the OISST_2022 data directly in geom_tile()
  geom_tile(data = filter(OISST_2022,
                          lon > 9, lon < 28, lat > 76, lat < 81), 
            aes(fill = temp)) +
  coord_cartesian(expand = F)

All-in-one

Better yet, just use ggOceanMaps
Click the package name to see the intro tutorial

basemap(limits = 60)

`ggOceanMaps`

Works everywhere

basemap(limits = c(100, 160, -20, 30), bathymetry = TRUE)

`ggOceanMaps`

Has glaciers, too

basemap(limits = 60, glaciers = TRUE, bathymetry = TRUE)

`ggOceanMaps`

Fundamentally designed to plot the Arctic

basemap(limits = c(-160, -80, 60, 85), rotate = TRUE)

Citing packages

When a package is very important to our workflow it is good to cite it
Especially when the author(s) requests that we do so

citation("ggOceanMaps")


To cite package 'ggOceanMaps' in publications use:

  Vihtakari M (2022). _ggOceanMaps: Plot Data on Oceanographic Maps
  using 'ggplot2'_. R package version 1.3.4,
  <https://CRAN.R-project.org/package=ggOceanMaps>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {ggOceanMaps: Plot Data on Oceanographic Maps using 'ggplot2'},
    author = {Mikko Vihtakari},
    year = {2022},
    note = {R package version 1.3.4},
    url = {https://CRAN.R-project.org/package=ggOceanMaps},
  }