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
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")
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)
)