Synopsis

In this document, I showcase the tmap and leaflet packages in R. The tmap package is primarily used to create static maps, that can be saved as images (.png). But an output object of tmap can also be saved as an interactive webmap (.html). In this case, leaflet is used for generating the interactive map. Leaflet is a powerful JavaScript library for creating interactive webmaps that can be integrated into R with the leaflet package.

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.4.3
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.4.3
## Warning: package 'sp' was built under R version 3.4.3
library(RColorBrewer)
library(tmap)
## Warning: package 'tmap' was built under R version 3.4.3
library(leaflet)
## Warning: package 'leaflet' was built under R version 3.4.3

Loading data from the tmap package and show it on static maps

The tmap package contains some sample data as sp objects, that is spatial objects with an associated data frame. With str() and the argument max.level = 2 you can get an overwiev of the data contained in an sp object. sp objects are of class S4 and contain slots that can be adressed with the @ sign.

# explore Europe data (polygon)
data(Europe)
summary(Europe)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##       min     max
## x 2500000 7350000
## y 1400000 5750000
## Is projected: TRUE 
## proj4string :
## [+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
## +ellps=GRS80 +units=m +no_defs]
## Data attributes:
##      iso_a3        name             sovereignt         continent 
##  ALA    : 1   Aland  : 1   United Kingdom: 4   Africa       : 3  
##  ALB    : 1   Albania: 1   Denmark       : 3   Asia         :14  
##  AND    : 1   Algeria: 1   Finland       : 2   Europe       :50  
##  ARM    : 1   Andorra: 1   Israel        : 2   North America: 1  
##  AUT    : 1   Armenia: 1   Albania       : 1                     
##  (Other):62   Austria: 1   Algeria       : 1                     
##  NA's   : 1   (Other):62   (Other)       :55                     
##               part               EU_Schengen      area         
##  Northern Europe:15   EU Schengen      :22   Min.   :       1  
##  Western Europe : 9   EU Schengen cand.: 4   1st Qu.:   11567  
##  Southern Europe:17   EU Non-Schengen  : 2   Median :   55960  
##  Eastern Europe :10   Schengen Non-EU  : 4   Mean   :  449248  
##  NA's           :17   NA's             :36   3rd Qu.:  216465  
##                                              Max.   :16376870  
##                                              NA's   :17        
##     pop_est           pop_est_dens         gdp_md_est     
##  Min.   :      832   Min.   :    3.059   Min.   :    355  
##  1st Qu.:   581978   1st Qu.:   63.690   1st Qu.:  10316  
##  Median :  4676305   Median :  109.418   Median :  93750  
##  Mean   : 15800734   Mean   :  533.517   Mean   : 382015  
##  3rd Qu.: 10561130   3rd Qu.:  191.157   3rd Qu.: 334650  
##  Max.   :140041247   Max.   :16482.500   Max.   :2918000  
##  NA's   :17          NA's   :17          NA's   :17       
##   gdp_cap_est                           economy  
##  Min.   :  2470   2. Developed region: nonG7:36  
##  1st Qu.: 16768   6. Developing region      : 9  
##  Median : 29616   1. Developed region: G7   : 4  
##  Mean   : 37890   3. Emerging region: BRIC  : 1  
##  3rd Qu.: 40042   4. Emerging region: MIKT  : 1  
##  Max.   :426683   (Other)                   : 0  
##  NA's   :17       NA's                      :17  
##                    income_grp    life_exp       well_being   
##  1. High income: OECD   :25   Min.   :68.50   Min.   :4.180  
##  2. High income: nonOECD:11   1st Qu.:74.45   1st Qu.:5.102  
##  3. Upper middle income :11   Median :78.80   Median :5.840  
##  4. Lower middle income : 4   Mean   :77.23   Mean   :6.011  
##  5. Low income          : 0   3rd Qu.:80.50   3rd Qu.:6.959  
##  NA's                   :17   Max.   :82.30   Max.   :7.771  
##                               NA's   :29      NA's   :29     
##       HPI       
##  Min.   :28.27  
##  1st Qu.:37.41  
##  Median :41.28  
##  Mean   :41.34  
##  3rd Qu.:46.26  
##  Max.   :54.05  
##  NA's   :29
str(Europe, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  68 obs. of  16 variables:
##   ..@ polygons   :List of 68
##   ..@ plotOrder  : int [1:68] 57 26 66 38 67 20 17 63 52 32 ...
##   ..@ bbox       : num [1:2, 1:2] 2500000 1400000 7350000 5750000
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
head(Europe@data)
##    iso_a3       name sovereignt continent            part EU_Schengen
## 5     ALB    Albania    Albania    Europe Southern Europe        <NA>
## 6     ALA      Aland    Finland    Europe Northern Europe        <NA>
## 7     AND    Andorra    Andorra    Europe Southern Europe        <NA>
## 10    ARM    Armenia    Armenia      Asia            <NA>        <NA>
## 17    AUT    Austria    Austria    Europe  Western Europe EU Schengen
## 18    AZE Azerbaijan Azerbaijan      Asia            <NA>        <NA>
##          area pop_est pop_est_dens gdp_md_est gdp_cap_est
## 5  27400.0000 3639453    132.82675      21810    5992.659
## 6    674.3813   27153     40.26357       1563   57562.700
## 7    470.0000   83888    178.48511       3660   43629.601
## 10         NA      NA           NA         NA          NA
## 17 82409.0000 8210281     99.62845     329500   40132.609
## 18         NA      NA           NA         NA          NA
##                       economy              income_grp life_exp well_being
## 5        6. Developing region  4. Lower middle income     76.9   5.268937
## 6  2. Developed region: nonG7    1. High income: OECD       NA         NA
## 7  2. Developed region: nonG7 2. High income: nonOECD       NA         NA
## 10                       <NA>                    <NA>       NA         NA
## 17 2. Developed region: nonG7    1. High income: OECD     80.9   7.346036
## 18                       <NA>                    <NA>       NA         NA
##         HPI
## 5  54.05118
## 6        NA
## 7        NA
## 10       NA
## 17 47.08514
## 18       NA
# display population density with default (pretty) classification
tm_shape(Europe) +
  tm_fill("pop_est_dens", id="name", title="Population Density") +
  tm_borders("grey20") + tm_text("iso_a3", size = 0.5) + tm_layout(legend.bg.color = "white")

# display population density with quantile classification
tm_shape(Europe) +
  tm_fill("pop_est_dens", id="name", title="Population Density", style = "quantile") +
  tm_borders("grey20") + tm_text("iso_a3", size = 0.5) + tm_layout(legend.bg.color = "white")

# display categocial economy
tm_shape(Europe) +
  tm_fill("economy", id="name", title="Economy Category") +
  tm_borders("grey20") + tm_text("iso_a3", size = 0.5) + tm_layout(legend.bg.color = "white")

# display line data
data(rivers)
summary(rivers)
## Object of class SpatialLinesDataFrame
## Coordinates:
##          min      max
## x -165.24394 176.3258
## y  -50.24014  73.3349
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs]
## Data attributes:
##        name                   type        scalerank     strokelwd     
##  Volga   :  24   Lake Centerline: 218   Min.   :0.0   Min.   : 1.500  
##  Paran   :  22   River          :1399   1st Qu.:3.0   1st Qu.: 2.000  
##  Colorado:  21   NA's           :   1   Median :5.0   Median : 2.500  
##  Niger   :  20                          Mean   :4.4   Mean   : 3.046  
##  Columbia:  19                          3rd Qu.:6.0   3rd Qu.: 3.500  
##  (Other) :1472                          Max.   :6.0   Max.   :12.000  
##  NA's    :  40
str(rivers, max.level = 2)
## Formal class 'SpatialLinesDataFrame' [package "sp"] with 4 slots
##   ..@ data       :'data.frame':  1618 obs. of  4 variables:
##   ..@ lines      :List of 1618
##   ..@ bbox       : num [1:2, 1:2] -165.2 -50.2 176.3 73.3
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
head(rivers@data)
##          name            type scalerank strokelwd
## 0        <NA>           River         5       2.0
## 1        Ebro           River         5       1.5
## 2        Ebro           River         5       2.0
## 3        Ebro           River         5       2.5
## 4        Ebro           River         5       3.0
## 5 Mississippi Lake Centerline         5       2.0
tm_shape(Europe) +
  tm_fill() + tm_borders("grey20") +
  tm_shape(rivers) + tm_lines(col = "dodgerblue3") +
  tm_add_legend(type = "line", col = "dodgerblue3", labels = "Rivers", title = "Europe River Map") +
  tm_layout(legend.bg.color = "white")

