Section 20 Maps with `leaflet`

For spatial data corresponding to geographical locations, it can be helpful to visualise the data on a map. There are various mapping tools available in R. Here, we will discuss the R package leaflet (Cheng, Karambelkar, and Xie 2019), which is an R implementation of the JavaScript library of the same name.

Maps produced using leaflet are best suited for web pages, as you can zoom and scroll a map as you would do with, say, Google maps. These maps can also be used inside shiny apps.

You will need to be online to use leaflet.

20.1 A basic map

We’ll first load the leaflet package.

library(leaflet)

To get a basic map, we will normally need to specify the longitude and latitude of the centre of our map, as well the initial zoom level (which may take some trial and error). For example, to obtain a map of Sheffield:

leaflet() %>%
  addTiles() %>%
  setView(lng = -1.473798, 
          lat = 53.38252, 
          zoom = 13)

The addTiles() specifies what is drawn on each map ‘tile’, with the default being streets and points of interest provided by OpenStreetMap. For example, to change to an aerial photo (as in Google Earth), we can instead use the addProviderTiles() command, and specify a different ‘tile provider’:

leaflet() %>%
  setView(lng = -1.473798, lat = 53.38252, zoom = 13) %>%
  addProviderTiles(providers$Esri.WorldImagery)

20.2 Indicating points on a map

If we have vectors of coordinates (latitude and longitude), we can add points with the addMarkers() function (there are variants such as addCircles() and addPopups().) We can include the argument popup so that some text is displayed if the mouse is clicked on the marker. Here’s an example where I’ve marked where I work, and my favourite pub.

latitude <- c(53.380909, 53.371566)
longitude <- c(-1.486197, -1.577604)
placeNames <- c("Hicks Building", "The Sportsman")
leaflet() %>%
  setView(lng = -1.525, lat = 53.38, zoom =12) %>%
  addTiles() %>%
  addMarkers(lat = latitude,
             lng = longitude,
             popup = placeNames)

Exercise 20.1 Use leaflet to make a map that displays a region and two points of interest to you, for example, where you live, and your nearest railway station.

20.3 Plotting region boundaries

We can plot region boundaries on maps, and use different colours to make the regions more distinctive, or to indicate different values of some variable in different regions (a “choropleth map”). This can be a little tricky, in that the data sets that indicate region boundaries can be quite complicated.

Here’s an example of plotting ‘county and unitary’ authorities in the UK. To start with, we’ll use a base map that uses plainer-looking tiles than OpenStreetMap:

leaflet() %>%
  setView(lng = -3.436, lat = 55.3781, zoom = 5) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas)

20.3.1 shapefiles

Boundary data can be stored in different formats. We will search for a format called a shapefile. If we try a Google search for “england county boundaries shapefile”, we can find this page at data.gov.uk, and a zip file containing the shapefile can be downloaded from here.

20.3.2 Importing and simplifying shapefiles

Shapefiles can be imported using the sf package (Pebesma and Bivand, 2023; Pebesma 2018) with the functions read_sf() and st_transform(). The read_sf() command takes the full file path, including the filename (with the .shp extension), of the shapefile as its argument. We must then use the st_transform() function to convert the obtained sf object so that it has the appropriate ‘coordinate reference system (crs)’ for map plotting with the leaflet package. We’ll suppose here that the path to the folder where the shapefile is stored (after being unzipped from the download) is "U:/MAS61004/EDA/Counties/": we do:

library(sf)
boundaries <- read_sf('U:/MAS61004/EDA/Counties/CTYUA_MAY_2023_UK_BGC.shp') %>%
                  st_transform(crs = '+proj=longlat +datum=WGS84')

This creates an object of class sf: it’s more complicated than an ordinary dataframe, and we won’t worry too much about its structure. But if we try

print(boundaries)

