Lecture-05 Examples

Introduction

This notebook covers interactive map production in R using the leflet library. We’ll review some of what we discussed during the first course meeting and then introduce some additional complexity with leaflet.

Dependencies

This notebook requires a variety of packages for working with spatial data:

# tidyverse packages
library(dplyr)        # data wrangling

# spatial packages
library(leaflet)      # interactive maps
library(mapview)      # preview spatial data
library(sf)           # spatial data tools

# other packages
library(here)         # file path management
library(measurements) # convert units
library(RColorBrewer) # color palettes

Load Data

This notebook requires the data stored in data/example-data/. We’ll automatically convert these to WGS84 for mapping, since this is the coordinate system they require:

city <- st_read(here("data", "example-data", "STL_BOUNDARY_City", "STL_BOUNDARY_City.shp"), stringsAsFactors = FALSE) %>%
  st_transform(crs = 4326)

Reading layer `STL_BOUNDARY_City' from data source `/Users/chris/GitHub/slu-soc5650/content/module-5-leaflet/data/example-data/STL_BOUNDARY_City/STL_BOUNDARY_City.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 17 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -90.32052 ymin: 38.53185 xmax: -90.16657 ymax: 38.77443
Geodetic CRS:  GRS 1980(IUGG, 1980)

nhoods <- st_read(here("data", "example-data", "STL_DEMOS_Nhoods", "STL_DEMOS_Nhoods.shp"), stringsAsFactors = FALSE) %>%
  st_transform(crs = 4326)

Reading layer `STL_DEMOS_Nhoods' from data source `/Users/chris/GitHub/slu-soc5650/content/module-5-leaflet/data/example-data/STL_DEMOS_Nhoods/STL_DEMOS_Nhoods.shp' using driver `ESRI Shapefile'
Simple feature collection with 79 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 733361.8 ymin: 4268512 xmax: 745417.9 ymax: 4295501
Projected CRS: UTM_Zone_15_Northern_Hemisphere

Interactive Mapping with leaflet (Module 1 review)

Basic Mapping of Geometric Objects

During our first course meeting, we discussed the creation of interactive maps using the leaflet package. Leaflet is a JavaScript library for creating interactive maps. It is primarily focused on web and mobile mapping. The leaflet package for R provides access to the JavaScript library.

Leaflet provides a number of basemaps for mapping. If you add map tiles using addTiles(), you’ll get the open street map basemap. Other basemaps can be added using addProviderTiles(). The names of available options can be found using:

names(providers)

  [1] "OpenStreetMap"                       "OpenStreetMap.Mapnik"                "OpenStreetMap.DE"                   
  [4] "OpenStreetMap.CH"                    "OpenStreetMap.France"                "OpenStreetMap.HOT"                  
  [7] "OpenStreetMap.BZH"                   "OpenSeaMap"                          "OpenPtMap"                          
 [10] "OpenTopoMap"                         "OpenRailwayMap"                      "OpenFireMap"                        
 [13] "SafeCast"                            "Thunderforest"                       "Thunderforest.OpenCycleMap"         
 [16] "Thunderforest.Transport"             "Thunderforest.TransportDark"         "Thunderforest.SpinalMap"            
 [19] "Thunderforest.Landscape"             "Thunderforest.Outdoors"              "Thunderforest.Pioneer"              
 [22] "Thunderforest.MobileAtlas"           "Thunderforest.Neighbourhood"         "OpenMapSurfer"                      
 [25] "OpenMapSurfer.Roads"                 "OpenMapSurfer.Hybrid"                "OpenMapSurfer.AdminBounds"          
 [28] "OpenMapSurfer.ContourLines"          "OpenMapSurfer.Hillshade"             "OpenMapSurfer.ElementsAtRisk"       
 [31] "Hydda"                               "Hydda.Full"                          "Hydda.Base"                         
 [34] "Hydda.RoadsAndLabels"                "MapBox"                              "Stamen"                             
 [37] "Stamen.Toner"                        "Stamen.TonerBackground"              "Stamen.TonerHybrid"                 
 [40] "Stamen.TonerLines"                   "Stamen.TonerLabels"                  "Stamen.TonerLite"                   
 [43] "Stamen.Watercolor"                   "Stamen.Terrain"                      "Stamen.TerrainBackground"           
 [46] "Stamen.TerrainLabels"                "Stamen.TopOSMRelief"                 "Stamen.TopOSMFeatures"              
 [49] "TomTom"                              "TomTom.Basic"                        "TomTom.Hybrid"                      
 [52] "TomTom.Labels"                       "Esri"                                "Esri.WorldStreetMap"                
 [55] "Esri.DeLorme"                        "Esri.WorldTopoMap"                   "Esri.WorldImagery"                  
 [58] "Esri.WorldTerrain"                   "Esri.WorldShadedRelief"              "Esri.WorldPhysical"                 
 [61] "Esri.OceanBasemap"                   "Esri.NatGeoWorldMap"                 "Esri.WorldGrayCanvas"               
 [64] "OpenWeatherMap"                      "OpenWeatherMap.Clouds"               "OpenWeatherMap.CloudsClassic"       
 [67] "OpenWeatherMap.Precipitation"        "OpenWeatherMap.PrecipitationClassic" "OpenWeatherMap.Rain"                
 [70] "OpenWeatherMap.RainClassic"          "OpenWeatherMap.Pressure"             "OpenWeatherMap.PressureContour"     
 [73] "OpenWeatherMap.Wind"                 "OpenWeatherMap.Temperature"          "OpenWeatherMap.Snow"                
 [76] "HERE"                                "HERE.normalDay"                      "HERE.normalDayCustom"               
 [79] "HERE.normalDayGrey"                  "HERE.normalDayMobile"                "HERE.normalDayGreyMobile"           
 [82] "HERE.normalDayTransit"               "HERE.normalDayTransitMobile"         "HERE.normalDayTraffic"              
 [85] "HERE.normalNight"                    "HERE.normalNightMobile"              "HERE.normalNightGrey"               
 [88] "HERE.normalNightGreyMobile"          "HERE.normalNightTransit"             "HERE.normalNightTransitMobile"      
 [91] "HERE.reducedDay"                     "HERE.reducedNight"                   "HERE.basicMap"                      
 [94] "HERE.mapLabels"                      "HERE.trafficFlow"                    "HERE.carnavDayGrey"                 
 [97] "HERE.hybridDay"                      "HERE.hybridDayMobile"                "HERE.hybridDayTransit"              
