Interactive Maps with Leaflet

Objectives

  1. Create interactive web maps using Leaflet and the leaflet R package
  2. Change base maps
  3. Symbolize and display point and polygon data
  4. Configure pop-ups
  5. Add layer lists and legends
  6. Render a web map to an HTML webpage

Overview

Leaflet is a free and open-source JavaScript library and Application Programming Interface (API) for developing interactive web maps. It is very powerful, and offers a free alternative to proprietary web mapping software, such as ArcGIS Online and the ArcGIS API for JavaScript. It is used by many organizations including NPR, The Washington Post, Facebook, GitHub, USA Today, and the European Commission. Here is a link to the official Leaflet webpage.

Since Leaflet is meant to be used in the web development environment, it makes use of JavaScript, the dominant client-side programming language on the web. So, to use it effectively you need to have some understanding of JavaScript. It also doesn’t hurt to have knowledge of Hypertext Markup Language (HTML) and Cascading Style Sheets (CSS). So, there is a bit of a learning curve. Fortunately, the leaflet package allows for the use of Leaflet without the need to learn JavaScript. Interactive web maps can be created using R code then rendered to an HTML webpage. Since this course focuses on R, I will demonstrate using Leaflet through the associated R package. If you have an interest in web development, I would recommend learning some HTML, CSS, and JavaScript; however, that is outside the scope of this course. Our Client-Side Web GIS course covers web development with HTML, CSS, and JavaScript along with the ArcGIS API for JavaScript and the Leaflet JavaScript API.

We will use a variety of packages in this module. If you don’t have these packages installed yet, you will need to do so before proceeding.

library(sf)
library(leaflet)
library(leaflet.extras)
library(leaflet.esri)
library(dplyr)
library(RColorBrewer)
library(htmlwidgets)
library(htmltools)

If you would like to install a development version of the leaflet package, you can do so using GitHub as demonstrated in the next code block, which is commented out.

#if (!require('devtools')) install.packages('devtools')
#devtools::install_github('rstudio/leaflet')

In this demonstration, I will produce an interactive web map to visualize data summarized by county for the high plains states in the United States. We will work with the high_plains_data.csv file. The elevation (“elev”), temperature (“temp”), and precipitation (“precip2”) data were extracted from raster grids provided by the PRISM Climate Group at Oregon State University. The elevation data are provided in feet, the temperature data represent 30-year annual normal mean temperature in degrees Celsius, and the precipitation data represent 30-year annual normal precipitation in inches. The original raster grids have a resolution of 4-by-4 kilometers and can be obtained here. I also summarized percent forest by county (“per_for”) from the 2011 National Land Cover Database (NLCD). NLCD data can be obtained from the Multi-Resolution Land Characteristics Consortium (MRLC). The link at the bottom of the page provides the example data and R Markdown file used to generate this module.

I begin by reading in the spatial data using the sf package.

#You must set your own working directory.
cities <- st_read("leaflet/city_points.shp")
precip <- st_read("leaflet/counties_precip.shp")
temp <- st_read("leaflet/county_temp.shp")
per_for <- st_read("leaflet/county_per_forest.shp")
states <- st_read("leaflet/hp_states_wm.shp")

Generally, web maps are generated using the WGS84 Web Mercator projection. So, data that will be used in a web map should be transformed to that projection. Or, they could be referenced to the WGS84 geographic datum. I have chosen to use the WGS84 Web Mercator projection, as defined by European Petroleum Survey Group code 4326. st_tranform() is used to make the transformation. I then change the column names to make them more interpretable.

cities_wm <- st_transform(cities, crs=4326)
per_for_wm <- st_transform(per_for, crs=4326)
precip_wm <- st_transform(precip, crs=4326)
temp_wm <- st_transform(temp, crs=4326)
states_wm <- st_transform(states, crs=4326)

names(per_for_wm) <- c("county", "per_for", "geometry")
names(precip_wm) <- c("county", "precip", "geometry")
names(temp_wm) <- c("county", "temp", "geometry")

The cities point data contains an ordinal field representing different population range categories as follows:

  • 5: 2,500 to 9,999,
  • 6: 10,000 to 49,999
  • 7: 50,000 to 99,999
  • 8: 100,000 to 499,999
  • 9: 500,000 to 999,999

In order to display the categories separately, I split the data into new objects using the filter() function from dplyr after converting the data to a factor and printing the levels.

