Creating Chloropleth Maps of Ile-de-France Departments using R: A Step-by-Step Guide

Author

Affiliation

Qatar University

Objectif

The objective of this tutorial is to demonstrate how to create a chloropleth map of the departments in the Ile-de-France region using R. Specifically, we will focus on representing the NEET indicator in eight Ile-de-France departments. It is important to note that NEET stands for Not in Education, Employment, or Training, which refers to individuals who are unemployed and not enrolled in any form of education or vocational training program.

What do you need for this tutorial?

To begin, you will need to download the shapefile for the Ile-de-France departments. You can obtain it by following this link:

https://www.data.gouv.fr/fr/datasets/r/2dedae24-eb12-4a92-803d-7ceddf08baf8

Once downloaded, create a new RStudio project and unzip the shapefiles. Then, place them in the folder corresponding to your RStudio project.

In addition, ensure that you have the following R packages installed:

require(rgdal)
require(ggplot2)
library(maptools)
library(sp)
library(shapefiles)
library(dplyr)

You may also require the following package and verify that the output of the provided R code displays TRUE. This package is necessary for constructing the polygons of your map.

library(gpclib)
gpclibPermit()
## [1] TRUE

Let’s Start

Step 1

Import the shapefile into R

dep_paris <- readShapePoly("geoflar-departements.shp")

Step 2

To proceed, add the variable of interest to the imported map data. In this case, the variable represents the NEET percentages in the eight Ile-de-France departments.

dep_paris@data$NEET=c(19,16,12.7,14.4,18.1,26.1,20,18.8)

Step 3

To begin, create the label variables for the maps. Since you will be creating two maps, one representing all departments and the other focusing on the departments in Paris, it is advisable to use different label variables. This is because the departments in Paris are smaller in size, and displaying the percentages directly on them may hinder readability of the map.

dep_paris@data$NEET_lb=paste0(dep_paris@data$NEET,"%")

dep_paris@data$NEET_lb_2=dep_paris@data$NEET_lb

dep_paris@data$NEET_lb[c(3,4,5,6)]=NA

Step 4

We will then use fortify to create the polygons

dep_paris_fr<- fortify(dep_paris,region = "id_geofla")

Step 5

Next, we will merge the variable that needs to be represented with the newly created data using the fortify function.

i=match(dep_paris_fr$id,dep_paris$id_geofla)
dep_paris_fr$NEET=dep_paris@data$NEET[i]

Stept 6

To display the percentages on the map, we need to create the centroids data. To accomplish this, we will first need to build a function called get.centroid.shp. Here is the implementation of the function:

get.centroid.shp=function(x){
    N <- length(x)  # Number of polygons
    # Initialise data.frame
    Centroids.shp <- data.frame(matrix(NA, N, 2, dimnames = list(NULL, c("long", "lat"))))
    for(i in 1:N){
        Centroids.shp[i,] <- x@polygons[[i]]@labpt
    }
    return(Centroids.shp)
}

The data with the centers is then:

dt.centers=cbind.data.frame(dep_paris@data,get.centroid.shp(dep_paris))
i=match(dt.centers$id_geofla,dep_paris_fr$id)
dt.centers$group=dep_paris_fr$id[i]

Step 7

Additionally, we will create another map that represents the Ile-de-France departments, with the ones that need to be zoomed in placed in the center of the map. For this purpose, we will create a separate polygon and centroids data specific to this map.

xx=unique(dep_paris$id_geofla)[c(3:6)]

i=which(dep_paris_fr$id%in%xx)
dep_paris_fr2=dep_paris_fr%>%slice(i)


dt.centers2=dt.centers[c(3:6),]
dt.centers2$long[2]=dt.centers2$long[2]-.05

Step 8

We are now prepared to create the maps. Let’s begin with the larger map that includes all departments, and then we will proceed to create the smaller, zoomed-in map. For better visualization, we will set the background to be transparent in both maps, allowing them to be overlapped seamlessly:

p2<-ggplot(dep_paris_fr, aes(x=long, y=lat, group=group)) +
  geom_polygon(aes(fill=NEET),color = "black")+
  labs(x="",y="")+ theme_bw()+
  scale_fill_gradientn(colours=c("lightblue","blue"),limits=c(12,28))+
  coord_fixed()
p2<-p2+theme(axis.line=element_blank(),
             axis.text.x=element_blank(),
             axis.text.y=element_blank(),
             axis.ticks=element_blank(),
             axis.title.x=element_blank(),
             axis.title.y=element_blank(),
             panel.grid.major = element_blank(),
             panel.grid.minor = element_blank(),
             panel.border = element_blank(),
             panel.background = element_rect(fill='transparent'),
             plot.background = element_rect(fill='transparent', color=NA))
p2<-p2+ theme(legend.position="none")
p2<-p2+ geom_label(data=dt.centers,aes(label = NEET_lb, x = long, y = lat),size=3)+
  labs(title = "NEET (%)")
p2

And now this the smallest map:

p2a<-ggplot(dep_paris_fr2, aes(x=long, y=lat, group=group)) +
  geom_polygon(aes(fill=NEET),color = "black")+
  labs(x="",y="")+ theme_bw()+
  scale_fill_gradientn(colours=c("lightblue","blue"),limits=c(12,28))+
  coord_fixed()
p2a<-p2a+theme(axis.line=element_blank(),
               axis.text.x=element_blank(),
               axis.text.y=element_blank(),
               axis.ticks=element_blank(),
               axis.title.x=element_blank(),
               axis.title.y=element_blank(),
               panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               panel.border = element_blank(),
               panel.background = element_rect(fill='transparent'),
               plot.background = element_rect(fill='transparent', color=NA))
p2a<-p2a+ theme(legend.position="none")
p2a<-p2a+ geom_label(data=dt.centers2,aes(label = NEET_lb_2, x = long, y = lat),size=2)

To combine the two maps together, we will utilize the patchwork package. This package enables us to arrange and merge multiple plots into a single composition. In our case, we will position the smaller map in the top-right corner of the larger map.

library(patchwork)
p2A<-p2 + inset_element(p2a, left = 0.7, bottom = 0.8, right = 1, top = 1.2)
p2A