Interactive plots - leaflet

library(rnaturalearth)
# need to install below if not already on your computer
# install.packages("rnaturalearthhires", repos = "http://packages.ropensci.org", type = "source")
library(sf)
library(htmltools)
library(leaflet)
library(tidyverse)

In the last decade there has been a proliferation of packages that allow one to generate JavaScript applications using R syntax and data. The interactive visualizations can be used in markdown reports, presentations, Shiny applications or as standalone webpages.

leaflet

Introduction

Package homepage

The basic construct

Of the various interactive map packages available in R, the {leaflet}-package is the most mature in terms of features and controls. Generating a leaflet map is as simple as:

leaflet() %>% 
  addTiles() %>% 
  addPopups(lng = -59.58779609159309,
            lat = 13.07416042366417, 
            popup = "Blue Horizon")

Here:

• a map widget is created by calling leaflet
• a background map is created by calling "addTiles
• a spatial point is added, here using addPopups

Leaflet works similar as {ggplot2}, in that one can compile multiple layers by calling add****-functions. The difference, besides difference in the layer-function names, is that in {ggplot2} we use the “+”, while in {leaflet} one uses “%>%”.

Exercise

Repeat the above command using your favorite location

Data input

{leaflet} can take as data input one of the following forms:

• data.frames containing coordinates
• spatial objects objects such as points, lines and polygons
• objects of class “map” as return from the map-package

Lets use the minke-dataset to demonstrate the basic leaflet components:

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

A simple plot of the sample locations can be obtained by using the addCircles-function:

leaflet() |> 
  addTiles() |>  
  addCircles(data = w,
             lng = ~lon,
             lat = ~lat,
             radius = 10000)  # radius in meters

If the coordinates in your dataframe are labelled “lon”, “lng”, “long”, or “longitude” and “lat” or “latitude” we do need to declare them in the arguement:

leaflet() |> 
  addTiles() |>  
  addCircles(data = w,
             radius = 10000)  # radius in meters

Also if the object is an sf-object:

w_sf <- 
  w |> 
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326)

we do not need to specify the xy-coordinates in the function call:

leaflet() %>% 
  addTiles() %>% 
  addCircles(data = w_sf,
             radius = 10000)  # radius in meters

Layer calls and arguements

The syntax used to call a variable in leaflet differs from that done in {ggplot2}, with {leaflet} one always needs the prefix “~” in front of the variable name. And, instead of having the x,y argument pairs in {ggplot2}, in {leaflet} the pairs are lng,lat.

The main layer functions in {leaflet} and corresponding sister argument in {ggplot2} are:

leaflet       ggplot2
-------       -------
addCircles    geom_point
addPolylines  geom_path
addPolygons   geom_polygon

The main argument in {leaflet} and corresponding sister argument in {ggplot2} are:

leaflet               ggplot2
-------               -------
lng                   x
lat                   y
color                 color, colour
fillColor             fill
weight                lwd
radius                size
opacity, fillOpacity  alpha       
label                 label
dashArray             type

Like in {ggplot2} these arguments can either be a fixed value or one may attempt to use a value of a variable in the data to dictate the aesthetics.

Exercise

Emulate the following plot where:

• The arguement radius is controlled by the values of the stomach.volume
• The arguement label is controlled by sex.
• The arguement weight (controls the appearance of the “outer ring”) is overridden/suppressed.
• The default colour value (“#03F”) is overridden (you can use standard English R colour names)

As usual it is best to read the help-file and/or check what are the arguments for the addCirles()-function (call args(addCircle)).

Mapping variables to colours

Lets try to distinguish the samples by sex using the argument color. A spoiled ggplotter may try this:

w %>% 
  leaflet() %>% 
  addTiles() %>% 
  addCircles(lng = ~lon,
             lat = ~lat,
             color = ~sex,
             radius = 10000)  # radius in meters

The above will however not work as anticipated because {leaflet} is not as user-friendly as {ggplot2} that provides some default colouring for free.

We hence need first to specify a colour palette (more on this later):

pal <- colorFactor(palette = c("navy", "red"), 
                   domain = c("Male", "Female"))

and then:

m <- 
  w %>% 
  leaflet() %>% 
  addTiles() %>% 
  addCircles(lng = ~lon,
             lat = ~lat,
             color = ~pal(sex),
             radius = 10000)
m

This plot however does not give us indication about which colour corresponds to which sex. Again, legends that come free in ggplot need to be specifically provided in leaflet using the addLegend-function:

m %>% 
  addLegend("bottomleft",
            title  = 'Sex',
            pal = pal,
            values = ~sex)

If the object is of class sf (or sp) one does not need to specify the lng or lat (so same as when generating sf-plots or geom_sf-plots):

wsf <- 
  w %>% 
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326,
           remove = FALSE)

