Web mapping — the easy way

First, let’s load our friends: sf (vector data), raster (raster data), and dplyr (general data tidying):

library(sf)
library(raster)
library(dplyr)

Then let’s load the mapview package, which makes the production of web maps in R very easy:

library(mapview)

Datasets

For this demonstration, we are using the cookfarm dataset from the GSIF package. This dataset is a list with 6 somponents storing various measurements that have been done on the farm.

data('cookfarm', package = "GSIF")

names(cookfarm)

## [1] "readings"    "profiles"    "bdensity"    "grids"       "weather"    
## [6] "proj4string"

Vector data

The profiles component stores soil profile data as a data.frame. We will turn this into a simple feature using the sf package:

profiles <- cookfarm$profiles %>% 
  st_as_sf(
    coords = c('Easting', 'Northing'), 
    crs = cookfarm$proj4string
  )

# simple plot using the sf package
plot(profiles)

Raster data

The grids component is storing gridded data. We will reansform these into a RasterStack using the rasterFromXYZ function. This function requires the x and y columns to be the first column of the data.frame (from the left) – a good opportunity to use the select function from the dplyr package:

grids <- cookfarm$grids %>% 
  select(x, y, DEM, TWI, Cook_fall_ECa, Cook_spr_ECa) %>%  # Re-order the columns, and select 4 interesting variables
  rasterFromXYZ(crs = cookfarm$proj4string) # The coordinate reference system is stored as one of the component in cookfarm

# simple plot using the raster package
plot(grids)

Interactive maps using mapview

Vector data

The mapview package is basically a wrapper around the leaflet package. It makes is very easy to create a leaflet map. You can basically just use it as you would use plot:

# couldn't be easier to create a map!
mapview(profiles)

But while the mapview function can be called just as is, there’s numerous details that can be tweaked. For example the colours:

# Create palettes 
library(RColorBrewer)
pal_continuous <- colorRampPalette(brewer.pal(7, "BrBG")) # For continuous data
pal_categorical <- colorRampPalette(brewer.pal(9, "Set1")) # For categorical data

# Pass the palette to mapview
mapview(
  profiles, 
  zcol = "TAXSUSDA", 
  col.regions = pal_categorical, 
  legend = TRUE
)

In particular, a wide range of background (web) maps is available. You can pick one or several. You can find a list of these background from the Leaflet-extras project website: http://leaflet-extras.github.io/leaflet-providers/preview/.

# Tweaking backgrounds
mapview(profiles, zcol = "BLD", col.regions = pal_continuous, legend = TRUE, map.types = "Esri.WorldImagery")

The burst option can be interesting when you visualise categorical data:

# Burst to separate soil classes
mapview(profiles, zcol = "TAXSUSDA", col.regions = pal_categorical, legend = TRUE, burst = TRUE)

Raster data

The mapview option also work for data loaded using the raster package. If you try to visualise a RasterStack (as opposed to RasterLayer), you can select and choose which layer to plot using the selection interface.

# Plot a RasterLayer
mapview(grids$TWI, col.regions = pal_continuous, na.color = "transparent", legend = TRUE)

# Plot a RasterStack
mapview(grids)

# The same sort of options are available
mapview(grids, col.regions = pal_continuous, na.color = "transparent",  legend = TRUE)

Extent of rasters

You can visualise the raster extent — rather than the data itself — using viewExtent:

viewExtent(grids)

Imagery data

For imagery, it is a little diffferent, because usually one wants to visualise a RGB composite of 3 layers, rather than the individual layers individually. In this case, the viewRGB is the way to go:

# Imagery specific functions
viewRGB(poppendorf)

# The combination of bands can be changed very easily
viewRGB(poppendorf, 4, 3, 2)

Adding mapviews together

Conveniently, maps can be aded to each other using the + operator:

# Create 2 maps
m1 <- mapview(grids$DEM, col.regions = pal_continuous, legend = TRUE) 
m2 <- mapview(profiles, zcol = "TAXSUSDA", col.regions = pal_categorical)

# Plot both together
m1 + m2

Visualising projected data

What mapview does behind the scenes is to change the projection system to EPSG:3857 (web mercator). In some cases, this is inconvenient and you might actually want to visualise your data in a local projected CRS. plainview is here to help:

plainview(grids$DEM)

Advanced topics

What happens behind the scenes

mapview is a wrapper around leaflet package, which is a R API for the popular Javascript library for web mapping called Leaflet. Leaflet has been designed with simplicity and rapidity in mind. For more power, you’ll have to learn a bit more about the leaflet package itself. Their website is a great starting point.

Changing default options

A bunch of options can be changed using the mapviewOptions function:

mapviewOptions(
  basemaps = c("Esri.WorldImagery", "Thunderforest.Landscape"),
  na.color = "transparent"
)

mapview(profiles, zcol = "BLD") + mapview(grids$DEM)

Sync’ing maps

You can associate and synchronise a set of maps using the sync function:

# Syncing several maps
m1 <- mapview(grids$DEM) 
m2 <- mapview(grids$TWI) 
m3 <- mapview(grids$Cook_fall_ECa)

sync(m1, m2, m3, ncol = 2, sync.cursor = TRUE)

This is an interactive analogue to the panelled graphs provided by ggplot2 or lattice.

Slideview

Another advanced visualisation tool is slideview, which is convenient to compare two maps:

img1 <- poppendorf[[1]]
img2 <- poppendorf[[5]]

slideview(
  img1, 
  img2,
  label1 = "Poppendorf-Layer-1",
  label2 = "Poppendorf-Layer-2",
  legend = TRUE
)

Popups

The popups can be either a table (popupTable, default behaviour), an image (popupImage), or a htmlwidget (popupGraph).

# Table
mapview(
  profiles, 
  popup = popupTable(profiles, zcol = 1:2)
)

# Image
mapview(
  profiles, 
  popup = popupImage('https://www.vcard.wur.nl/WebServices/GetMedia.ashx?id=37263')
)

Garnish map

The leaflet and leaflet.extras are providing a LOT of different map widgets. The garnishMaps function facilitates their integration with mapview:

library(leaflet)

m <- mapview(profiles)
garnishMap(
  m,
  addMouseCoordinates,
  addGraticule,
  addScaleBar
)

Save map

# Create map
m  <- mapview(profiles)

# Save interactive HTML page
mapshot(m, url = 'my_map.html')

# Save static image (PNG, JPEG, or PDF)
mapshot(m, file = 'my_image.png')

---

Super advanced mapview 🎉

library(xts)
library(dygraphs)

profiles$SOURCEID <- as.character(profiles$SOURCEID)

records <- cookfarm$readings
records$SOURCEID <- as.character(records$SOURCEID)

ids <- unique(records$SOURCEID)

# Subset sensors
ids <- sample(ids, size = 5) 
  
idx_sensors <- which(profiles$SOURCEID %in% ids)
sensors <- profiles[idx_sensors,]

make_ts <- function(id) {
  
  records %>% 
    filter(SOURCEID == id) %>% 
    dplyr::select(-SOURCEID) %>%
    dplyr::select(Date, ends_with('VW')) %>% 
    xts(.$Date)
}

make_dygraph <- function(id){
  ts <- make_ts(id)
  dygraph(ts)
} 

l_graphs <- lapply(
  ids,
  make_dygraph
)

make_dygraph(ids[1])

mapview(sensors, popup = popupGraph(graphs = l_graphs, width = 300, height = 300))