cities_wm$POP_CLASS <- as.factor(cities_wm$POP_CLASS)
levels(cities_wm$POP_CLASS)
[1] "5" "6" "7" "8" "9"
cities_5 <- cities_wm %>% filter(POP_CLASS==5)
cities_6 <- cities_wm %>% filter(POP_CLASS==6)
cities_7 <- cities_wm %>% filter(POP_CLASS==7)
cities_8 <- cities_wm %>% filter(POP_CLASS==8)
cities_9 <- cities_wm %>% filter(POP_CLASS==9)

Base Maps and Extents

Generally, you start designing a web map by defining a base map and default extent and zoom level. In the first example, I am simply creating a map with the default base map, extent, and zoom level by initializing the map with the leaflet() function then adding the base map using addTiles(). Raster Tile Layers represent pictures or images of raw data as opposed to the original geospatial data files. They are stored as tiles that change based on the zoom level. Most base maps are stored as Raster Tile Layers.

leaflet() %>%
addTiles()

In the example above, I did not explicitly state the desired extent and zoom level, and by default the map loaded to a global extent. The setView() function is used to define the extent using latitude and longitude coordinates. The zoom level is then defined. Larger numbers indicate more zoomed in. A zoom level of 1 indicates a global extent. Take some time to change the coordinates and zoom level to see how this impacts the map. Note that positive values indicate north latitude or east longitude while negative values indicate south latitude or west longitude.

leaflet() %>%
addTiles() %>%
setView(lat= 37.1, lng=-95.7, zoom=4)

It is also possible to limit the user’s ability to zoom and pan using leafletOptions(). In the next code block I am defining a default extent and a default zoom level. Using minZoom and maxZoom I limit the allowed zoom levels. Setting dragging to FALSE will not allow the user to pan the map.

We will not use these options to produce our maps; however, I wanted to demo them, as this can be useful in many situations.

leaflet(options=leafletOptions(center=c(54.5, 15.2), zoom =3, minZoom=2, maxZoom=5, dragging=FALSE)) %>%
addTiles()

There are many base maps available other than the default OpenStreetMap base map. To print a list of available base maps, you can use the code shown below. I have had issues in the past getting all of these layers to work depending on the computer and/or browser I am using. Unfortunately, I have not found a workaround for this issue. So, if you cannot get a base map to work, simply use a different one.