# NOT RUN
wsf %>% 
  leaflet() %>% 
  addTiles() %>% 
  # note: no arguement lng nor lat used
  addCircles(color = ~pal(sex),
             radius = 10000) %>% 
  addLegend("bottomleft",
            title  = 'Sex',
            pal = pal, 
            values = ~sex)

As with ggplot, one can specify the data used from within a layer-call:

# NOT RUN
leaflet() %>% 
  addTiles() %>% 
  addCircles(data = w,
             color = ~pal(sex),
             radius = 10000) %>% 
  addLegend("bottomleft",
            title  = 'Sex',
            pal = pal, 
            values = c("Male", "Female"))

Note however that one needs to specify explicitly the “domain” within the addLegend-function (NOTE: check if there is not an alternative to this).

Basemaps

The addTiles-function provides a reference background to spatial maps, the default being openstreetmap-tiles. That background may not be the most suitable for oceanographic data we are most likely confronted with. The addProviderTiles-function provides a quick access to some free third-party basemaps that are encapsulated in the “providers”-list, a data object in the leaflet-package.

Leaflet build-in alternatives

Some potential ocean background could be (for a full sneak-preview see here):

m <- leaflet() %>% setView(lng = -59, lat = 13, zoom = 5)

m %>% addProviderTiles(providers$Esri.WorldImagery)

m %>% addProviderTiles(providers$Esri.OceanBasemap)

Alternatives

Some recent maps tiles of potential use as backgrounds for marine data can be found on NOAA map homepage. Example for calling these tiles are:

# Create functions - just a wrapper around a long url
addTiles_noaa2 <- function(map) {
  map %>%
    addTiles(url = "https://tiles.arcgis.com/tiles/C8EMgrsFcRFL6LrL/arcgis/rest/services/GEBCO_2019_basemap_ncei/MapServer/tile/{z}/{y}/{x}")
}
addTiles_noaa3 <- function(map) {
  map %>%
    addTiles(url = "https://tiles.arcgis.com/tiles/C8EMgrsFcRFL6LrL/arcgis/rest/services/web_mercator_gebco_2019_contours/MapServer/tile/{z}/{y}/{x}")
}
m <- leaflet() %>% setView(-59, lat = 13, zoom = 5)
m %>% 
  addTiles_noaa2() %>% 
  addTiles_noaa3() |> 
  addPolygons(data = 
                 ne_countries(scale = 50, returnclass = "sf") |> 
                 st_cast("POLYGON"),
              label = ~admin,
              color = "red")

Own tiles - and layer control

If one has access to a server one could even attempt to provide ones own tiles, like this bottom topography shaded relief map:

m <- 
  m %>% 
  addProviderTiles(providers$Esri.WorldImagery,
                   group = "World Image") %>% 
  # Multiple calls to disparate tile sets
  addTiles(urlTemplate = "http://www.hafro.is/~einarhj/tiles2/olex-rayshaded/{z}/{x}/{y}.png",
           options = tileOptions(minZoom = 5, maxZoom = 12),
           group = "Crowd sourced (Olex)") %>%
  addTiles(urlTemplate = "http://www.hafro.is/~einarhj/tiles2/mb-rayshaded/{z}/{x}/{y}.png",
           options = tileOptions(minZoom = 5, maxZoom = 12),
           group = "MRI multibeam") %>%
  addTiles(urlTemplate = "http://www.hafro.is/~einarhj/tiles2/mb2-rayshaded/{z}/{x}/{y}.png",
           options = tileOptions(minZoom = 5, maxZoom = 12),
           group = "MRI multibeam")

The above script includes some options, here the minimum and the maximum zoom of 5 and 12 specifies the zoom level that the tiles are made visible to the user. The argument group comes in handy when allowing the user to control what layers are shown/hidden. For that we use the addLayersControl-function:

m %>% 
  addLayersControl(overlayGroups = c("World Image", "Crowd sourced (Olex)", "MRI multibeam"),
                   options = layersControlOptions(collapsed = FALSE)) %>% 
  setView(lng = -19, lat = 65, zoom = 5)

Polygons

iceland <- 
  ne_countries(country = "Iceland",
               scale = 10,
               returnclass = "sf")
class(iceland)

## [1] "sf"         "data.frame"

leaflet() %>% 
  addPolygons(data = iceland,
              color = "grey90",
              weight = 0) %>% 
  addCircles(data = w)

