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.
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 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 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.3 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.4 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 Familiarisation avec TDLM
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\)).
# Préparons les données d'entrées
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)Nous allons ensuite estimer les flux à l’aide d’une loi gravitaire
exponentielle (law = "NGravExp") avec un modèle doublement contraint
(model = "DCM").
# Loi gravitaire exponentielle avec un modèle doublement contraint
od_sim <- run_law_model(law = "NGravExp",
mass_origin = mi,
mass_destination = mj,
distance = distance_mtp,
opportunity = NULL,
param = 0.01,
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = FALSE,
nbrep = 3,
check_names = TRUE)Nous pouvons faire une comparaison rapide pour comprendre les résultats.
# Comparaison rapide des flux simulés et observés
od_sim$replication_1[1:10,1:10]
od_sim$replication_2[1:10,1:10]
od_sim$replication_3[1:10,1:10]
od_mtp[1:10,1:10]Pour aller plus loin, la fonction gof() (goodness-of-fit) permet de
calculer plusieurs indicateurs pour comparer les matrices simulées à la
matrice observée.
# Goodness-of-fit entre flux simulés et observés
gof(od_sim, od_mtp,
measures = "all",
distance = distance_mtp)Vous pouvez également comparer les distances moyennes simulées à la distance moyenne observée à l’aide de votre fonction ou du code ci-dessous.
d_mean_1 <- sum(od_sim$replication_1 * distance_mtp) / sum(od_sim$replication_1)
d_mean_1
d_mean_2 <- sum(od_sim$replication_2 * distance_mtp) / sum(od_sim$replication_1)
d_mean_2
d_mean_3 <- sum(od_sim$replication_3 * distance_mtp) / sum(od_sim$replication_1)
d_mean_3
d_mean_obs <- sum(od_mtp * distance_mtp) / sum(od_mtp)
d_mean_obs2.2.2 Calibration de paramètre
Nous avons observé que notre modèle tendait à surestimer la distance moyenne de navettage. Cela signifie que la valeur du paramètre \(\beta\), qui contrôle l’effet de la distance dans la loi gravitaire, est probablement trop faible.
Pour améliorer l’ajustement, nous allons tester plusieurs valeurs de \(\beta\) afin d’évaluer la qualité des ajustements via différents indicateurs de goodness-of-fit et de comparer la distance moyenne simulée avec celle observée.
# Loi gravitaire exponentielle avec un modèle doublement contraint
od_sim <- run_law_model(law = "NGravExp",
mass_origin = mi,
mass_destination = mj,
distance = distance_mtp,
opportunity = NULL,
param = c(0.01, 0.05, 0.1, 0.2), # Valeurs testées pour β
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = FALSE,
nbrep = 1,
check_names = TRUE)# Comparaison entre flux simulés et flux observés
gof(od_sim, od_mtp,
measures = "all",
distance = distance_mtp)# Calcul de la distance moyenne pour chaque paramètre β testé
d_mean_1 <- sum(od_sim$parameter_1$replication_1 * distance_mtp) /
sum(od_sim$parameter_1$replication_1)
d_mean_1
d_mean_2 <- sum(od_sim$parameter_2$replication_1 * distance_mtp) /
sum(od_sim$parameter_2$replication_1)
d_mean_2
d_mean_3 <- sum(od_sim$parameter_3$replication_1 * distance_mtp) /
sum(od_sim$parameter_3$replication_1)
d_mean_3
d_mean_4 <- sum(od_sim$parameter_4$replication_1 * distance_mtp) /
sum(od_sim$parameter_4$replication_1)
d_mean_4
# Distance observée
d_mean_obs <- sum(od_mtp * distance_mtp) / sum(od_mtp)
d_mean_obs2.2.3 Exercices
1. Trouver la valeur optimale du paramètre (selon le CPC), arrondie à trois chiffres après la virgule.
Solution (cliquer pour afficher)
# Loi gravitaire exponentielle avec un modèle doublement contraint
od_sim <- run_law_model(law = "NGravExp",
mass_origin = mi,
mass_destination = mj,
distance = distance_mtp,
opportunity = NULL,
param = c(0.124, 0.125, 0.126, 0.127, 0.128, 0.129),
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = FALSE,
nbrep = 1,
check_names = TRUE)
gof(od_sim, od_mtp,
measures = "all",
distance = distance_mtp)2. Utiliser la fonction extract_opportunities() pour calculer le nombre
d’opportunités entre les origines et destinations en vous basant sur la
population (i.e. \(m_i\)).
Solution (cliquer pour afficher)
3. Chercher le paramètre optimal (selon le CPC) de la loi radiative
(law = "RadExt").
Solution (cliquer pour afficher)
# Loi radiative avec un modèle doublement contraint
od_sim_rad <- run_law_model(law = "RadExt",
mass_origin = mi,
mass_destination = mj,
distance = NULL,
opportunity = opportunity_mtp,
param = c(0.126, 0.127, 0.128, 0.129),
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = FALSE,
nbrep = 1,
check_names = TRUE)
gof(od_sim_rad, od_mtp,
measures = "all",
distance = distance_mtp)4. Quelle est la meilleure loi sur ce cas d’étude : gravitaire ou radiative ?
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ées et simulées
sous un format “réseaux” (i.e. trois colonnes) à l’aide de la fonction
mat_to_net() du package bioregion.
# Conversion des matrices en format "réseau"
net_obs <- mat_to_net(od_mtp, weight=TRUE)
net_sim <- mat_to_net(od_sim$parameter_2$replication_1, weight=TRUE)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 et osmdata :
# 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.