names(providers)
  [1] "OpenStreetMap"                      
  [2] "OpenStreetMap.Mapnik"               
  [3] "OpenStreetMap.DE"                   
  [4] "OpenStreetMap.CH"                   
  [5] "OpenStreetMap.France"               
  [6] "OpenStreetMap.HOT"                  
  [7] "OpenStreetMap.BZH"                  
  [8] "OpenSeaMap"                         
  [9] "OpenPtMap"                          
 [10] "OpenTopoMap"                        
 [11] "OpenRailwayMap"                     
 [12] "OpenFireMap"                        
 [13] "SafeCast"                           
 [14] "Thunderforest"                      
 [15] "Thunderforest.OpenCycleMap"         
 [16] "Thunderforest.Transport"            
 [17] "Thunderforest.TransportDark"        
 [18] "Thunderforest.SpinalMap"            
 [19] "Thunderforest.Landscape"            
 [20] "Thunderforest.Outdoors"             
 [21] "Thunderforest.Pioneer"              
 [22] "Thunderforest.MobileAtlas"          
 [23] "Thunderforest.Neighbourhood"        
 [24] "OpenMapSurfer"                      
 [25] "OpenMapSurfer.Roads"                
 [26] "OpenMapSurfer.Hybrid"               
 [27] "OpenMapSurfer.AdminBounds"          
 [28] "OpenMapSurfer.ContourLines"         
 [29] "OpenMapSurfer.Hillshade"            
 [30] "OpenMapSurfer.ElementsAtRisk"       
 [31] "Hydda"                              
 [32] "Hydda.Full"                         
 [33] "Hydda.Base"                         
 [34] "Hydda.RoadsAndLabels"               
 [35] "MapBox"                             
 [36] "Stamen"                             
 [37] "Stamen.Toner"                       
 [38] "Stamen.TonerBackground"             
 [39] "Stamen.TonerHybrid"                 
 [40] "Stamen.TonerLines"                  
 [41] "Stamen.TonerLabels"                 
 [42] "Stamen.TonerLite"                   
 [43] "Stamen.Watercolor"                  
 [44] "Stamen.Terrain"                     
 [45] "Stamen.TerrainBackground"           
 [46] "Stamen.TerrainLabels"               
 [47] "Stamen.TopOSMRelief"                
 [48] "Stamen.TopOSMFeatures"              
 [49] "TomTom"                             
 [50] "TomTom.Basic"                       
 [51] "TomTom.Hybrid"                      
 [52] "TomTom.Labels"                      
 [53] "Esri"                               
 [54] "Esri.WorldStreetMap"                
 [55] "Esri.DeLorme"                       
 [56] "Esri.WorldTopoMap"                  
 [57] "Esri.WorldImagery"                  
 [58] "Esri.WorldTerrain"                  
 [59] "Esri.WorldShadedRelief"             
 [60] "Esri.WorldPhysical"                 
 [61] "Esri.OceanBasemap"                  
 [62] "Esri.NatGeoWorldMap"                
 [63] "Esri.WorldGrayCanvas"               
 [64] "OpenWeatherMap"                     
 [65] "OpenWeatherMap.Clouds"              
 [66] "OpenWeatherMap.CloudsClassic"       
 [67] "OpenWeatherMap.Precipitation"       
 [68] "OpenWeatherMap.PrecipitationClassic"
 [69] "OpenWeatherMap.Rain"                
 [70] "OpenWeatherMap.RainClassic"         
 [71] "OpenWeatherMap.Pressure"            
 [72] "OpenWeatherMap.PressureContour"     
 [73] "OpenWeatherMap.Wind"                
 [74] "OpenWeatherMap.Temperature"         
 [75] "OpenWeatherMap.Snow"                
 [76] "HERE"                               
 [77] "HERE.normalDay"                     
 [78] "HERE.normalDayCustom"               
 [79] "HERE.normalDayGrey"                 
 [80] "HERE.normalDayMobile"               
 [81] "HERE.normalDayGreyMobile"           
 [82] "HERE.normalDayTransit"              
 [83] "HERE.normalDayTransitMobile"        
 [84] "HERE.normalDayTraffic"              
 [85] "HERE.normalNight"                   
 [86] "HERE.normalNightMobile"             
 [87] "HERE.normalNightGrey"               
 [88] "HERE.normalNightGreyMobile"         
 [89] "HERE.normalNightTransit"            
 [90] "HERE.normalNightTransitMobile"      
 [91] "HERE.reducedDay"                    
 [92] "HERE.reducedNight"                  
 [93] "HERE.basicMap"                      
 [94] "HERE.mapLabels"                     
 [95] "HERE.trafficFlow"                   
 [96] "HERE.carnavDayGrey"                 
 [97] "HERE.hybridDay"                     
 [98] "HERE.hybridDayMobile"               
 [99] "HERE.hybridDayTransit"              