# display point data as bubbles
data(metro)
summary(metro)
## Object of class SpatialPointsDataFrame
## Coordinates:
##                 min       max
## Longitude -123.1193 174.76667
## Latitude   -37.8140  60.16925
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0]
## Number of points: 436
## Data attributes:
##      name            name_long            iso_a3         
##  Length:436         Length:436         Length:436        
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##     pop1950            pop1960            pop1970        
##  Min.   :    3127   Min.   :    7192   Min.   :   16540  
##  1st Qu.:  123885   1st Qu.:  208500   1st Qu.:  339692  
##  Median :  282090   Median :  432850   Median :  632775  
##  Mean   :  659393   Mean   :  917510   Mean   : 1236404  
##  3rd Qu.:  696197   3rd Qu.:  948813   3rd Qu.: 1284255  
##  Max.   :12338471   Max.   :16678821   Max.   :23297503  
##     pop1980            pop1990            pop2000        
##  Min.   :   38657   Min.   :  138179   Min.   :  388603  
##  1st Qu.:  528870   1st Qu.:  763620   1st Qu.: 1066005  
##  Median :  901310   Median : 1129506   Median : 1437566  
##  Mean   : 1593220   Mean   : 2017844   Mean   : 2584675  
##  3rd Qu.: 1652315   3rd Qu.: 2095953   3rd Qu.: 2648557  
##  Max.   :28548512   Max.   :32530003   Max.   :34449908  
##     pop2010            pop2020            pop2030        
##  Min.   : 1003320   Min.   :  910051   Min.   :  821440  
##  1st Qu.: 1329874   1st Qu.: 1598478   1st Qu.: 1843849  
##  Median : 1806583   Median : 2361189   Median : 2744112  
##  Mean   : 3201460   Mean   : 3921452   Mean   : 4574673  
##  3rd Qu.: 3289534   3rd Qu.: 3965543   3rd Qu.: 4888469  
##  Max.   :36833979   Max.   :38323229   Max.   :37190489
str(metro, max.level = 2)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  436 obs. of  12 variables:
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:436, 1:2] 69.17 3.04 13.23 -58.4 -64.18 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ bbox       : num [1:2, 1:2] -123.1 -37.8 174.8 60.2
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
head(metro@data)
##            name             name_long iso_a3 pop1950 pop1960 pop1970
## 2         Kabul                 Kabul    AFG  170784  285352  471891
## 8       Algiers El Djazair  (Algiers)    DZA  516450  871636 1281127
## 13       Luanda                Luanda    AGO  138413  219427  459225
## 16 Buenos Aires          Buenos Aires    ARG 5097612 6597634 8104621
## 17      Cordoba               Cordoba    ARG  429249  605309  809794
## 25      Rosario               Rosario    ARG  554483  671349  816230
##    pop1980  pop1990  pop2000  pop2010  pop2020  pop2030
## 2   977824  1549320  2401109  3722320  5721697  8279607
## 8  1621442  1797068  2140577  2432023  2835218  3404575
## 13  771349  1390240  2591388  4508434  6836849 10428756
## 16 9422362 10513284 12406780 14245871 15894307 16956491
## 17 1009521  1200168  1347561  1459268  1562509  1718192
## 25  953491  1083819  1152387  1298073  1453814  1606993
tm_shape(Europe) +
  tm_polygons("pop_est_dens", id="name", title="Population Density", style = "quantile") +
  tm_shape(metro) +  tm_bubbles("pop2030", col = "blue", title.size = "City Population 2030") +
  tm_text("name", size = 0.5, legend.size.show = FALSE, root = 8, size.lowerbound = .7, auto.placement = TRUE) +
  tm_layout(legend.bg.color = "white")

