previous | blogs | next

tmap: quick and easy thematic mapping in R

In this tutorial we are going to have a look at a very quick and very easy way to create a choropleth map with the Thematic Maps package tmap in R. This package is created by Martijn Tennekes.

You can recreate the map above, including the inset, using only a few lines of code.

As input for this exercise you will use a shapefile with municipal data which can be downloaded from the website of the French Institut national de l’information géographique et forestière (IGN). On this page you will find the dataset: GEOFLA® 2016 v2.2 Communes France Métropolitaine.

Download and unzipg this file (e.g. using 7-Zip). Browse down the many subfolders in this archive until you find the shapefile COMMUNE.SHP (together with it’s corresponding *.DBF *.PRJ and *.SHX files). Set the folder containing these files as your R working directory.

You will need the three packages below. If necessary you should install them first, using install.packages().

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library(tmap)
library(grid)

Some first explorations

Let’s start by loading the data:

FR_Communes_2016 <- st_read("COMMUNE.shp")
## Reading layer `COMMUNE' from data source `C:\IGN\Donnees\GEOFLA_2-2_COMMUNE_SHP_LAMB93_FXX_2016-06-28\GEOFLA\1_DONNEES_LIVRAISON_2016-06-00236\GEOFLA_2-2_SHP_LAMB93_FR-ED161\COMMUNE\COMMUNE.shp' using driver `ESRI Shapefile'
## Simple feature collection with 35798 features and 17 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 99217.1 ymin: 6049646 xmax: 1242417 ymax: 7110480
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

You might want to have a quick look at the data by plotting them:

plot(st_geometry(FR_Communes_2016))

When you are using RStudio you can have a look at the attribute data using the Data Viewer:

View(FR_Communes_2016)

Yes, as you can see, historically there are quite a few municipalities in France - more than 35,000 of them, some with very small population numbers (and hence low densities).

Now we can calculate the total population of metropolitan France (i.e. not including overseas departments and territories):

sum(FR_Communes_2016$POPULATION)
## [1] 63697865

A population of almost 64 million inhabitants.

The dataset does contain 17 columns with attribute data (plus a geometry column), but no column with the population density. You do have POPULATION and SUPERFICIE (i.e. area), in hectares. To get the area in square kilometers you have to divide it by 100. Create a new column DENS_POP using this command:

FR_Communes_2016$DENS_POP <- round((FR_Communes_2016$POPULATION / (FR_Communes_2016$SUPERFICIE / 100)))
summary(FR_Communes_2016$DENS_POP)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    18.0    40.0   177.6    94.0 41929.0

Density ranges from 0 to 42,000 inhabitants per square kilometre.

Now you can create your thematic map. After playing with the different options we came up with the code below:

tm_shape(FR_Communes_2016) +
  tm_fill("DENS_POP", style = "kmeans", n = 10, title=expression("Inhabitants/" * km^2)) +
  tm_credits("data from http://www.ign.fr/", position = c("center", "bottom")) +
  tm_layout(title = "Population Density\nFrance, 2016", bg.color="lightgray", legend.title.size=.8)

Not bad for a first attempt, is it? Depending on your knowledge of the French topography you might be able to point cities like Bordeaux, Toulouse, Lyon and a few others on the map. And of course: Paris! Yes, the City of Light is clearly recognizable in this plot.

Further steps

Based on the insights gained from this attempt there was the wish to have a closer look at the region around Paris, the region Île-de-France. We created a subset ile_de_france:

ile_de_france <- FR_Communes_2016[FR_Communes_2016$NOM_REG == "ILE-DE-FRANCE",]

And we compared the population of this region with the whole of France. Almost 20 percent of the French population lives in Paris or it’s direct surroundings:

sum(ile_de_france$POPULATION)
## [1] 11959807
sum(ile_de_france$POPULATION)/sum(FR_Communes_2016$POPULATION) * 100
## [1] 18.77584

Hence the wish to create an inset on the map to highlight this population concentration.

To create space for this addition we had to modify the layout of the map. Apart from that we also changed the color scheme and the layout of the legend (for all the details: please refer to the code sample below).

And to mark the Paris region on the map we needed a rectangle indicating it’s location. We created this Boundig Box using the following line of code:

bb_idf <- st_sf(st_as_sfc(st_bbox(ile_de_france)))

And we manually rounded the range values suggested by kmeans (in the first plot):

breaks_pop <- c(0,200,500,1200,2200,3700,6000,10000,15000,23000,42000)

With some trial and error we came up with the following code for the entire layout with the map of France (with bounding box), leaving space for our inset. Please note that we do not plot the map immediately, but store it in a tmap object (map_france) first:

map_france <- tm_shape(FR_Communes_2016) +
  tm_fill("DENS_POP", style = "fixed", breaks = breaks_pop,
          title=expression("Inhabitants/" * km^2), palette = "Greens") +
  tm_shape(bb_idf) +
  tm_polygons(alpha = 0, border.col = "grey10") +
  tm_credits("data from http://www.ign.fr/", col = "grey10", position = c("center", "bottom")) +
  tm_layout("Population Density\nFrance, 2016", bg.color="grey75", legend.title.size=.8,
            legend.position = c("left", "top"), legend.format = c(text.align = "right", text.separator = "-"),
            outer.margins=c(.05,0,.05,0), inner.margins=c(.02,.25,.02,.02), asp=0, frame = FALSE)

Next we created a viewport for our inset:

viewport_ile_de_france <- viewport(x = .2, y = .25, width = .45, height = .45)

And we created the map of Île-de-France. Yes, we made sure to use the same range values (breaks) to be able to refer to the same legend:

map_ile_de_france <- tm_shape(ile_de_france) +
  tm_fill("DENS_POP", style="fixed", breaks = breaks_pop, palette = "Greens") +
  tm_layout("Île-de-France\n(Paris region)", bg.color="grey75", legend.show = FALSE, frame = "grey10",
            title.size = 0.7, title.position = c("left", "bottom"))

And then we were ready to print the entire map and add the inset. The result is the same as the map shown before, at the top of this page:

map_france
## Warning in legend.format$scientific <- FALSE: Coercing LHS to a list
print(map_ile_de_france, vp = viewport_ile_de_france)


Egge-Jan Pollé