## Simple feature collection with 218 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -8.650007 ymin: 49.8648 xmax: 1.76368 ymax: 60.86074
## Geodetic CRS:  +proj=longlat +datum=WGS84
## # A tibble: 218 × 9
##    CTYUA23CD CTYUA23NM            CTYUA23NMW  BNG_E  BNG_N   LONG   LAT GlobalID
##  * <chr>     <chr>                <chr>       <dbl>  <dbl>  <dbl> <dbl> <chr>   
##  1 E06000001 Hartlepool           <NA>       447160 531474 -1.27   54.7 3bf4312…
##  2 E06000002 Middlesbrough        <NA>       451141 516887 -1.21   54.5 f97de8d…
##  3 E06000003 Redcar and Cleveland <NA>       464361 519597 -1.01   54.6 ef76d83…
##  4 E06000004 Stockton-on-Tees     <NA>       444940 518183 -1.31   54.6 05fc23f…
##  5 E06000005 Darlington           <NA>       428029 515648 -1.57   54.5 5a21326…
##  6 E06000006 Halton               <NA>       354246 382146 -2.69   53.3 f87a98b…
##  7 E06000007 Warrington           <NA>       362744 388456 -2.56   53.4 51f4cdf…
##  8 E06000008 Blackburn with Darw… <NA>       369490 422806 -2.46   53.7 11e5f05…
##  9 E06000009 Blackpool            <NA>       332819 436635 -3.02   53.8 8ec1c78…
## 10 E06000010 Kingston upon Hull,… <NA>       511894 431650 -0.304  53.8 3fef2d1…
## # ℹ 208 more rows
## # ℹ 1 more variable: geometry <MULTIPOLYGON [°]>

we can see we have 218 regions (which sounds right!).

The shapefile is a large object; it may help to simplify it, so that our code will run a little more quickly. We can do this with the rmapshaper package (Teucher and Russell 2020).

simplifiedBoundaries <- rmapshaper::ms_simplify(boundaries)

(It’s worth checking that we get the same-looking map, whether we use boundaries or simplifiedBoundaries.)

20.3.3 Plotting the boundaries

We can now add the boundaries to our map, using the command addPolygons():

leaflet() %>%
  setView(lng = -3.436, lat = 55.3781, zoom = 5) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  addPolygons(data = simplifiedBoundaries,
              color = "blue",
              fillOpacity = 0,
              weight = 1,
              popup = ~CTYUA23NM)

For the arguments to addPolygons, we include

fillOpacity: a number between 0 and 1. Set to 0 for no colouring within each region;
weight: the thickness of the boundary lines;
popup = ~CTYUA23NM: if the mouse is clicked inside a region, display the value in the CTYUA23NM column in the simplifiedBoundaries data frame.

20.3.4 Using different colours for regions

We can use a different colour for each region, depending on some other data value. Here, we will colour the counties by the UK country they are in, where England, Scotland, Wales and Northern Ireland counties are coloured orange, blue, red and green respectively. If we inspect simplifiedBoundaries$CTYUA23NM, we can see the order of the counties and we can create a colour vector and assign the appropriate colours for the regions manually:

country = rep("orange", times = length(simplifiedBoundaries$CTYUA23NM))
country[154:164] = "green"; country[165:196] = "blue"; country[197:218] = "red"

We’ll redraw the map, increasing fillOpacity slightly above 0:

leaflet() %>%
  setView(lng = -3.436, lat = 55.3781, zoom = 5) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  addPolygons(data = simplifiedBoundaries,
              fillOpacity = 0.05,
              weight = 1,
              popup = ~CTYUA23NM,
              fillColor = country,
              color = country)

20.4 Further reading

The leaflet R package is well documented: a manual is available here.
Although not necessary for using the R package, if you want to understand more about leaflet itself, documentation and tutorials are available at The JavaScript leaflet homepage.

20.5 Acknowledgement

The shapefile was provided by the Office for National Statistics and obtained from https://www.data.gov.uk/dataset/85228aec-fe0e-49bf-9455-df000d61e731/counties-and-unitary-authorities-may-2023-boundaries-uk-bgc, accessed 27/08/2024. Contains public sector information licensed under the Open Government Licence v3.0.

20.6 References:

Pebesma, E., & Bivand, R. (2023). Spatial Data Science: With Applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016

Pebesma, E., 2018. Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal 10 (1), 439-446, https://doi.org/10.32614/RJ-2018-009

Teucher, Andy, and Kenton Russell. 2020. Rmapshaper: Client for ’Mapshaper’ for ’Geospatial’ Operations. https://CRAN.R-project.org/package=rmapshaper.

References

Cheng, Joe, Bhaskar Karambelkar, and Yihui Xie. 2019. Leaflet: Create Interactive Web Maps with the JavaScript ’Leaflet’ Library. https://CRAN.R-project.org/package=leaflet.

Section 20 Maps with leaflet