TP2 Estimation des flux domicile-travail
Avant de commencer le TP, je vais vous présenter une introduction à la modélisation des flux domicile-travail à l’aide des modèles d’interactions spatiales, ainsi que le contenu du TP. Pour aller plus loin, ce TP et le package TDLM s’appuient sur les travaux présentés dans Lenormand et al. (2016).
2.1 Prise en main des données
Le package R TDLM met à disposition des données d’exemple de flux domicile-travail entre les 342 communes du département de l’Hérault en 2020.
- od_mtp : matrice Origine-Destination (OD) des flux observés.
- distance_mtp : matrice des distances.
- mass_mtp : masses et marges (population, flux sortants, flux entrants).
- mtp : données spatiales associées (nous aurons besoin de sf pour les manipuler).
Nous allons dans un premier temps calculer quelques statistiques de base pour explorer ces données.
Créez un script TP2.R et enregistrez-le dans votre répertoire de travail, que vous pouvez également nommer TP2.
2.1.2 Taille et aperçu des objets
dim(od_mtp) # Dimensions et aperçu de la matrice OD
od_mtp[1:10,1:10]
dim(distance_mtp) # Dimensions et aperçu de la matrice de distance
distance_mtp[1:10,1:10]
mass_mtp[1:10,] # Aperçu des masses
mtp # Aperçu des entités spatiales
class(mtp) # Objet sf
plot(mtp) # Représentation spatiale de mtp
plot(mtp[,2]) # Représentation spatiale de l'attribut Longitude2.1.3 Réseau et matrice Origine-Destination
Avant de commencer à analyser les données, prenons un instant pour aborder la notion de réseau dans le cadre de ce TP et la façon dont elle est prise en charge par R. Il existe deux principales façons de représenter un réseau dirigé et pondéré comme celui que nous étudions ici.
- Le format matriciel : une matrice Origine-Destination carrée où chaque
case représente la valeur du flux entre une origine (ligne) et une destination
(colonne). C’est le format que nous utilisons avec od_mtp.
- Le format en trois colonnes : un tableau (data frame) comportant une colonne pour l’origine, une pour la destination, et une pour le poids (flux).
Ce second format est souvent appelé edge list (liste d’arêtes). Il est très pratique pour l’analyse de réseaux, car il permet d’intégrer les flux dans d’autres outils spécialisé dans l’analyse de graphe (comme le package igraph (Csárdi et al., 2025) par exemple).
Particularité importante : les liens inexistants (ou les flux nuls) ne sont généralement pas représentés. Autrement dit, seules les paires origine–destination où le flux est strictement positif apparaissent.
Nous pouvons facilement convertir notre matrice od_mtp en un format
“trois colonnes” à l’aide de la fonction mat_to_net() du package
bioregion.
# Conversion matrice -> "réseau"
net_mtp <- mat_to_net(od_mtp, weight = TRUE)
net_mtp[1:10,]
# Nombre de paires "non nulles"
dim(net_mtp)Et inversement avec la fonction net_to_mat().
2.1.4 Statistiques de base
# Population moyenne
mean(mass_mtp$Population)
# Entité avec le plus d'habitants
mass_mtp[which.max(mass_mtp$Population), ]
# Flux totaux sortants et entrants
sum(mass_mtp$Out_Commuters)
sum(mass_mtp$In_Commuters)
# Comparaison
plot(mass_mtp$Population, mass_mtp$Out_Commuters,
xlab="Population", ylab="Flux sortants",
main="Population vs Flux sortants")
# Comparaison log-log
plot(mass_mtp$Population, mass_mtp$Out_Commuters,
log="xy",
xlab="Population", ylab="Flux sortants",
main="Population vs Flux sortants")
# Vérifier que les flux sont cohérents
sum(mass_mtp$Out_Commuters == apply(od_mtp, 1, sum))
sum(mass_mtp$In_Commuters == apply(od_mtp, 2, sum))2.1.5 Exercices
1. Calculer la distance moyenne de navettage à partir des matrices
od_mtp et distance_mtp. Calculer également la distance moyenne de
navettage par commune (à l’aide de la fonction apply()). Ces métriques vont
nous resservir. N’hésitez pas à développer des fonctions dédiées dans un script
et à les charger ensuite avec source(), comme vu lors du premier TP.
Solution (cliquer pour afficher)
2. Visualiser les flux sortants par commune sur une carte. Pour cela,
ajoutez les valeurs de Oi (flux sortants) à l’objet spatial mtp et
affichez-les avec la fonction plot().
Solution (cliquer pour afficher)
# Ajouter les flux sortants à l'objet spatial
mtp$Oi <- mass_mtp$Out_Commuters
# Visualiser avec une carte simple
plot(mtp["Oi"], main="Flux sortants par commune")
# Visualisation par quantile
library(viridis)
classes <- cut(mtp$Oi,
breaks = quantile(mtp$Oi, probs = seq(0,1,length=11),
na.rm=TRUE),
include.lowest = TRUE,
labels = FALSE)
colors <- viridis(10)
plot(mtp["Oi"], main="Flux sortants par commune", col=colors[classes])
legend("right",
legend = paste0("Q", 1:10),
bty = "n",
xpd = NA,
fill = colors)3. Visualiser le ratio flux entrant / flux sortant par commune sur une carte.
Solution (cliquer pour afficher)
2.2 Estimation des flux avec TDLM
2.2.1 Préparation de données
Nous allons commencer par préparer les données en formatant les masses (\(m_i\) et \(m_j\)) et les marges (\(D_i\) et \(D_j\)).
mi <- mass_mtp$Population
names(mi) <- rownames(mass_mtp)
mj <- mass_mtp$Population
names(mj) <- rownames(mass_mtp)
Oi <- mass_mtp$Out_Commuters
names(Oi) <- rownames(mass_mtp)
Dj <- mass_mtp$In_Commuters
names(Dj) <- rownames(mass_mtp)mi et mj vont nous permettre d’estimer les probabilités de déplacement à partir d’une loi (de gravité ou de radiation), tandis que Oi et Dj vont nous permettre de générer les flux à partir de cette loi à l’aide d’un modèle doublement contraint.
La loi de gravité ou loi gravitaire se base sur la distance entre deux communes (distance_mtp) pour estimer le coût de déplacement, contrairement à la loi de radiation, qui se base sur le nombre d’opportunités entre celles-ci (opportunity_mtp, que nous calculerons plus tard). Dans les deux cas, il y a un paramètre à calibrer pour ajuster l’effet du coût.
Dans ce TP, les paramètres ont déjà été calibrés, mais vous pourrez tester différentes valeurs si vous êtes en avance.
Le modèle doublement contraint a été choisi pour comparer les lois de manière transparente. Ce modèle permet de générer les flux à partir d’une loi en respectant les marges (Oi et Dj).
Il s’agit du modèle le plus performant car il contraint les deux types de marge en même temps, mais n’hésitez pas à tester d’autres modèles si vous êtes en avance.
Il est important de noter que lee processus de génération des flux est stochastique (i.e. tirage aléatoire).
Dans ce TP, nous générerons qu’une seule matrice mais vous pouvez comparer plusieurs tirages aléatoires si vous êtes en avance.
2.2.2 La loi gravitaire
Nous allons ici utiliser la loi gravitaire
exponentielle (law = "NGravExp") avec un modèle doublement contraint
(model = "DCM").
Il existe plusieurs types de lois gravitaires implémentées dans TDLM, nous ne rentrerons pas dans le détail ici mais n’hésitez pas à jouer avec les différentes lois et modèles.
La fonction run_law_model() ci-dessous permet de générer une matrice
od_grav.
od_grav <- run_law_model(law = "NGravExp",
mass_origin = mi,
mass_destination = mj,
distance = distance_mtp,
opportunity = NULL,
param = 0.126,
model = "DCM",
out_trips = Oi,
in_trips = Dj,
nbrep = 1,
check_names = TRUE)
od_grav <- od_grav$replication_1Nous pouvons faire une comparaison rapide entre la matrice observée od_mtp et la matrice simulée od_grav pour comprendre les résultats.
Pour aller plus loin, la fonction gof() (goodness-of-fit) permet de
calculer plusieurs indicateurs pour comparer la matrice simulée à la
matrice observée. Nous nous concentrerons ici sur le proportion de navetteurs
en commun (measures = "CPC"), calculant comme son nom l’indique la proportion
d’information en commun entre la matrice observée et la matrice simulée. Cet
indicateur est compris entre 0 et 1. Une valeur de 0 signifie qu’il n’y a aucun
navetteur en commun entre les deux matrices. Une valeur de 1 siginifie que les
deux matrices sont identiques.
Vous pouvez également comparer les distances moyennes simulées à la distance moyenne observée à l’aide de votre fonction ou du code ci-dessous.
2.2.3 La loi de radiation
Nous allons utiliser maintenant la loi de radiation étendue
(law = "RadExt") avec un modèle doublement contraint (model = "DCM").
Il existe plusieurs types de lois de radiation implémentées dans TDLM, nous ne rentrerons pas dans le détail ici mais n’hésitez pas à jouer avec les différentes lois et modèles.
Nous allons avoir besoin ici du nombre d’opportunités entre les origines et
destinations pour estimer la matrice de coût utilisée dans les lois de radiation.
Nous nous baserons sur la population (i.e. \(m_i\)) en utilisant la fonction
extract_opportunities().
opportunity_mtp <- extract_opportunities(opportunity = mi,
distance = distance_mtp,
check_names = TRUE)La fonction run_law_model() ci-dessous permet de générer une matrice
od_rad.
od_rad <- run_law_model(law = "RadExt",
mass_origin = mi,
mass_destination = mj,
distance = NULL,
opportunity = opportunity_mtp,
param = 0.127,
model = "DCM",
out_trips = Oi,
in_trips = Dj,
nbrep = 1,
check_names = TRUE)
od_rad <- od_rad$replication_1Nous pouvons faire une comparaison rapide pour comprendre les résultats.
Ici aussi, vous pouvez utiliser la fonction gof() (goodness-of-fit) pour
comparer la matrice simulée et la matrice observée.
Vous pouvez également comparer les distances moyennes simulées à la distance moyenne observée à l’aide de votre fonction ou du code ci-dessous.
2.2.4 Exercices
1. Utilisez la fonction cor() pour calculer le coefficient de corrélation
de Pearson entre les distances et nombres d’opportunités entre origines et
destination. Visualisez le nuage de points avec la fonction plot()
Solution (cliquer pour afficher)
2. Utiliser la fonction run_law_model() avec la loi uniforme
(law = "Unif") pour générer une matrice de référence.
Solution (cliquer pour afficher)
3. Utiliser la fonction gof() avec une list() pour comparer le CPC
entre od_grav, od_rad et od_unif (calculer précédemment) et
od_mtp.
Solution (cliquer pour afficher)
4. Quelle est la meilleure loi sur ce cas d’étude : gravitaire ou radiation ?
2.3 Flux et itinéraires
Nous allons maintenant rentrer un peu plus dans le détail de l’analyse des flux et je vais vous proposer une façon simple de calculer rapidement des itinéraires.
2.3.1 Identifier les flux
Commençons par mettre nos matrices origine-destination observée et simulée
(la meilleure des trois) sous un format “réseaux” (i.e. trois colonnes) à l’aide
de la fonction mat_to_net() du package bioregion.
Nous pouvons ensuite les trier par ordre décroissant de taille de flux et comparer le top 10 des plus gros flux dans les deux matrices.
# Tri des flux par ordre décroissant
net_obs <- net_obs[order(net_obs$Weight, decreasing=TRUE),]
net_sim <- net_sim[order(net_sim$Weight, decreasing=TRUE),]
# Top 10 des flux observés et simulés
net_obs[1:10,]
net_sim[1:10,]Nous observons que le plus gros flux d’individus va de la commune ayant pour code INSEE 34057 (Castelnau-le-Lez) vers le code INSEE 34172 (Montpellier). Nous allons maintenant essayer de déterminer un itinéraire permettant d’aller d’une commune à l’autre.
2.3.2 Calcul d’itinéraires
Pour le calcul d’itinéraire, nous allons avoir besoin d’un réseau routier autour de Montpellier. Vous pouvez télécharger ce réseau ici.
Il s’agit d’un objet road stocké dans un fichier road.rda que vous
pouvez placer dans votre répertoire de travail et importer avec la commande
load() comme vu lors du premier TP.
Si cela vous intéresse, voici le bout de code que j’ai utilisé pour obtenir road à l’aide des packages dodgr (Padgham, M., 2019) et osmdata (Padgham, M., 2019):
# Charger les packages
library(dodgr)
library(osmdata)
# Charger, filtrer, pondérer et simplifier le réseau routier
road_dl <- dodgr_streetnet(getbb("Montpellier, France"), quiet = TRUE)
road <- road_dl[road_dl$highway %in% c("motorway", "trunk", "primary",
"secondary", "tertiary", "residential"),]
road <- weight_streetnet(road, wt_profile = "motorcar")
road <- dodgr_contract_graph(road)
# Exporter road
save(road, file="road.rda")Pour calculer l’itinéraire entre les centroides des deux communes, nous allons avoir besoin du package dodgr et du package ggplot2 pour visualiser facilement cet itinéraire.
# Charger les packages
library(dodgr)
library(ggplot2)
# Charger le réseau routier
load("road.rda")
# Définir l'origine et la destination (coordonnées des centroides)
from <- as.numeric(mtp[mtp$ID == 34057, 2:3])[-3] # Castelnau-le-Lez
to <- as.numeric(mtp[mtp$ID == 34172, 2:3])[-3] # Montpellier
# Trouver les sommets les plus proches
vertices <- dodgr_vertices(road)
from_id <- vertices$id[which.min((vertices$x - from[1])^2 + (vertices$y - from[2])^2)]
to_id <- vertices$id[which.min((vertices$x - to[1])^2 + (vertices$y - to[2])^2)]
# Calculer le chemin le plus court
path_nodes <- dodgr_paths(road, from = from_id, to = to_id, vertices = TRUE)
# Convertir en sf LINESTRING
coords <- do.call(rbind, lapply(path_nodes[[1]][[1]], function(id) {
c(vertices$x[vertices$id == id], vertices$y[vertices$id == id])
}))
path_line <- st_sfc(st_linestring(coords), crs = 4326)
# Tracer le chemin
ggplot() +
geom_sf(data = mtp[mtp$ID==34057 | mtp$ID==34172,], color = "grey80") +
geom_sf(data = path_line, color = "#CC6666", size = 1.2) +
ggtitle("Castelnau-le-Lez → Montpellier")2.3.3 Exercice
Essayez de bien vous approprier le code ci-dessus et calculez un itinéraire pour le second plus gros flux. Vous pouvez ensuite représenter les deux itinéraires (le plus gros et le second plus gros flux) sur une même carte afin de les comparer visuellement.
Solution (cliquer pour afficher)
# Définir l'origine et la destination (coordonnées des centroides)
from <- as.numeric(mtp[mtp$ID == 34129, 2:3])[-3] # Lattes
to <- as.numeric(mtp[mtp$ID == 34172, 2:3])[-3] # Montpellier
# Trouver les sommets les plus proches
vertices <- dodgr_vertices(road)
from_id <- vertices$id[which.min((vertices$x - from[1])^2 + (vertices$y - from[2])^2)]
to_id <- vertices$id[which.min((vertices$x - to[1])^2 + (vertices$y - to[2])^2)]
# Calculer le chemin le plus court
path_nodes <- dodgr_paths(road, from = from_id, to = to_id, vertices = TRUE)
# Convertir en sf LINESTRING
coords <- do.call(rbind, lapply(path_nodes[[1]][[1]], function(id) {
c(vertices$x[vertices$id == id], vertices$y[vertices$id == id])
}))
path_line_2 <- st_sfc(st_linestring(coords), crs = 4326)
# Tracer les chemins
ggplot() +
geom_sf(data = mtp[mtp$ID==34057 | mtp$ID==34129 | mtp$ID==34172,], color = "grey80") +
geom_sf(data = path_line, color = "#CC6666", size = 1.2) +
geom_sf(data = path_line_2, color = "steelblue2", size = 1.2)2.4 Application R Shiny
Pour finir, je vais vous montrer un exemple d’application plus élaborée que celle du TP1 que nous pouvons développer à partir de l’objet spatial mtp.
Les fichiers sources de cette application sont disponibles ici. Vous pouvez dézipper le dossier dans votre répertoire de travail.
Dans le TP1, nous avions construit une application dans un seul fichier app.R contenant à la fois l’interface utilisateur (ui) et la logique du serveur (server). Dans ce TP, l’application est découpée en plusieurs fichiers, ce qui rend le code plus structuré et plus facile à maintenir lorsque l’application devient plus complexe.
L’application est composée de quatre fichiers principaux :
global.R contient les éléments communs à l’application, chargés une seule fois au lancement. On y met par exemple le chargement des packages, des données ou la définition de fonctions utiles.
ui.R définit l’interface utilisateur (User Interface). C’est ici que l’on décrit la mise en page, les onglets, les boutons, les cartes ou graphiques que l’utilisateur pourra voir et manipuler.
server.R Contient la logique de l’application, c’est-à-dire la manière dont les entrées de l’utilisateur sont traitées et reliées aux sorties affichées (graphiques, cartes, tableaux…).
style.css est un fichier de mise en forme optionnel. Il permet de modifier l’apparence de l’application (couleurs, polices, marges, tailles, etc.).
Prenez le temps de parcourir ces fichiers et d’expérimenter en changeant le titre de l’application, en renommant l’onglet du tableau de bord, en modifiant la palette de couleurs de la cartes ou encore personnaliser les couleurs dans style.css.
Ces ajustements simples vous aideront à mieux comprendre l’articulation des différents fichiers et à prendre en main la structure d’une application Shiny plus avancée.
2.4.1 Exercice
Si vous êtes en avance essayez d’ajouter la variable Ratio dans la table attributaire et dans le menu déroulant et afficher les valeurs associées sur la carte.