More ggplot2

library(tidyverse)
library(patchwork)

Poor-man’s gridding

ICES rectangles

If ICES rectangles is not a variable in your dataset there are some functions floating around that allows one to assign a geographical position to a rectangle without converting things first to an sf-object. One sweep of such functions is located in the geo-package. The package has been thrown off cran, mostly because improper bookkeeping of c++-functions given the latest standard. So to install it we need:

library(devtools)
install_github("hafro/geo")

library(geo)

Importing data used:

minke <- 
  read_csv("ftp://ftp.hafro.is/pub/data/csv/minke.csv")

To assign a position to an ICES rectangle we use the d2ir-function.
To get the midpoint of the ICES rectangle we use the ir2d-function.

minke <- 
  minke %>% 
  select(lon, lat) %>% 
  mutate(ir = d2ir(lat, lon),
         ir_lon = ir2d(ir)$lon,
         ir_lat = ir2d(ir)$lat) 
minke %>% glimpse()

## Rows: 190
## Columns: 5
## $ lon    <dbl> -21.42350, -21.39183, -19.81333, -21.57500, -15.61167, -18.6700…
## $ lat    <dbl> 65.66183, 65.65350, 66.51167, 65.67333, 66.29000, 66.17833, 65.…
## $ ir     <chr> "60C8", "60C8", "62D0", "60C8", "61D4", "61D1", "60C8", "61C7",…
## $ ir_lon <dbl> -21.5, -21.5, -19.5, -21.5, -15.5, -18.5, -21.5, -22.5, -17.5, …
## $ ir_lat <dbl> 65.75, 65.75, 66.75, 65.75, 66.25, 66.25, 65.75, 66.25, 66.75, …

Now we can do some summary statistics by ICES rectangles (actually the mid-points of lon and lat) and visualize:

minke.sum <- 
  minke %>% 
  group_by(ir_lon, ir_lat) %>% 
  summarise(n = n()) %>% 
  ungroup() 

minke.sum %>% 
  ggplot() +
  geom_tile(aes(x = ir_lon, y = ir_lat, fill = n)) +
  geom_text(aes(x = ir_lon, y = ir_lat, label = n), angle = 45, colour = "red") +
  coord_quickmap() +
  labs(x = NULL, y = NULL)

Home-made gridding function

It is easiest to use functions in the raster-package for setting spatial data on a grid. There are however other means, here we demonstrate using a home-made function (so use with caution).

#' grade
#'
#' @param x A numerical vector to set on a grid
#' @param dx The resolution (NOTE: not tested for values greater than 1)
#'
#' @return A vector of grid midpoint values

grade <- function(x, dx) {

  if(dx > 1) warning("Not tested for grids larger than one")
  brks <- seq(floor(min(x)), ceiling(max(x)),dx)
  ints <- findInterval(x, brks, all.inside = TRUE)
  x <- (brks[ints] + brks[ints + 1]) / 2
  return(x)
}

This function in action:

minke %>% 
  select(lon, lat) %>% 
  mutate(glon = grade(lon, 0.50),
         glat = grade(lat, 0.25)) %>% 
  glimpse()

## Rows: 190
## Columns: 4
## $ lon  <dbl> -21.42350, -21.39183, -19.81333, -21.57500, -15.61167, -18.67000,…
## $ lat  <dbl> 65.66183, 65.65350, 66.51167, 65.67333, 66.29000, 66.17833, 65.66…
## $ glon <dbl> -21.25, -21.25, -19.75, -21.75, -15.75, -18.75, -21.25, -22.75, -…
## $ glat <dbl> 65.625, 65.625, 66.625, 65.625, 66.375, 66.125, 65.625, 66.125, 6…

Effectively we have binned the:

longitude to 0.50 degrees and provided the midpoint (.25 or .75)
latitude to 0.25 degress and provided the midpoint (.125, .375, …)

Lets repeat the tile-plot we did in the first ggplot but using double as high resolution as the ICES rectangle:

smb2019 <- 
  read_csv("ftp://ftp.hafro.is/pub/data/csv/is_smb.csv") %>% 
  filter(year == 2019)
smb2019.g <- 
  smb2019 %>% 
  mutate(glon = grade(lon1, 0.50),
         glat = grade(lat1, 0.25)) %>% 
  group_by(glon, glat) %>% 
  summarise(Cod = mean(cod_n)) %>% 
  ungroup()
p <- 
  smb2019.g %>% 
  ggplot() +
  geom_tile(aes(x = glon, y = glat, fill = Cod)) +
  coord_quickmap() +
  labs(x = NULL, y = NULL)

p +
  scale_fill_viridis_c(option = "B", direction = -1)

While we are at it, just as we can control the values on the x- and y-axis we can set the label breaks on colours:

p +
  scale_fill_viridis_c(option = "B", direction = -1,
                       breaks = seq(0, 1500, by = 250))

Some more gritty details on poor-mans gridding can be found here

Mapping colours to values

Built in transformations

In ggplot there are a bunch of inbuilt transformation (“asn”, “atanh”, “boxcox”, “date”, “exp”, “hms”, “identity”, “log”, “log10”, “log1p”, “log2”, “logit”, “modulus”, “probability”, “probit”, “pseudo_log”, “reciprocal”, “reverse”, “sqrt” and “time”) that may be useful when plotting spatial data mapped to colour. Take e.g. these two examples:

p1 <- 
  p +
  scale_fill_viridis_c(option = "B", direction = -1,
                       trans = "log")
p2 <- 
  p +
  scale_fill_viridis_c(option = "B", direction = -1,
                       trans = "sqrt")
p1 + p2

Using discreet value