Using leaflet to create interactive webmaps

The leaflet function addTiles() adds the default OpenStreetMap base layer to a map. With addProviderTiles() you can add other base layers. With the argument “group” you can assign a layer to a group. These layer groups can be referenced by the addLayersControl() function to generate radio buttons or check boxes to select layers to show or hide in the map.

The interactive leaflet maps expect geo-data with a WGS84 projection. The data set Europe of the tmap package has another projection and needs to be projected to WGS84. For choropleth maps in leaflet you need to specify a color palette (color ramp) for each displayed layer. You can use predefined color palettes of ColorBrewer or Viridis. I will show the countries of Europe colored based on their population density.

# transform projection of the europe sp to wgs84
proj4string(metro)
## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
proj4string(Europe)
## [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
europe_wgs84 <- spTransform(Europe, "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

# show all ColorBrewer ramps
display.brewer.all()

# create a color palette for population density
pal_pop_dens <- colorQuantile("viridis", europe_wgs84$pop_est_dens, n = 7)

Again, I want to show the estimated city population in 2030 of the largest cities in the world as circles with varying size. The argument “radius” of addCircles() expects a measurement in meters. I use the sqrt() of pop2030 to make the areas of the circles proportional to pop2030 and multiply by 50 to get an appropriate range of sizes. To make sure, large circles don’t obscure smaller ones, I first sort the data in descending order by pop2030.

To showcase the ability of leaflet to display clustered markers I also plot some earth quake point data. Just for fun, I also incude a convex hull of all the earth quake points as an additional layer. The earth quake data is not a spatial object but a simple data frame. It can still be displayed on a map, because it contains the latitude and longitude in two columns.

# sort the metro data descending by pop2030 to display small circles on top of large ones
metro_sorted <- metro[order(-metro$pop2030), ]

# use earth quake point data as example for clustered markers
str(quakes)
## 'data.frame':    1000 obs. of  5 variables:
##  $ lat     : num  -20.4 -20.6 -26 -18 -20.4 ...
##  $ long    : num  182 181 184 182 182 ...
##  $ depth   : int  562 650 42 626 649 195 82 194 211 622 ...
##  $ mag     : num  4.8 4.2 5.4 4.1 4 4 4.8 4.4 4.7 4.3 ...
##  $ stations: int  41 15 43 19 11 12 43 15 35 19 ...
# generate an outline surrounding all earth quake point data
quakes_outline <- quakes[chull(quakes$long, quakes$lat), ]

Putting it all together: leaflet works with the pipe notation and add layer one at the time. First, I add several base layer, then the discussed data layers. For each data layer, you can add a separate legend. With the addLayersControl() function you specify which layer groups can be selected or deselected.

leaflet(width = "100%", height = "800px") %>%
  # base layers
  addTiles(group = "OpenStreeMap.Default") %>%
  addProviderTiles("OpenStreetMap.BlackAndWhite", group = "OpenStreetMap.BlackAndWhite") %>%
  addProviderTiles("OpenStreetMap.DE", group = "OpenStreetMap.DE") %>%
  addProviderTiles("OpenStreetMap.HOT", group = "OpenStreetMap.HOT") %>%
  addProviderTiles("OpenTopoMap", group = "OpenTopoMap") %>%
  addProviderTiles("Thunderforest.Transport", group = "Thunderforest.Transport") %>%
  addProviderTiles("OpenMapSurfer.Roads", group = "OpenMapSurfer.Roads") %>%
  addProviderTiles("OpenMapSurfer.Grayscale", group = "OpenMapSurfer.Grayscale") %>%
  addProviderTiles("Stamen.Toner", group = "Stamen.Toner") %>%
  addProviderTiles("Stamen.TonerBackground", group = "Stamen.TonerBackground") %>%
  addProviderTiles("Stamen.TonerLite", group = "Stamen.TonerLite") %>%
  addProviderTiles("Stamen.TerrainBackground", group = "Stamen.TerrainBackground") %>% 
  addProviderTiles("Esri.WorldStreetMap", group = "Esri.WorldStreetMap") %>%
  addProviderTiles("Esri.WorldTopoMap", group = "Esri.WorldTopoMap") %>%
  addProviderTiles("Esri.WorldGrayCanvas", group = "Esri.WorldGrayCanvas") %>%
  addProviderTiles("Esri.WorldImagery", group = "Esri.WorldImagery") %>%
  addProviderTiles("CartoDB.Positron", group = "CartoDB.Positron") %>%
  addProviderTiles("CartoDB.PositronNoLabels", group = "CartoDB.PositronNoLabels") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "CartoDB.DarkMatter") %>%
  addProviderTiles("CartoDB.DarkMatterNoLabels", group = "CartoDB.DarkMatterNoLabels") %>%
  # data layers
  # countries of Europe with population density
  addPolygons(data = europe_wgs84, weight = 1, color = "black", smoothFactor = 0.3,
              fillOpacity = 0.8, fillColor = ~pal_pop_dens(pop_est_dens),
              label = ~paste0(name, ": ", formatC(pop_est_dens, big.mark = "'")),
              group = "Population Density (Europe)") %>%
  addLegend(pal = pal_pop_dens, values = europe_wgs84$pop_est_dens, opacity = 0.8,
            labFormat = labelFormat(big.mark = "'"),
            title = "Population Density") %>%
  # largest cities in the World with estimated population in 2030
  addCircles(data = metro_sorted, radius = ~sqrt(pop2030)*50, stroke = FALSE, col = "red",
             label = ~paste0(name, ": ", formatC(pop2030/1000000, big.mark = "'"), " Mio."),
             group = "Cities (Population 2030)") %>%
  addLegend(labels = "Pop. Est. 2030", color = "red", title = "City Population 2030") %>%
  # earth quake data with magnitude
  addMarkers(data = quakes, lng = ~long, lat = ~lat,
             popup = ~paste0("Magnitude: ", mag),
             clusterOptions = markerClusterOptions(),
             group = "Earthquakes") %>%
  addPolygons(data = quakes_outline, lng = ~long, lat = ~lat, fill = FALSE,
              weight = 3, color = "#FFFFCC", group = "Earthquakes Outline") %>%
  # layers control
  addLayersControl(
    baseGroups = c("OpenStreeMap.Default", "OpenStreetMap.BlackAndWhite", "OpenStreetMap.DE",
                   "OpenStreetMap.HOT", "OpenTopoMap", "Thunderforest.Transport", "OpenMapSurfer.Roads",
                   "OpenMapSurfer.Grayscale", "Stamen.Toner", "Stamen.TonerBackground",
                   "Stamen.TonerLite", "Stamen.TerrainBackground", "Esri.WorldStreetMap",
                   "Esri.WorldTopoMap", "Esri.WorldImagery", "CartoDB.Positron",
                   "CartoDB.PositronNoLabels", "CartoDB.DarkMatter", "CartoDB.DarkMatterNoLabels"),
    overlayGroups = c("Population Density (Europe)", "Cities (Population 2030)",
                      "Earthquakes", "Earthquakes Outline"),
    position = "topleft",
    options = layersControlOptions(collapsed = FALSE)
  )