[100] "HERE.hybridDayGrey"                 
[101] "HERE.hybridDayTraffic"              
[102] "HERE.pedestrianDay"                 
[103] "HERE.pedestrianNight"               
[104] "HERE.satelliteDay"                  
[105] "HERE.terrainDay"                    
[106] "HERE.terrainDayMobile"              
[107] "FreeMapSK"                          
[108] "MtbMap"                             
[109] "CartoDB"                            
[110] "CartoDB.Positron"                   
[111] "CartoDB.PositronNoLabels"           
[112] "CartoDB.PositronOnlyLabels"         
[113] "CartoDB.DarkMatter"                 
[114] "CartoDB.DarkMatterNoLabels"         
[115] "CartoDB.DarkMatterOnlyLabels"       
[116] "CartoDB.Voyager"                    
[117] "CartoDB.VoyagerNoLabels"            
[118] "CartoDB.VoyagerOnlyLabels"          
[119] "CartoDB.VoyagerLabelsUnder"         
[120] "HikeBike"                           
[121] "HikeBike.HikeBike"                  
[122] "HikeBike.HillShading"               
[123] "BasemapAT"                          
[124] "BasemapAT.basemap"                  
[125] "BasemapAT.grau"                     
[126] "BasemapAT.overlay"                  
[127] "BasemapAT.highdpi"                  
[128] "BasemapAT.orthofoto"                
[129] "nlmaps"                             
[130] "nlmaps.standaard"                   
[131] "nlmaps.pastel"                      
[132] "nlmaps.grijs"                       
[133] "nlmaps.luchtfoto"                   
[134] "NASAGIBS"                           
[135] "NASAGIBS.ModisTerraTrueColorCR"     
[136] "NASAGIBS.ModisTerraBands367CR"      
[137] "NASAGIBS.ViirsEarthAtNight2012"     
[138] "NASAGIBS.ModisTerraLSTDay"          
[139] "NASAGIBS.ModisTerraSnowCover"       
[140] "NASAGIBS.ModisTerraAOD"             
[141] "NASAGIBS.ModisTerraChlorophyll"     
[142] "NLS"                                
[143] "JusticeMap"                         
[144] "JusticeMap.income"                  
[145] "JusticeMap.americanIndian"          
[146] "JusticeMap.asian"                   
[147] "JusticeMap.black"                   
[148] "JusticeMap.hispanic"                
[149] "JusticeMap.multi"                   
[150] "JusticeMap.nonWhite"                
[151] "JusticeMap.white"                   
[152] "JusticeMap.plurality"               
[153] "Wikimedia"                          
[154] "GeoportailFrance"                   
[155] "GeoportailFrance.parcels"           
[156] "GeoportailFrance.ignMaps"           
[157] "GeoportailFrance.maps"              
[158] "GeoportailFrance.orthos"            
[159] "OneMapSG"                           
[160] "OneMapSG.Default"                   
[161] "OneMapSG.Night"                     
[162] "OneMapSG.Original"                  
[163] "OneMapSG.Grey"                      
[164] "OneMapSG.LandLot"

The next set of code blocks demonstrate loading in different base maps by providing the layer’s name within the addProviderTiles() function. Copyright and source information are automatically added to the interactive map.

Take some time to experiment with different base maps by altering the example code.

leaflet() %>%
addProviderTiles("Stamen.Watercolor") %>%
setView(lat= 54.5, lng=15.2, zoom=3)
leaflet() %>%
addProviderTiles("Esri.NatGeoWorldMap") %>%
setView(lat= 54.5, lng=15.25, zoom=3)

What if you would like the user to be able to switch between some pre-defined base maps? This can be accomplished by loading in multiple layers using addTiles() or addProviderTiles(). You will then need to use addLayersControl() to create a widget to switch between the base maps. In order to reference the layers in the addLayersControl() function you will need to assign a unique name to each one using group. These names will be used in the widget. A collapse argument can be added to the addLayersControl() function to indicate whether or not it should collapse.

leaflet() %>%
  
addTiles(group = "OSM") %>%
  
addProviderTiles("Esri.NatGeoWorldMap", group="ESRI") %>%
  
addProviderTiles("CartoDB.DarkMatter", group= "CartoDB") %>%
  
addLayersControl(baseGroups = c("OSM", "ESRI", "CartoDB")) %>%

setView(lat= 42, lng=-105, zoom=5)

Add Points

I now have a map with three base map layers that the user can choose between. I have also defined an initial extent and zoom level. However, I haven’t added any data yet.

To begin, I will add the smallest cities to the map using addMarkers(). Within the function, I must specify the data to be used. I am also defining the content of the pop-up using the popup argument. This will simply print “Some Text” in each pop-up.

Throughout the remaining examples, you will see the htmlEscape() function from the htmltools package used often. This function is used to make sure that pop-up content is not incorrectly interpreted as HTML tags or code. Also, it should be used to avoid any security risks. In short, you should use this function when creating pop-ups.

leaflet() %>%
  
addTiles(group = "OSM") %>%
  
addProviderTiles("Esri.NatGeoWorldMap", group="ESRI") %>%
  
addProviderTiles("CartoDB.DarkMatter", group= "CartoDB") %>%
  
addMarkers(data=cities_5, popup = ~htmlEscape(paste0("Some Text"))) %>%
  
addLayersControl(baseGroups = c("OSM", "ESRI", "CartoDB")) %>%
  
setView(lat= 42, lng=-105, zoom=5)

Other symbology options are available for point data. In this example I am using addCircleMarkers() as opposed to addMarkers(). I am also now displaying all the city population groups using multiple addCircleMarkers() calls. Using the radius and color arguments, I am specifying sizes and colors for each symbol. I am setting fillOpacity to 1 so that no transparency is applied and the stroke to FALSE so that no outline is included. I am also providing a group name for each point layer.

