1 A brief introduction to the tmap package

The tmap package is a brand new easy way to plot thematic maps in R. Thematic maps are geographical maps in which spatial data distributions are visualized. This package offers a flexible, layer-based, and easy to use approach to create thematic maps, such as choropleths and bubble maps. The syntax for creating plots is similar to that of ggplot2.

1.1 Two useful websites

1.2 A first example with tmap

library(pacman)
pacman::p_load(tmap, statR, tidyverse, sp, cartogram)
data(metro, World)

Let’s start with a simple first example included in the package. All we need ist the package tmap and the spatial data “World”, which is a SpatialPolygonsDataFrame and “metro”, which is a SpatialPointsDataFrame.

class(metro)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
class(World)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

The “World”-DataFrame includes informations on the population, the economy and so on of all countries in the world.

World %>% as.data.frame() %>% head()

It is a SpatialPolygonsDataFrame, so we can simply draw a map that contains this information. For instance “The countries of the world by happiness index”

tm_shape(World)+tm_polygons("HPI", palette="-Blues", contrast=.7, id="name", title="Happiness Index")

If we want to add more information to this map, we can do this simply by drawing more layers over this base map.

Let’s have a look into the metro-DataFrame. It shows the population of the metropolitan regions in the world, for example all metropolitan regions of Switzerland, i.e. Zurich.

metro %>% as.data.frame() %>% dplyr::filter(iso_a3=="CHE") %>% head()
tm_shape(subset(metro,metro$iso_a3=="CHE"))+tm_bubbles("pop2010")

OK, that’s not very meaningful. So let’s look at the whole world.

tm_shape(metro)+tm_bubbles("pop2010")

A much better picture can be obtained if the information on the population is presented in a somewhat better graphic form and supplemented by a growth rate of the population figures. So, let’s add some color and a nicer looking legend.

metro$growth <- (metro$pop2020 - metro$pop2010) / (metro$pop2010 * 10) * 100

    tm_shape(metro) +
    tm_bubbles("pop2010", col = "growth", 
 breaks=c(-Inf, seq(0, 6, by=2), Inf),
               palette="-RdYlBu", contrast=1, 
               title.size="Metro population", 
               title.col="Growth rate (%)", id="name") + 
    tm_format_World_wide()

Now we can simply put the two maps (the happiness index world map an the population growth bubbles) on top of each other and get a simple yet informative map.

data(World, metro)
metro$growth <- (metro$pop2020 - metro$pop2010) / (metro$pop2010 * 10) * 100

mapWorld <- tm_shape(World) +
    tm_polygons("income_grp", palette="-Blues", contrast=.7, id="name", title="Income group") +
    tm_shape(metro) +
    tm_bubbles("pop2010", col = "growth", 
               border.col = "black", border.alpha = .5, 
               style="fixed", breaks=c(-Inf, seq(0, 6, by=2), Inf),
               palette="-RdYlBu", contrast=1, 
               title.size="Metro population", 
               title.col="Growth rate (%)", id="name") + 
    tm_style_gray() + tm_format_World_wide()

mapWorld

So this is our first static map. You can’t see details, but you can already see a lot of things from an overview point of view. But if you want to look at details of individual regions or metropolises, you need something more, namely an interactive version, which allows you to zoom into your point of interest.

This is very easy with the tmap package. The only thing to do is to change from the plot mode to the view mode. Afterwards we easyly plot the existing map a second time. Nad abracadabra you have an interactive leaflet map, where you can zoom in and some additional informations in a tooltip.

# set mode to view:
tmap_mode("view")
## tmap mode set to interactive viewing
mapWorld

1.3 An example for Zurich (Numbers of sold houses between 2007 und 2017)

This example includes all sales of detached houses since 2007. Because the data cannot be published for data protection reasons, I focus on the number of sold houses and artificial prices. Here is a first view on the desired result:

The aim of the map is to provide a simple overview of prices in the canton. Up to now, the amount of sold houses has been presented at regional level and, if necessary, at municipal level. However, a more differentiated overview is obtained when the data is aggregated into grids or hexbins.