One may often want to generate a plot that shows a range of some continuous values as discrete interval. For that one could e.g. use the cut-function. Here is an example were based on intervals that such that the number in each interval are roughly the same:

p <-
  smb2019.g %>% 
  mutate(Cod = cut(Cod,
                    breaks =  c(0, 23, 62, 145, 1600),
                    include.lowest = TRUE)) %>% 
  ggplot() +
  geom_tile(aes(x = glon, y = glat, fill = Cod)) +
  coord_quickmap() +
  labs(x = NULL, y = NULL) 
p +
  scale_fill_viridis_d(option = "B", direction = -1)

# Try a softer colour palette, check also tcrenv

Formatting axis values

If you only have pure tbl-object, i.e. your object is not sf and/or you are not using the geom_sf-function) you can format the axis using these home-made functions (if there is a bug send me an e-mail):

# https://stackoverflow.com/questions/33302424/format-latitude-and-longitude-axis-labels-in-ggplot
scale_longitude <- function(min = -180, max = 180, step = 1, ...) {
  breaks <- seq(min, max, step)
  labels <- 
    ifelse(breaks < 0,
           paste0(breaks, "\u00B0", "W"),
           paste0(breaks, "\u00B0", "E"))
  return(scale_x_continuous(name = NULL, breaks = breaks, labels = labels, ...))
}
scale_latitude <- function(min = -90, max = 90, step = 0.5, ...) {
  breaks <- seq(min, max, step)
  labels <- 
    ifelse(breaks < 0,
           paste0(breaks, "\u00B0", "S"),
           paste0(breaks, "\u00B0", "N"))
  return(scale_y_continuous(name = NULL, breaks = breaks, labels = labels, ...))
}

p1 <- 
  ggplot(geo::island, aes(lon, lat)) +
  geom_path() +
  coord_quickmap()
p2 <-
  p1 +
  scale_longitude(step = 2) +
  scale_latitude(step = 1)

p1 + p2

ICES rectangles

Here is an example if you want the ICES rectangles as a label on the x- and the y-axis:

scale_longitude_ices <- function(min = -44, max = 68.5, step = 1, ...) {
  breaks <- seq(min + 0.5, max - 0.5, step)
  labels <- geo::d2ir(60, breaks) %>% str_sub(3)
  return(scale_x_continuous(name = NULL, breaks = breaks, labels = labels, ...))
}
scale_latitude_ices <- function(min = 36, max = 84.5, step = 0.5, ...) {
  breaks <- seq(min + 0.25, max - 0.25, step)
  labels <- geo::d2ir(breaks, 0) %>% str_sub(1, 2)
  return(scale_y_continuous(name = NULL, breaks = breaks, labels = labels, ...))
}

ggplot(geo::island, aes(lon, lat)) +
  geom_path() +
  coord_quickmap() +
  geom_path() +
  scale_longitude_ices() +
  scale_latitude_ices()

A detailed formatting example:

ir <- 
  sf::read_sf("ftp://ftp.hafro.is/pub/data/shapes/ices_rectangles.gpkg") %>% 
  filter(between(west, -25, -14),
         between(south, 63, 66))
ggplot() +
  #geom_sf(data = ir) +
  geom_path(data = geo::island, aes(lon, lat),
            colour = "black") +
  geom_text(data = ir,
            aes(x = west + 0.5, y = south + 0.25,
                label = icesname),
            angle = 45, 
            colour = "red",
            size = 3) +
  scale_longitude_ices() +
  scale_latitude_ices() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_line(size = 1),
        axis.ticks = element_blank()) +
  coord_quickmap()

ggplot2 details

If you dig deep into ggplot you can set up subplot within your spatial plot, here an example that uses pure tibble-approach. What is shown is the abundance trend by year (1985-2019) of the cod in the Icelandic spring survey:

smb <- 
  read.csv("ftp://ftp.hafro.is/pub/data/csv/is_smb.csv",
           stringsAsFactors = FALSE) %>% 
  as_tibble()

library(GGally)
by.rect <- 
  smb %>% 
  #filter(ir_lon <= -22, ir_lat >= 66) %>% 
  rename(n = cod_n) %>% 
  mutate(n = ifelse(n > quantile(n, 0.99), quantile(n, 0.99), n)) %>% 
  group_by(year, ir_lon, ir_lat) %>% 
  summarise(n = mean(n))
n.glyph <-
  by.rect %>% 
  #mutate(n = ifelse(n > 0 & n > quantile(n, 0.40), quantile(n, 0.40), n)) %>% 
  glyphs(x_major = "ir_lon", 
         y_major = "ir_lat",
         x_minor = "year", 
         y_minor = "n", 
         width = 1, 
         height = 0.5)


tows <- 
  smb %>% 
  filter(year == 2019)

n.glyph %>% 
  mutate(pos = ifelse(n != 0, TRUE, FALSE),
         base = ir_lat - 0.25,
         gy = ifelse(n == 0, gy + 0.005, gy)) %>% 
  ggplot() +
  theme_bw() +
  geom_linerange(aes(x = gx, ymin = base, ymax = gy),
                 colour = "black") +
  geom_segment(data = tows,
                aes(x = lon1, y = lat1, xend = lon2, yend = lat2),
                colour = "red") +
  geom_path(data = geo::island, aes(lon, lat)) +
  coord_quickmap() +
  scale_longitude_ices() +
  scale_latitude_ices() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_line(size = 1),
        axis.ticks = element_blank())

More ggplot2

Various tidbits

Poor-man’s gridding

ICES rectangles

Home-made gridding function

Mapping colours to values

Built in transformations

Using discreet value

Formatting axis values

ICES rectangles

ggplot2 details