Making Interactive Maps in R

A Rapid Introduction to Interactive Mapping in R for Research Scientists: All/some things Leaflet.

One of the main issues with doing mapping in R over other software is our ability to zoom in and out of maps to assess different scales, from whole ocean basin to small harbour? Likewise, when planning fieldwork or even just assessing whether the coordinates you recorded are in the right location with no missing minus signs or wrong decimal location, having a tool you can quickly upload points onto a map that then allows you to zoom in and out on is really handy.

Leaflet

Leaflet maps are a Javascript (I think) mapping library that are used very often for just the purposes I highlighted above! Used as widgets on websites they allow uses to zoom and out of maps rather than having static boring maps. Thankfully there is a leaflet package in r. And it is SIMPLE to use.

#install.packages("leaflet")
#install.packages("leaflet.extras")
#install.packages("leaflet.esri")
#install.packages("leaflet.providers")

library(leaflet)
library(leaflet.extras)
library(leaflet.esri)
library(leaflet.providers)

leaflet()

Layering

Base Map

Like with ggplot2 the first blank arguement creates a blank layer, we can then layer preloaded and local data into the map. Unlike ggplot2 we use the pipe ( %>% ) rather than the plus (+).

leaflet() %>% 
  addTiles()

The base map is zoomed to 2 ish globes and with the standard Open Street Map. We can change these element by choosing some other background tile and setting the zoom. You can find a whole list of all the available basemaps here

leaflet() %>%
  addProviderTiles(
    "OpenTopoMap"
  )%>%
    setView(lng = 6.8652, lat = 45.8326, zoom = 6) ## Location of Mont Blanc

We can set one we like, or we can allow our user to choose. Here we will save the map as m but we also want to plot it too so we can put brackets around the whole thing, this means we can just run the line to save and plot.

(m<-leaflet()%>%
  addProviderTiles(
    "OpenTopoMap"
  ) )

Local Data

Point Data

So lets assume we are interested in the tallest mountains in each continent. Lets make a dataframe then add it to the leaflet map.

Mountains<-data.frame(
  Mountain=c("Everest", "Aconcagua", "Denali", "Kilimanjaro", "Vinson", "Mont Blanc", "Mount Wilhelm"),
  Longitude=c(86.9250,-70.0109,-151.0070, 37.3556,-85.2135, 6.8652,145.0297),
  Latitude=c(27.9881,-32.6532,63.0692,-3.0674,-78.6341,45.8326,-5.7800),
  Height=c(8848,6961,6194,5895,4892,4810,4509),
  Continent=c("Asia","South America","North America","Africa","Antarctica","Europe","Oceania")
)

m%>% 
  addMarkers(data=Mountains,lng=~Longitude,lat=~Latitude)

The base markers aren’t very nice so lets add some colour and style plus some info.

pal1 <- colorFactor(c("navy", "red", "green"),
                   domain = unique(Mountains$Continent))

m %>% 
  addCircleMarkers(data=Mountains,
                    lng=~Longitude,
                    lat=~Latitude,
                    popup=~paste0(Mountain," - ",Height),
                    label=~as.character(Mountain),
                   color=~pal1(Continent),
                   fillOpacity = 0.8
                   )

Spatial Data

Lets download and plot some bathymetry data for the Caribbean and Easter Tropical seas. We need to convert our bathy object to be a SpatRaster from the Terra package (also accepts raster objects but the raster package is being deprecated to switch to terra!!!!) Then we can add a colour palette along the bathymetry data we have. This is then just a couple arguments to make an alright legend.

library(marmap)

Registered S3 methods overwritten by 'adehabitatMA':
  method                       from
  print.SpatialPixelsDataFrame sp  
  print.SpatialPixels          sp  

Attaching package: 'marmap'

The following object is masked from 'package:grDevices':

    as.raster

library(terra)

terra 1.7.71

library(tidyterra)

Attaching package: 'tidyterra'

The following object is masked from 'package:stats':

    filter

bat_panama <- getNOAA.bathy(lon1=-100,lon2=-54,lat1=32,lat2=-10, res = 1)

Querying NOAA database ...

This may take seconds to minutes, depending on grid size

Building bathy matrix ...

Bathy_panama <- as.xyz(bat_panama) %>% 
  rename(Longitude=V1,Latitude=V2,Depth=V3) %>% 
  filter(Depth<0) %>% 
  rast()

crs(Bathy_panama)<-"+proj=longlat" 

pal2 <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(Bathy_panama),
  na.color = "transparent")

leaflet()%>%
  addTiles()%>%
  addRasterImage(Bathy_panama,colors=pal2,opacity=0.7) %>%
  addLegend(pal = pal2, values = values(Bathy_panama),
    title = "Depth (m)")