The pop-ups are more complicated than the above example. I am having the pop-up provide the city name and the city population. The “b” tags are HTML tags used to bold text while the “/br” tag indicates a line break. Note that these components are not wrapped in htmlEscape() since I do want them to be interpreted as HTML tags as opposed to text. I am using the base R paste0() function to combine all arguments into a single string. I am also indicating that I want comma separators included in the population numbers.

You may have noticed the use of ~ in the code. This allows for the data object to be passed on. So, I can just use the column name in the pop-up syntax and do not need to define the data object.

leaflet() %>%
  
addTiles(group = "OSM") %>%
  
addProviderTiles("Esri.NatGeoWorldMap", group="ESRI") %>%
  
addProviderTiles("CartoDB.DarkMatter", group= "CartoDB") %>%

addCircleMarkers(data=cities_5, group="2,500 to 9,999", radius=3, color="#ffffd4", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>", htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircleMarkers(data=cities_6, group="10,000 to 49,999", radius=5, color="#fed98e", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%

addCircleMarkers(data=cities_7, group="50,000 to 90,999", radius=7, color="#fe9929", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircleMarkers(data=cities_8, group="100,000 to 499,999", radius=9, color="#d95f0e", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>", htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircleMarkers(data=cities_9, group="500,000 to 999,999", radius=11, color="#993404", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addLayersControl(baseGroups = c("OSM", "ESRI", "CartoDB")) %>%
  
setView(lat= 42, lng=-105, zoom=5)

Similar to the base maps, you may also want the user to be able to turn the point layers on and off. This can be accomplished by adding an overlayGroups argument to addLayerConrols().

leaflet() %>%
  
addTiles(group = "OSM") %>%
  
addProviderTiles("Esri.NatGeoWorldMap", group="ESRI") %>%
  
addProviderTiles("CartoDB.DarkMatter", group= "CartoDB") %>%
  
addCircleMarkers(data=cities_5, group="2,500 to 9,999", radius=3, color="#ffffd4", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>", htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircleMarkers(data=cities_6, group="10,000 to 49,999", radius=5, color="#fed98e", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%

addCircleMarkers(data=cities_7, group="50,000 to 90,999", radius=7, color="#fe9929", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircles(data=cities_8, group="100,000 to 499,999", radius=1100, color="#d95f0e", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>", htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircleMarkers(data=cities_9, group="500,000 to 999,999", radius=11, color="#993404", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addLayersControl(overlayGroups= c("2,500 to 9,999", "10,000 to 49,999", "50,000 to 90,999", "100,000 to 499,999", "500,000 to 999,999"), baseGroups = c("OSM", "ESRI", "CartoDB")) %>%
  
setView(lat= 42, lng=-105, zoom=5)

Add Polygons

I will now add some polygon layers to the map, which can be accomplished using the addPolygons() function. In this first example, I am adding the county-level percent forest data. I am not defining the symbology, so the default is used, which isn’t great for this use case.

leaflet() %>%
  
addTiles(group = "OSM") %>%
  
addProviderTiles("Esri.NatGeoWorldMap", group="ESRI") %>%
  
addProviderTiles("CartoDB.DarkMatter", group= "CartoDB") %>%
  
addPolygons(data=per_for_wm) %>%
  
addCircleMarkers(data=cities_5, group="2,500 to 9,999", radius=3, color="#ffffd4", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>", htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircleMarkers(data=cities_6, group="10,000 to 49,999", radius=5, color="#fed98e", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%

addCircleMarkers(data=cities_7, group="50,000 to 90,999", radius=7, color="#fe9929", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircleMarkers(data=cities_8, group="100,000 to 499,999", radius=9, color="#d95f0e", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>", htmlEscape(format(POP2010, big.mark=",")))) %>%
  
addCircleMarkers(data=cities_9, group="500,000 to 999,999", radius=11, color="#993404", stroke=FALSE, fillOpacity=1, popup = ~paste0("<b>", htmlEscape(NAME), "</b>","<br/>",htmlEscape(format(POP2010, big.mark=",")))) %>%

addLayersControl(overlayGroups= c("2,500 to 9,999", "10,000 to 49,999", "50,000 to 90,999", "100,000 to 499,999", "500,000 to 999,999"), baseGroups = c("OSM", "ESRI", "CartoDB")) %>%
  
setView(lat= 42, lng=-105, zoom=5)