[100] "HERE.hybridDayGrey"                  "HERE.hybridDayTraffic"               "HERE.pedestrianDay"                 
[103] "HERE.pedestrianNight"                "HERE.satelliteDay"                   "HERE.terrainDay"                    
[106] "HERE.terrainDayMobile"               "FreeMapSK"                           "MtbMap"                             
[109] "CartoDB"                             "CartoDB.Positron"                    "CartoDB.PositronNoLabels"           
[112] "CartoDB.PositronOnlyLabels"          "CartoDB.DarkMatter"                  "CartoDB.DarkMatterNoLabels"         
[115] "CartoDB.DarkMatterOnlyLabels"        "CartoDB.Voyager"                     "CartoDB.VoyagerNoLabels"            
[118] "CartoDB.VoyagerOnlyLabels"           "CartoDB.VoyagerLabelsUnder"          "HikeBike"                           
[121] "HikeBike.HikeBike"                   "HikeBike.HillShading"                "BasemapAT"                          
[124] "BasemapAT.basemap"                   "BasemapAT.grau"                      "BasemapAT.overlay"                  
[127] "BasemapAT.highdpi"                   "BasemapAT.orthofoto"                 "nlmaps"                             
[130] "nlmaps.standaard"                    "nlmaps.pastel"                       "nlmaps.grijs"                       
[133] "nlmaps.luchtfoto"                    "NASAGIBS"                            "NASAGIBS.ModisTerraTrueColorCR"     
[136] "NASAGIBS.ModisTerraBands367CR"       "NASAGIBS.ViirsEarthAtNight2012"      "NASAGIBS.ModisTerraLSTDay"          
[139] "NASAGIBS.ModisTerraSnowCover"        "NASAGIBS.ModisTerraAOD"              "NASAGIBS.ModisTerraChlorophyll"     
[142] "NLS"                                 "JusticeMap"                          "JusticeMap.income"                  
[145] "JusticeMap.americanIndian"           "JusticeMap.asian"                    "JusticeMap.black"                   
[148] "JusticeMap.hispanic"                 "JusticeMap.multi"                    "JusticeMap.nonWhite"                
[151] "JusticeMap.white"                    "JusticeMap.plurality"                "Wikimedia"                          
[154] "GeoportailFrance"                    "GeoportailFrance.parcels"            "GeoportailFrance.ignMaps"           
[157] "GeoportailFrance.maps"               "GeoportailFrance.orthos"             "OneMapSG"                           
[160] "OneMapSG.Default"                    "OneMapSG.Night"                      "OneMapSG.Original"                  
[163] "OneMapSG.Grey"                       "OneMapSG.LandLot"                   

As you can see, there are a ton of choices! We’ll use CartoDB.Positron here, but feel free to pick one that you like for assignments if we don’t specify what you should use. Make sure that your other cartographic selections, such as color, do not clash with your basemap.

The basic leaflet workflow involves piping functions together (the %>% operator). Each time to see the pipe, think of the word “then.” For example, the following code chunk would read:

1.Take the city object, then 2. use it as the basis for creating a leaflet object with leaflet(), then 3. add a basemap using the CartoDB.Positron tiles, then 4. add polygons and create a pop-up.

city %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(popup = ~NAME)

Adding Pop-ups

We can get more into the weeds with the neighborhood data since they have additional features. We can create more detailed pop-ups using the base::paste() function and some html tags. The most important html tags to know are:

• <b>text</b> - bold text
• <em>text</em> - italicized text
• <br> - line break

nhoods %>%  
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(popup = paste("<b>Name:</b> ", nhoods$NHD_NAME, "<br>",
                            "<b>2017 Population:</b> ", round(nhoods$pop17, digits = 0)))