Let’s start the new simple way

pacman::p_load(rgeos, rgdal, ggmap, leaflet)
tmap_mode("plot")
## tmap mode set to plotting

A first glance at the data set

class(hexbin_orig)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
class(hex_geb_fake)
## [1] "data.frame"
head(hex_geb_fake)

1.3.1 A simple map

And now let’s produce our map.

mapZH <- tm_shape(hexbin_zh) +
    tm_polygons(value, breaks=brekks, palette=rev(zhpal$zhlake),  colorNA =NULL, contrast=.8, title=anzeige,  border.col = "lightgrey")

mapZH

OK, now we have a map with an entry for every area where a detached house has been sold. However, it is not yet clear in which region of the world we are operating. For this you need a background map. Here we have different maps as shape-files available. So let’s take a map with the current municipal boundaries.

load(file = "ZHshape.RData")

#ZHshape <- rgdal::readOGR('L:/STAT/08_DS/03_GIS/Geodaten/Administrative_Grenzen/Gemeinden/Politische_Gemeinden_ZH/Shapes_Stände/Shape_generalisiert_2018/GEN_A4_GEMEINDEN_2018_F.shp')

And this is what this cantonal map now looks like:

qtm(ZHshape, fill="#b7c6e5")

We see a great little feature of tmap, the Quick Thematic Map (qtm), which simply plots a map. We can now add more layers to this simple command as we did with ggplot2.

qtm(ZHshape) + tm_compass() + tm_style_classic()+ tm_scale_bar()

Now let’s add the information we created earlier, so that we know in which areas in the canton of Zurich how many houses were sold. To make everything look good, I’ll change the legend a little bit and introduce you to the classic version of the tmap.

qtm(ZHshape, fill= "#b7c6e5")  +tm_compass() + tm_scale_bar()+
  tm_layout(scale = .5,
            legend.position = c(.78,.72), 
            legend.outside.size = 0.1,
            legend.title.size = 1.6,
            legend.height = 0.9,
            legend.text.size = 1.2, 
            legend.hist.size = 0.6)+
  mapZH

1.3.2 Interactive map with leaflet and tmap

So far so good. Another strength of tmap is the direct integration of the leaflet. We saw that earlier. Let’s try that for the cantonal map, too. To be able to display the informative tooltips here as well, I have to extend the function with which we produced the map earlier.

mapZH <- tm_shape(hexbin_zh) +
    tm_polygons(value, breaks=brekks, palette=rev(zhpal$zhlake), colorNA =NULL, contrast=.8, id="GEMEINDENA", title=anzeige, border.col = "lightgrey", border.alpha = NA, popup.vars=c("Numbers of Sales"="anz", "Price in CHF (Median)"="med_preis",  "ID"="GRID_ID"))
tmap_EFH <- tmap_leaflet(mapZH, mode = "view", show = FALSE) %>% setView(1249033,2682288,zoom=11)%>% fitBounds(-72, 40, -70, 43)%>% clearBounds() 
p_load(htmlwidgets)
#saveWidget(tmap_EFH, file="tmap_EFH_anz.html")
tmap_EFH

1.3.3 The cartogram function

And finally another great feature of the tmap package: The cartogram function. Many of you probably know these maps with color-filled polygons. At least I have used them more often in my last publications.

qtm(ZHshape2, fill="med_preis")

As you can see, this is also very easy to produce with the qtm function.

With the additional “cartogram” package you can add another dimension of information to this simple informative map. Here you see a map, whose coloring represents the price and where the size of the polygons additionally represents the number of sold objects. Especially in Winterthri many EInfamilienhäuser were sold, relative to the size of the municipality Winterthur.

library(cartogram)
ZHcarto<- cartogram(ZHshape2, "anz", itermax=5)
## Warning in cartogram(ZHshape2, "anz", itermax = 5): NA not allowed in
## weight vector. Features will be removed from Shape.
tm_shape(ZHcarto) + tm_fill("med_preis", style="jenks") + 
  tm_borders() + tm_layout(frame=F)