4. Visualizing spatial data with the tmap package
Foreword on working directory, data and packages
The Data for this tutorial are provided via Github. In order to reproudce the code with your own data, replace the url with your local filepath
If you work locally, you may define the working directory
setwd("YOUR/FILEPATH")
# this will not work. if you want to use it, define a path on your computerThis should be the place where your data and other folders are stored. more Information can be found here
Code chunks for saving data and results are disabeled in this tutorial (using a #). Please replace (“YOUR/FILEPATH”) with your own path.
The code is written in such a way, that packages will be automatically installed (if not already installed) and loaded. Packages will be stored in the default location of your computer.
Intro
The tmap package (short for thematic map) offers a wide variety of approaches to create thematic maps. The folowwing tutorial introduces static and dynamic maps created with the in-packaga data. It is inspired by the Introduction from the CRAN website.
Load the tmap package
if (!require(tmap)){install.packages("tmap"); library(tmap)}We often use the tmap package for so-called “quick and dirty mapping” during data exploration.
Get tips and tricks. Random examples with in-package data
tmap_tip()## Use size.max to set the value mapped to maximum symbol/bubble size. Use scale to increase (or decrease) all symbols/bubbles.
## New since tmap 1.1
##
## data(metro)
## tm_shape(metro) + tm_bubbles("pop2020", size.max = 5e7, scale = 2) + tm_legend(outside = TRUE, outside.position = "bottom")## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
Basic static Maps
Load data from tmap package.
data("World")Check the structure and format of the data set
class(World)## [1] "sf" "data.frame"str(World)## Classes 'sf' and 'data.frame': 177 obs. of 16 variables:
## $ iso_a3 : Factor w/ 177 levels "AFG","AGO","ALB",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ name : Factor w/ 177 levels "Afghanistan",..: 1 4 2 166 6 7 5 56 8 9 ...
## $ sovereignt : Factor w/ 171 levels "Afghanistan",..: 1 4 2 159 6 7 5 52 8 9 ...
## $ continent : Factor w/ 8 levels "Africa","Antarctica",..: 3 1 4 3 8 3 2 7 6 4 ...
## $ area : Units: [km^2] num 652860 1246700 27400 71252 2736690 ...
## $ pop_est : num 28400000 12799293 3639453 4798491 40913584 ...
## $ pop_est_dens: num 43.5 10.3 132.8 67.3 15 ...
## $ economy : Factor w/ 7 levels "1. Developed region: G7",..: 7 7 6 6 5 6 6 6 2 2 ...
## $ income_grp : Factor w/ 5 levels "1. High income: OECD",..: 5 3 4 2 3 4 2 2 1 1 ...
## $ gdp_cap_est : num 784 8618 5993 38408 14027 ...
## $ life_exp : num 59.7 NA 77.3 NA 75.9 ...
## $ well_being : num 3.8 NA 5.5 NA 6.5 4.3 NA NA 7.2 7.4 ...
## $ footprint : num 0.79 NA 2.21 NA 3.14 2.23 NA NA 9.31 6.06 ...
## $ inequality : num 0.427 NA 0.165 NA 0.164 ...
## $ HPI : num 20.2 NA 36.8 NA 35.2 ...
## $ geometry :sfc_MULTIPOLYGON of length 177; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:69, 1:2] 61.2 62.2 63 63.2 64 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:15] "iso_a3" "name" "sovereignt" "continent" ...World## Simple feature collection with 177 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## First 10 features:
## iso_a3 name sovereignt continent
## 1 AFG Afghanistan Afghanistan Asia
## 2 AGO Angola Angola Africa
## 3 ALB Albania Albania Europe
## 4 ARE United Arab Emirates United Arab Emirates Asia
## 5 ARG Argentina Argentina South America
## 6 ARM Armenia Armenia Asia
## 7 ATA Antarctica Antarctica Antarctica
## 8 ATF Fr. S. Antarctic Lands France Seven seas (open ocean)
## 9 AUS Australia Australia Oceania
## 10 AUT Austria Austria Europe
## area pop_est pop_est_dens economy
## 1 652860.000 [km^2] 28400000 4.350090e+01 7. Least developed region
## 2 1246700.000 [km^2] 12799293 1.026654e+01 7. Least developed region
## 3 27400.000 [km^2] 3639453 1.328268e+02 6. Developing region
## 4 71252.172 [km^2] 4798491 6.734519e+01 6. Developing region
## 5 2736690.000 [km^2] 40913584 1.495003e+01 5. Emerging region: G20
## 6 28470.000 [km^2] 2967004 1.042151e+02 6. Developing region
## 7 12259213.973 [km^2] 3802 3.101341e-04 6. Developing region
## 8 7257.455 [km^2] 140 1.929051e-02 6. Developing region
## 9 7682300.000 [km^2] 21262641 2.767744e+00 2. Developed region: nonG7
## 10 82523.000 [km^2] 8210281 9.949082e+01 2. Developed region: nonG7
## income_grp gdp_cap_est life_exp well_being footprint inequality
## 1 5. Low income 784.1549 59.668 3.8 0.79 0.42655744
## 2 3. Upper middle income 8617.6635 NA NA NA NA
## 3 4. Lower middle income 5992.6588 77.347 5.5 2.21 0.16513372
## 4 2. High income: nonOECD 38407.9078 NA NA NA NA
## 5 3. Upper middle income 14027.1261 75.927 6.5 3.14 0.16423830
## 6 4. Lower middle income 6326.2469 74.446 4.3 2.23 0.21664810
## 7 2. High income: nonOECD 200000.0000 NA NA NA NA
## 8 2. High income: nonOECD 114285.7143 NA NA NA NA
## 9 1. High income: OECD 37634.0832 82.052 7.2 9.31 0.08067825
## 10 1. High income: OECD 40132.6093 81.004 7.4 6.06 0.07129351
## HPI geometry
## 1 20.22535 MULTIPOLYGON (((61.21082 35...
## 2 NA MULTIPOLYGON (((16.32653 -5...
## 3 36.76687 MULTIPOLYGON (((20.59025 41...
## 4 NA MULTIPOLYGON (((51.57952 24...
## 5 35.19024 MULTIPOLYGON (((-65.5 -55.2...
## 6 25.66642 MULTIPOLYGON (((43.58275 41...
## 7 NA MULTIPOLYGON (((-59.57209 -...
## 8 NA MULTIPOLYGON (((68.935 -48....
## 9 21.22897 MULTIPOLYGON (((145.398 -40...
## 10 30.47822 MULTIPOLYGON (((16.97967 48...names(World)## [1] "iso_a3" "name" "sovereignt" "continent" "area"
## [6] "pop_est" "pop_est_dens" "economy" "income_grp" "gdp_cap_est"
## [11] "life_exp" "well_being" "footprint" "inequality" "HPI"
## [16] "geometry"Summarize data
summary(World)## iso_a3 name sovereignt
## AFG : 1 Afghanistan: 1 France : 3
## AGO : 1 Albania : 1 Denmark : 2
## ALB : 1 Algeria : 1 Israel : 2
## ARE : 1 Angola : 1 United Kingdom : 2
## ARG : 1 Antarctica : 1 United States of America: 2
## ARM : 1 Argentina : 1 Afghanistan : 1
## (Other):171 (Other) :171 (Other) :165
## continent area pop_est pop_est_dens
## Africa :51 Min. : 2590 Min. :1.400e+02 Min. : 0.0003
## Asia :47 1st Qu.: 42390 1st Qu.:3.442e+06 1st Qu.: 20.8203
## Europe :39 Median : 183630 Median :9.036e+06 Median : 65.3438
## North America:18 Mean : 814568 Mean :3.827e+07 Mean : 105.5498
## South America:13 3rd Qu.: 622980 3rd Qu.:2.595e+07 3rd Qu.: 116.8924
## Oceania : 7 Max. :16376870 Max. :1.339e+09 Max. :1198.8237
## (Other) : 2
## economy income_grp
## 1. Developed region: G7 : 7 1. High income: OECD :32
## 2. Developed region: nonG7:32 2. High income: nonOECD:17
## 3. Emerging region: BRIC : 4 3. Upper middle income :44
## 4. Emerging region: MIKT : 4 4. Lower middle income :47
## 5. Emerging region: G20 :19 5. Low income :37
## 6. Developing region :66
## 7. Least developed region :45
## gdp_cap_est life_exp well_being footprint
## Min. : 300.5 Min. :48.91 Min. :2.867 Min. : 0.610
## 1st Qu.: 2215.4 1st Qu.:65.04 1st Qu.:4.575 1st Qu.: 1.425
## Median : 7394.8 Median :73.24 Median :5.200 Median : 2.604
## Mean : 14977.1 Mean :70.80 Mean :5.412 Mean : 3.223
## 3rd Qu.: 19193.8 3rd Qu.:76.93 3rd Qu.:6.300 3rd Qu.: 4.482
## Max. :200000.0 Max. :83.24 Max. :7.800 Max. :15.820
## NA's :1 NA's :41 NA's :41 NA's :41
## inequality HPI geometry
## Min. :0.04322 Min. :12.78 MULTIPOLYGON :177
## 1st Qu.:0.13931 1st Qu.:21.21 epsg:4326 : 0
## Median :0.21293 Median :26.29 +proj=long...: 0
## Mean :0.23427 Mean :26.48
## 3rd Qu.:0.32932 3rd Qu.:31.73
## Max. :0.50734 Max. :44.71
## NA's :41 NA's :41Plot population estimate and check histogram and summarize to get an overview of the data.
plot(World$pop_est)
hist(World$pop_est)
summary(World$pop_est)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.400e+02 3.442e+06 9.036e+06 3.827e+07 2.595e+07 1.339e+09Plot a simple Map of HPI (Happy Planet Index)
tm_shape(World) +
tm_polygons("HPI")
Plot a simple Map of life expectancy
tm_shape(World) +
tm_polygons("life_exp")
Interactive Maps
The mode of plotting can be set with the function tmap_mode –> static image “plot” mode
tmap_mode("plot")## tmap mode set to plottingtm_shape(World) +
tm_polygons("life_exp")
Or viewed interactively using “view” modes –> uses leaflet TO DO
tmap_mode("view")
tm_shape(World) +
tm_polygons("life_exp")# l <- tmap_mode("view") +
# tm_shape(World) +
# tm_polygons("life_exp")
# frameWidget(l)Plot Multiple shapes and layers
We now combine several layers which are included in the tmap package:
- ‘metro’ population times series from 1950 to 2030. All metro areas >1 mio inhabitants
- ‘Land’ Spatial data of global land cover
- ‘rivers’ major rivers worldwide
We now want to plot them all together
data(World, metro, rivers, land)
tmap_mode("plot")## tmap mode set to plotting## tmap mode set to plotting
tm_shape(land) +
tm_raster("elevation", palette = terrain.colors(10)) +
tm_shape(World) +
tm_borders("white", lwd = .5)+
tm_text("iso_a3", size = "AREA")+
tm_shape(metro) +
tm_symbols(col = "red", size = "pop2020", scale = .5) +
tm_legend(show = FALSE)## old-style crs object detected; please recreate object with a recent sf::st_crs()## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
Facets - multiple plots in one window
There are multiple ways to do so:
By assigning multiple variable names to one aesthetic:
In this example we display “HPI” and “Economy” of the data set together and interactively (interactivity is turned off on the website due to package restrictions. It should work on your computer though).
tmap_mode("view")## tmap mode set to interactive viewingtm_shape(World) +
tm_polygons(c("HPI", "economy")) +
tm_facets(sync = TRUE, ncol = 2)By splitting the spatial data with the by argument of tm_facets
In this example we calculate the percentage of the male population in the prvinces of the Netherlands. The data are provided by the tmap package.
tmap_mode("plot")## tmap mode set to plotting## tmap mode set to plotting
data(NLD_muni)
NLD_muni$perc_men <- NLD_muni$pop_men / NLD_muni$population * 100
tm_shape(NLD_muni) +
tm_polygons("perc_men", palette = "RdYlBu") +
tm_facets(by = "province")## old-style crs object detected; please recreate object with a recent sf::st_crs()## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
By using the tmap_arrange function
In this example we calculate the population density of the communities of the Netherlands, simply by using the convert2density argument in the tm_polygons function. In it the area size is either approximated from the shape object, or given by the argument area.
Also, we use the tm_bubbles function to create proportionate bubbles.
tmap_mode("plot")## tmap mode set to plotting## tmap mode set to plotting
data(NLD_muni)
tm1 <- tm_shape(NLD_muni) + tm_polygons("population", convert2density = TRUE)
tm2 <- tm_shape(NLD_muni) + tm_bubbles(size = "population")
tmap_arrange(tm1, tm2)## old-style crs object detected; please recreate object with a recent sf::st_crs()## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()## Legend labels were too wide. Therefore, legend.text.size has been set to 0.32. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
Basemaps and overlay tile maps
Tiled basemaps can be added with the layer function tm_basemap. Semi-transparent overlay maps (for example annotation labels) can be added with tm_tiles.
tmap_mode("plot")## tmap mode set to plottingtm_basemap("Stamen.Watercolor") +
tm_shape(metro) + tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")## old-style crs object detected; please recreate object with a recent sf::st_crs()## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
To see all available basemaps for tm_tiles check link http://leaflet-extras.github.io/leaflet-providers/preview/ for example basemaps from ESRI
tm_basemap("Esri.WorldImagery") +
tm_shape(metro) + tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
#O ptions and style with tm_layout
tmap_mode("plot")## tmap mode set to plotting## tmap mode set to plotting
tm_shape(World) +
tm_polygons("HPI") +
tm_layout(bg.color = "skyblue", inner.margins = c(0, .02, .02, .02))
tmap_options(bg.color = "black", legend.text.color = "white")tm_shape(World) +
tm_polygons("HPI", legend.title = "Happy Planet Index")
tmap_style("natural")## tmap style set to "natural"## other available styles are: "white", "gray", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"tm_shape(World) +
tm_polygons("HPI", legend.title = "Happy Planet Index")
See what options have been changed
tmap_options_diff()## current tmap options (style "natural") that are different from default tmap options (style "white"):## $bg.color
## [1] "lightskyblue1"
##
## $aes.color
## fill borders symbols dots
## "darkolivegreen3" "black" "tomato2" "firebrick"
## lines text na null
## "steelblue" "black" "white" "grey70"
##
## $aes.palette
## $aes.palette$seq
## [1] "YlGn"
##
## $aes.palette$div
## [1] "RdYlGn"
##
## $aes.palette$cat
## [1] "Set3"
##
##
## $earth.boundary
## [1] TRUE
##
## $space.color
## [1] "white"
##
## $legend.frame
## [1] TRUE
##
## $legend.bg.color
## [1] "grey90"
##
## $basemaps
## [1] "Esri.NatGeoWorldMap"
##
## $basemaps.alpha
## [1] 1Reset the options to the default values
tmap_options_reset()## tmap options successfully resetExporting maps
tm <- tm_shape(World) +
tm_polygons("HPI", legend.title = "Happy Planet Index")save an image (“plot” mode) to the working directory.
tmap_save(tm, filename = "world_map.png")save as stand-alone HTML file (“view” mode)
tmap_save(tm, filename = "world_map.html")Quick and dirty mapping
Use the qtm function if you want a simply map fast
qtm(World, fill = "HPI", fill.pallete = "RdYlGn")