Mapping Quantities with leaflet

If we want to turn this into a thematic choropleth map, we can add some additional parameters to the addPolygons() function. Before we do that, however, we should create normalized versions of our two population variables, pop50 and pop17. We can use the AREA column, which represents the area of each neighborhood in square meters, as a basis for this.

nhoods %>%
  mutate(sq_km = conv_unit(AREA, from = "m2", to = "km2"), .after = "AREA") %>%
  mutate(pop50_den = pop50/sq_km, .after = "pop50") %>%
  mutate(pop17_den = pop17/sq_km, .after = "pop17") -> nhoods

Now that we have these normalized, we can get mapping! The additional cartographic options we’ll mention are:

• color - outline (“stroke”) color for each polygon
• weight - stroke width
• opacity - stroke opacity
• smoothFactor - allows leaflet to simplify polygons depending on zoom
• fillOpacity - fill opacity
• fillColor - creates the fill itself
• highlightOptions - creates effect when mouse drags over specific polygons

What I have here are good default settings for most of these options, but feel free to experiment!

When we created our pop-up, we want to round our values so that we don’t see the very long real number associated with our data. By using base::round(var, digits = 0), we round to the nearest integer. digits = 2 would give us two decimal places in contrast.

# create color palette
npal <- colorNumeric("YlOrRd", nhoods$pop17_den)

# create leaflet object
nhoods %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    color = "#444444", 
    weight = 1, 
    opacity = 1.0, 
    smoothFactor = 0.5,
    fillOpacity = 0.5,
    fillColor = ~npal(pop17_den),
    highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE),
    popup = paste("<b>Name:</b> ", nhoods$NHD_NAME, "<br>",
                  "<b>2017 Population:</b> ", round(nhoods$pop17, digits = 0), "<br>",
                  "<b>2017 Population per Square Kilometer:</b> ", round(nhoods$pop17_den, digits = 2))) 

Next, we should add a legend to make the map easier to interpret. This is done with the addLegend() argument. The opacity argument in addLegend() should match the fillOpacity argument in addPolygons()!

# create color palette
npal <- colorNumeric("YlOrRd", nhoods$pop17_den)

# create leaflet object
nhoods %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    color = "#444444", 
    weight = 1, 
    opacity = 1.0, 
    smoothFactor = 0.5,
    fillOpacity = 0.5,
    fillColor = ~npal(pop17_den),
    highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE),
    popup = paste("<b>Name:</b> ", nhoods$NHD_NAME, "<br>",
                  "<b>2017 Population:</b> ", round(nhoods$pop17, digits = 0), "<br>",
                  "<b>2017 Population per Square Kilometer:</b> ", round(nhoods$pop17_den, digits = 2))) %>%
    addLegend(pal = npal, values = ~pop17_den, opacity = .5, title = "Population Density (2017)")

As a review, the color palette we’re using comes from the RColorBrewer package. We can use RColorBrewer::display.brewer.all() to identify other color ramps:

display.brewer.all(type = "seq")

One final complication we’ll discuss is adding the city’s outline on top of the addPolygons() function. We can’t use a second instance of addPolygons() in our call, and so we’ll use the addPolylines() function instead. This allows us to display the boundary as a line instead of a polygon feature.

# create color palette
npal <- colorNumeric("Blues", nhoods$pop17_den)

# create leaflet object
nhoods %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    color = "#444444", 
    weight = 1, 
    opacity = 1.0, 
    smoothFactor = 0.5,
    fillOpacity = 0.5,
    fillColor = ~npal(pop17_den),
    highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE),
    popup = paste("<b>Name:</b> ", nhoods$NHD_NAME, "<br>",
                  "<b>2017 Population:</b> ", round(nhoods$pop17, digits = 0), "<br>",
                  "<b>2017 Population per Square Kilometer:</b> ", round(nhoods$pop17_den, digits = 2))) %>%
    addPolylines(
      data = city,
      color = "#000000",
      weight = 3
    ) %>%
    addLegend(pal = npal, values = ~pop17_den, opacity = .5, title = "Population Density (2017)")

For our final leaflet map, write your own code to map the 1950 population density of neighborhoods:

# create color palette
npal <- colorNumeric("RdPu", nhoods$pop50_den)

# create map
nhoods %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    color = "#444444", 
    weight = 1, 
    smoothFactor = 0.5,
    opacity = 1.0, 
    fillOpacity = 0.5,
    fillColor = ~npal(pop50_den),
    highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE),
    popup = paste("<b>Name:</b> ", nhoods$NHD_NAME, "<br>",
                  "<b>1950 Population:</b> ", round(nhoods$pop50, digits = 0), "<br>",
                  "<b>1950 Population per Square Kilometer:</b> ", 
                      round(nhoods$pop50_den, digits = 2))) %>%
    addPolylines(
      data = city,
      color = "#000000",
      weight = 3
    ) %>%
    addLegend(pal = npal, values = ~pop50_den, opacity = .5, 
              title = "Population Density (1950)")