Note that in the above no call was made to background tile. This may often be the best option, because colourful background may create a distraction from the main data that supposedly contains the message for the reader.

Lines

Lets take the cruise track from the Icelandic spring survey in 2019 as an example.

trail <- 
  read.csv("ftp://ftp.hafro.is/pub/data/csv/is_smb_vms2019.csv",
           stringsAsFactors = FALSE) %>% 
  as_tibble() %>% 
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326) %>% 
  group_by(vessel) %>% 
  summarise(do_union = FALSE) %>% 
  st_cast("LINESTRING")

pal <-
  colorFactor(palette = RColorBrewer::brewer.pal(n = 5, name = "Set1"),
              domain = trail$vessel)
trail %>% 
  leaflet() %>% 
  addTiles_noaa3() |> 
  addPolylines(#data = trail,
    color = ~pal(vessel),
    group = ~vessel,
    opacity = 1,
    weight = 3) %>% 
  addLegend("bottomleft",
            title  = 'Vessels',
            pal = pal,
            values = ~vessel) %>% 
  addLayersControl(overlayGroups = c(trail$vessel),
                   options = layersControlOptions(collapsed = FALSE))

Line-segments

Lets use all the tows in the DATRAS database from 2018 as and example:

hh <- 
  read.csv("ftp://ftp.hafro.is/pub/data/csv/datras_2018_haul.csv",
           stringsAsFactors = FALSE) %>% 
  as_tibble() 
hh %>% 
  mutate(id = paste(id = paste0(survey, "_", id))) %>% 
  select(id,
         lon_start = shootlong,
         lat_start = shootlat,
         lon_end = haullong,
         lat_end = haullat) %>% 
  pivot_longer(-id, 
               names_to = c(".value", "action"), 
               names_sep = "_") %>% 
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326) %>% 
  group_by(id) %>% 
  summarise(do_union = FALSE) %>% 
  st_cast("LINESTRING") %>% 
  leaflet() %>% 
  addTiles() %>% 
  addPolylines(weight = 3,
               opacity = 1,
               label = ~id)

Exercise

Use the data of your liking and experiment with various basemaps of choice. Give some thought on:

• The use of colours for the data, given the different backgrounds chosen.
• What backgrounds are useful vs those that act as a distraction from the message you are trying to convey.

Multiple layers

Example is the abundance mackerel from the DATRAS 2018 survey plotted over ecoregions.

egos <- 
  read_sf("ftp://ftp.hafro.is/pub/data/shapes/ices_ecoregions.gpkg")
hl <- 
  read.csv("ftp://ftp.hafro.is/pub/data/csv/datras_2018_length.csv",
           stringsAsFactors = FALSE) %>% 
  as_tibble()

# some common species
mysp <- c("Scomber scombrus")
hh %>% 
  select(id, lon = shootlong, lat = shootlat) %>% 
  left_join(hl %>% filter(latin == mysp)) %>% 
  group_by(id, lon, lat) %>% 
  summarise(n = sum(n, na.rm = TRUE)) %>% 
  mutate(nr = ifelse(n > 1000, 1000, n)) %>% 
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326) %>% 
  leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = egos,
              color = ~colorFactor("viridis", ecoregion)(ecoregion),
              weight = 1, 
              group = ~ecoregion,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.3) %>% 
  addCircles(radius = 1000,
             color = "blue",
             weight = 0,
             fillOpacity = 1) %>% 
  addCircles(radius = ~nr * 100,
             label = ~as.character(round(n)),
             weight = 0,
             color = "red") %>% 
  addLayersControl(overlayGroups = c(egos$ecoregion))

Raster image

nephrops  <- raster::raster("ftp://ftp.hafro.is/pub/data/rasters/nephrops.tif")
inf <- viridis::inferno(12, alpha = 1, begin = 0, end = 1, direction = -1)
pal <- leaflet::colorNumeric(inf, raster::values(nephrops), na.color = "transparent")
m <- 
  leaflet() %>%
  setView(lng = -15, lat = 64, zoom = 7) %>%
  addProviderTiles(providers$Esri.WorldImagery) %>% 
  addRasterImage(nephrops, colors = pal, opacity = 1, group = "Nephrops trawl",
                 maxBytes = Inf)
m

The real thing

Saving as html

library(htmltools)
saveWidget(m, "nephrops.html")

• The html-document can be put on your institutes webserver for others to view
• Or you could send is like any other document as an attachment

Other packages

E.g.:

• {mapview}
  • Based on {leaflet} with helpful defaults.
• {mapedit}
• …