Ce document présente un usage possible des packages R qui interfacent OpenStreetMap (OSM): géocoder une adresse, import de tuiles et de couches géographiques, import de points d’intérêt, routage et calcul d’itinéraires.
Il propose aussi des méthodes spatiales (carroyages, KDE, segments de lignes) et des visualisations associées pour restituer graphiquement la structure spatiale de semis de points ou de lignes.
En utilisant exclusivement des informations issues d’OSM, l’exercice consiste à identifier les secteurs accessibles à pied dans le voisinage géographique d’un lieu défini par une adresse.
Reproduire l’analyse
Pour reproduire l’analyse, télécharger l’archive ZIP associée au dépôt https://github.com/riatecom/osm-elementr-2025-application, ouvrir le fichier proj.Rproj puis le fichier script.R qui reproduit en tout point le contenu du support de formation. N’oubliez pas d’installer les packages en amont via la commande inscrite dans le README.
Cette analyse est reproductible en n’importe quel lieu, avec une adresse (adr) et un nom de ville (lib). Néanmoins les temps de calcul utilisant le serveur de démonstration d’OSRM ont été paramétrés pour ne pas le surcharger. L’exécution de ces blocs de code prennent du temps.
En conséquence, 6 jeux de données ont été préparés en amont à partir de l’adresse des gares SNCF de Lille, Toulouse, Grenoble, Nantes, Clermont-Ferrand et Compiègne, ainsi que de la station de métro Front Populaire. Les blocs de code relevant des calculs d’itinéraires / temps de trajets sont placés en commentaires dans le fichier script.R, mais sont fonctionnels.
lib adr
1 Grenoble 1 Place de la Gare, 38000 Grenoble, France
2 Clermont 40 Avenue de l'Union Soviétique, 63000 Clermont-Ferrand, France
3 Toulouse Gare SNCF Matabiau, 31500 Toulouse, France
4 Compiegne Gare de Compiègne, 60200 Compiègne, France
5 Condorcet Front Populaire, 93210 Saint-Denis, France
6 Nantes 27 Boulevard de Stalingrad, 44000 Nantes, France
7 Lille Lille Flandres, 59000 Lille, France
Nous présentons ici les résultats depuis la gare SNCF de Compiègne.
# 1. Packages qui interfacent des données OSM#install.packages('maposm', repos = 'https://riatelab.r-universe.dev')library(tidygeocoder) # géocodagelibrary(maptiles) # import de tuiles OSM (raster)library(osmdata) # import de données OSM (vecteur)library(maposm) # import de données OSM (couches géo)library(osrm) # calcul d'itinéraires# 2. Autres packages utilitaireslibrary(sf) # manipulation de données vectorielleslibrary(terra) # manipulation de données rasterlibrary(mapsf) # cartographie thématiquelibrary(mapiso) # discrétiser des données continueslibrary(maplegend) # légendeslibrary(spatstat) # analyse de semis de pointslibrary(mapview) # cartographie interactivelibrary(stplanr) # segmentiser des lignes
1. Géocodage
Le package tidygeocoder(Cambon et al. 2025) permet d’utiliser des services de géocodage en ligne. Par défaut, le package utilise Nominatim qui se base sur les données OSM.
Passing 1 address to the Nominatim single address geocoder
Query completed in: 1 seconds
pt
Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 2.82352 ymin: 49.42206 xmax: 2.82352 ymax: 49.42206
Geodetic CRS: WGS 84
# A tibble: 1 × 2
address geometry
* <chr> <POINT [°]>
1 Gare de Compiègne, 60200 Compiègne, France (2.82352 49.42206)
2. Import de tuiles OSM
Le package maptiles(Giraud 2025) permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au format raster.20 niveaux de zoom sont disponibles, selon la taille de l’objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.
Plusieurs modèles de tuiles sont disponibles ici, certains nécessitent une inscription.
emprise <-st_buffer(pt, 3000) |>st_bbox()tiles <-get_tiles(x = emprise,project =FALSE, # crs en output : 3857crop =TRUE,cachedir ="cache",zoom =14)pt_reproj <-st_transform(pt, crs ="EPSG:3857")mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt =paste0("Gare de ", lib))mf_title("Tuiles OSM - Défaut")
# Autre exemple de fournisseurtiles <-get_tiles(x = emprise,project =FALSE,crop =TRUE,provider ="OpenStreetMap.HOT",zoom =14)mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt =paste0("Gare de ", lib))mf_title("Tuiles OSM - OpenStreetmap.HOT")
3. Import de couches géographiques OSM
Pour la cartographie, le package maposm (pas sur le CRAN) se base sur les fonctions du package osmdata(Padgham et al. 2023) (partie suivante) pour télécharger les données OSM, au format vectoriel.
La fonction principale, om_get(), permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc) autour d’un point, dans un rayon donné. Le résultat prend la forme d’une liste d’objets sf.
res <-om_get(x =c(st_coordinates(pt)[1],st_coordinates(pt)[2]),r =2000)
mf_map(res$zone, col ="#f2efe9", border =NA, add =FALSE)mf_map(res$urban, col ="#e0dfdf", border ="#e0dfdf", lwd = .5, add =TRUE)mf_map(res$green, col ="#c8facc", border ="#c8facc", lwd = .5, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$building, col ="#d9d0c9", border ="#c6bab1", lwd = .5, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Autour de la gare")mf_arrow()mf_scale(size =500, scale_units ="m")mf_credits("OpenStreetMap contributors, 2025", pos ="bottomleft")
4. Import de tags avec osmdata
Le package osmdata permet d’extraire des données OSM à partir de l’API Overpass. Il nécessite dans un premier temps d’établir une bounding box au sein de laquelle seront extraites les données, puis de l’exécuter (fonction add_osm_feature()). Le résultat est une liste, qu’il faut ensuite convertir en objet sf.
Les objets OSM prennent la forme de nodes (points), ways (lignes ou polygones), ou relations (relation entre plusieurs objets). Ces éléments sont décrits par un ensemble de tags, recensées sous forme de clé-valeurs (key-values). Ici, nous recherchons dans la clé amenity les points (nodes) correspondant aux restaurants, bars et cafés (value = restaurant, bar, cafe). En dehors des clés-valeurs principales, d’autres attributs peuvent être présents, comme en l’occurence le type de cuisine, ou bien les horaires.
# Définition d'une bounding box (emprise Compiègne)q <-opq(bbox = emprise, osm_types ="node") # seulement les points# Extraction des restaurants, bars et cafésreq <-add_osm_feature(opq = q,key ='amenity',value =c("bar", "cafe", "restaurant"))poi <-osmdata_sf(req)poi <- poi$osm_points# Sélectionner les variablespoi <- poi[, c("osm_id", "name", "amenity", "cuisine")]# Intersection et reprojectionpoi <-st_intersection(poi, st_transform(res$zone, crs =4326))poi_reproj <-st_transform(poi , crs ="EPSG:3857")
La fonction om_map() affiche les couches géographiques récupérées.
# Affichageom_map(res, title ="Autour de la gare")mf_map(poi_reproj, type ="typo", var ="amenity", cex = .7, pch =15, leg_pos ="topright", leg_val_cex = .9, leg_title_cex = .9, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, lwd =2, col ="darkblue", add =TRUE)
Il est possible d’utiliser le package mapview(Appelhans et al. 2023) pour afficher les points ainsi que leurs attributs, de manière interactive. Dans le cas de données issues d’OSM, il est intéressant d’observer le nombre d’attributs disponibles par point.
mapview(res$zone, alpha.regions =0, layer.name ="Zone") +mapview(poi, layer.name ="Manger et boire")
osmextract
Le package osmextract(Gilardi et Lovelace 2024) permet d’extraire des données directement depuis une base de données OSM, sans avoir à passer par une API. Le plus petit niveau géographique auquel il est possible de télécharger les données correspond en France aux anciennes régions. Pour travailler sur Compiègne, il nous faudrait d’abord télécharger l’ensemble des données OSM de la Picardie, disponibles ici, et pesant plus de 240 Mb.
Le code ci-dessous effectue les mêmes manipulations que précédemment.
Une fois les points d’intérêt récupérés (aussi appelés POI, pour Points of Interest), il est possible de les agréger dans une maille.
# Création d'une grille régulière hexagonale de 300m de résolutiongrid <-st_make_grid(res$zone, cellsize =300, square =FALSE)grid <-st_intersection(grid, res$zone)grid <-st_sf(ID =1:length(grid), geom = grid)# Compter les POIs dans les carreauxinter <-st_intersects(grid, poi_reproj, sparse =TRUE)grid$n_poi <-lengths(inter)# Cartographieom_map(res)mf_map(grid, var ="n_poi", type ="choro", breaks =c(0,1,3,5,10,15),alpha = .6, border =NA, leg_title ="Nombre \nd'aménités",leg_val_rnd =0, leg_pos ="topright", add =TRUE)mf_title("Aménités gustatives - données carroyées")
6. Lissage (KDE)
Le package spatstat(Baddeley, Turner, et Rubak 2025) regroupe de nombreuses fonctions pour analyser des semis de points. Il génère des objets de classe ppp(planar point pattern), dans des fenêtres d’observation, de classe owin(observation window).
La fonction density.ppp() calcule des valeurs à partir d’une estimation par noyau (KDE). Elle prend en entrée une portée (sigma), ainsi qu’une résolution au sein de laquelle sera calculée la densité, en pixels (eps).
p <-as.ppp(X =st_coordinates(poi_reproj), W =as.owin(res$zone))ds <-density.ppp(x = p, sigma =100, eps =10, positive =TRUE)
Une fois les densités calculées, nous créons un objet de type SpatRaster. La densité étant calculée sur un mètre carré, nous la multiplions par 10.000 afin de récupérer une densité par hectare.
# Calcul densité de restaurants par hectarer <-rast(ds) *100*100crs(r) <-st_crs(poi_reproj)$wktr
class : SpatRaster
dimensions : 400, 400, 1 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : 312312.9, 316312.9, 6344781, 6348781 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
source(s) : memory
name : lyr.1
min value : 2.225074e-304
max value : 1.968529e+00
pal <-c( "#418e40", "#a5be6a", "#e2cf6a", "#c19849", "#ab6c32", "#842f1a")mf_map(res$zone, col ="#f2efe9", border =NA, expandBB =c(0,.3,0,0))mf_raster(r, leg_title ="N. restaurants / ha.", leg_pos ="left", pal = pal, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(poi_reproj, pch =20, cex = .3, col ="red", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("La montagne culinaire - Compiègne")
Le package mapiso(Giraud 2023) permet de transformer des grilles régulières (ici notre raster) en contours de polygones (format vectoriel).
# Discrétiser le rasterbks <- bks <-c(0, 0.1, 0.5, 1, minmax(r)[2])isos <-mapiso(x = r, var ="lyr.1 ", breaks = bks)mf_map(res$zone, col ="#f2efe9", border =NA)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5,add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(isos, "isomin", "choro", border =NA, alpha =0.5, breaks = bks,leg_title ="Densité\n (n/ha)", add =TRUE)mf_map(poi_reproj, pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Densités de restaurants - vectorisé")
7. Temps d’accès
Le package osrm(Giraud 2024) interface l’engin de routage OSRM et propose différentes fonctions pour le calcul de temps de trajet et d’itinéraires. Ces calculs s’effectuent sur la base de profils .lua qui pour chaque mode de transport (voiture, vélo, à pied) définissent des algorithmes permettant la sélection de telle ou telle route, une vitesse par défaut selon le type de route, etc. Le profil piéton est disponible ici.
D’après la documentation de l’API OSRM, le nombre de requêtes sur le serveur de démo est limité à 512 par connection à l’API, espacées d’au moins 5 secondes. Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d’installer une instance OSRM sur son propre serveur (voir ici).
7.1. Depuis la gare
Pour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la gare de Compiègne et l’ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction osrmTable(), qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d’origine et de destination, et d’un vecteur correspondant à la métrique calculée. Il est possible de calculer les distances et temps de trajet entre plusieurs couples origine / destination en une seule requête.
Nous construisons d’abord une matrice O/D composée de coordonnées X/Y des points d’origine et de destination (CRS 4326).
Utilisation de l’API OSRM
Par souci de respect d’utilisation de l’API, une partie du code ci-dessous est passé en commentaire et non joué. Le résultat est enregistré sous format RDS pour être utilisé par la suite. Dans les morceaux de code suivants, le code apparaît dans des chunks sous option eval = FALSE.
# Création de la matrice O/Do <-data.frame(X =st_coordinates(pt)[, 1],Y =st_coordinates(pt)[, 2])row.names(o) <-1:nrow(o)poi <- poi[poi$amenity =="restaurant",]d <-data.frame(X =st_coordinates(poi)[, 1],Y =st_coordinates(poi)[, 2])# 1 requête : 1 point d'origine, 50 points de destination# mat_dist <- osrmTable(src = o,# dst = d,# measure = "distance",# osrm.profile = 'foot')# saveRDS(mat_dist, "data/mat_dist.rds")# Intervalle de 5 secondes entre les requêtes# Sys.sleep(5)# mat_dur <- osrmTable(src = o,# dst = d,# measure = "duration",# osrm.profile = 'foot')# saveRDS(mat_dur, "data/mat_dur.rds")mat_dist <-readRDS("data/mat_dist.rds")mat_dur <-readRDS("data/mat_dur.rds")dist <-t(mat_dist$distances)dur <-t(mat_dur$durations)poi$dist <- distpoi$dur <- durapply(poi$dur, 2, summary)
1
Min. 1.100
1st Qu. 7.000
Median 9.700
Mean 9.222
3rd Qu. 11.675
Max. 17.600
apply(poi$dist, 2, summary)
1
Min. 82.00
1st Qu. 522.50
Median 738.50
Mean 696.26
3rd Qu. 878.75
Max. 1501.00
Il est possible d’atteindre 20 restaurants, cafés ou bars en moins de 10 minutes à pied depuis la gare de Compiègne, et les trois quarts se situent à moins de 880 mètres à pied.
Les restaurants à moins de 5 minutes à pied de la gare sont les suivants :
poi$name[poi$dur <5]
[1] "So Gusto" "Bistrot de Flandre" "Brasserie du Nord"
[4] "L'Entre Mets" "Le Mood" "Taj Mahal"
[7] "Pizza Time's" "Mosaïque Pizza"
La fonction osrmIsochrone() du package osrm permet de calculer des isochrones depuis un point.
isos <-st_read("data/mat.gpkg", layer ="isos", quiet =TRUE)isos <-st_transform(isos, crs ="EPSG:3857")isos <-st_intersection(isos, res$zone)om_map(res, title ="Isochrones autour de la gare de Compiègne")mf_map(isos, type ="typo", var ="isomax", alpha = .6, pal ="Geyser",border ="white", leg_pos =NA, add =TRUE)mf_map(st_transform(poi, crs ="EPSG:3857"), var ="dur", type ="choro", breaks =sort(unique(c(isos$isomin, isos$isomax))),pal ="Geyser", pch =21, border ="white", leg_pos =NA,add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, lwd =2, col ="darkblue", add =TRUE)leg(type ="choro", val =sort(unique(c(isos$isomin, isos$isomax))), pal ="Geyser", title ="Temps\n(minutes)", pos ="left")
7.2. Depuis la grille
Afin de construire des indicateurs d’accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville à pied depuis chaque point de grille.
nrow(grid) *nrow(d)
[1] 9450
Nous avons en tout 9450 couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos 189 points d’origine et un point de destination, sous forme de boucle.
Le code pour effectuer les calculs est présenté ci-dessous, mais a été joué en local. Les résultats sont enregistrés dans le dossier data.
# Coordonnées des points d'origine : centroïdes des cellules de la grilleo <-st_centroid(grid) |>st_transform(crs =4326)o$X <-st_coordinates(o)[, 1]o$Y <-st_coordinates(o)[, 2]o <-st_set_geometry(o[, c("X", "Y")], NULL)
osrmNearest()
Pour que l’engin de routage calcule un itinéraire, les points d’origine et de destination doivent être localisés sur une route. Si le point d’entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l’itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d’espace (urbain ou rural).
La carte ci-dessous est un exemple d’indicateur d’accessibilité basé sur le temps de trajet : il s’agit du temps à parcourir pour atteindre le service le plus proche. Cet indicateur, que l’on peut décliner pour calculer le deuxième voire le troisième service le plus proche, et ainsi calculer la différence entre le premier et le second service, ouvre la question du choix de la destination : un individu est-il prêt à parcourir une plus grande distance pour bénéficier d’un meilleur / différent service ?
mat <-read.csv("data/matfoot.csv")# Restaurant le plus proche pour chaque point de grillepoi_min <-aggregate(dur_foot ~ ID_ORIG, mat, FUN = min, na.rm =TRUE)# Récupérer l'identifiant de destination associé# (Important pour la question du choix de la destination)poi_min <-do.call(rbind,lapply(split(mat, mat$ID_ORIG),function(x) x[which.min(x$dur_foot), ]))# Cartographierpoi_min_sf <-merge(grid, poi_min,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)om_map(res, title ="Temps d'accès au restaurant le plus proche")mf_map(poi_min_sf, var ="dur_foot", type ="choro", breaks =seq(0,15,3),alpha = .6, border =NA, leg_title ="Temps\n(minutes)",col_na ="#8B000060", leg_no_data ="> 15",leg_val_rnd =0, pal ="Greens", leg_pos ="topright", add =TRUE)
24 values are outside the class limits and will be classified as 'NA'.
Valhalla
Valhalla est un engin de routage open source basé sur les données OSM, au même titre que OSRM. Il fonctionne sur un système de tuiles contenant différents niveaux de routes, ainsi que les graphes routiers associés (les tuiles de niveau 0 contiennent les autoroutes, celles de niveau 2 les routes départementales). Les routes sont calculées au moment de la requête (tandis qu’elles sont pré-calculées sur OSRM) via différents algorithmes (A* pour les plus courts chemins, Contraction Hierarchies pour les longues distances). Valhalla propose de nombreux profils (vélo de location, multimodalité, taxi), ainsi que la possibilité de prendre en compte les pentes.
Un autre exemple d’indicateur d’accessibilité : les opportunités cumulées : combien de restaurants sont disponibles sous un certain seuil de temps ? Parmi eux, combien de restaurants offrant une cuisine spécifique ?
# moins de 5 minutesmat5 <- mat[mat$dur_foot <5 ,]mat5$n <-1mat5 <-aggregate(list(n_resto = mat5$n),list(ID_ORIG = mat5$ID_ORIG),FUN = sum)grid <-merge(grid[, "ID"], mat5,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)grid$n_resto[is.na(grid$n_resto)] <-0om_map(res, title ="Restaurants à moins de 5 minutes à pied")mf_map(grid, var ="n_resto", type ="choro", breaks =c(0,1,5,10,20,max(grid$n_resto)),alpha = .6, border =NA, leg_title ="Nombre de\nrestaurants",leg_val_rnd =0, pal ="Reds", leg_pos ="topright", add =TRUE)
8. Itinéraires
Pour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l’aide de la fonction osrmRoute(). Attention, l’API route d’OSRM permet de calculer des itinéraires seulement entre deux points. Une requête correspond donc à un couple origine / destination.
Le package stplanr(Lovelace, Ellison, et Morgan 2024) permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l’utilisons ici à des fins cartographiques, pour découper des lignes en segments.
routes <-st_read("data/mat.gpkg", layer ="routes", quiet =TRUE)routes <-st_transform(routes, crs ="EPSG:3857")routes$n <-1routes$distance <- routes$distance *1000routes <- routes[order(routes$distance, decreasing =TRUE), ]l_seg <-list()# La fonction line_segment découpe un objet sf LINESTRING en segments réguliersfor(i in1:nrow(routes)){ l_seg[[i]] <-line_segment( routes[i, ],segment_length =max(routes$distance) /20)}pal <-mf_get_pal(nrow(l_seg[[1]]), rev =TRUE, palette ="Blues 3")om_map(res, title ="Routes", theme ="dark")for(i innrow(routes):1){mf_map(l_seg[[i]], col = pal, lwd =1.6, add =TRUE)}leg(type ="cont", val =c(min(routes$distance), 1000, max(routes$distance)),pal = pal, pos ="topright", title ="Distance\n(mètres)")
# La fonction overline permet de compter le nombre d'occurence des tronçonsseg <-overline(routes, attrib ="n", fun = sum)om_map(res, title ="Routes")mf_map(seg, var ="n", type ="grad", col="darkred", breaks ="fisher", nbreaks =5, lwd =c(1,3,5,9,13), leg_val_rnd =0,leg_pos ="topright", leg_title ='Fréquence', add =TRUE)
Bonus : la tournée des pizzerias
Je fais une escale d’une demie-heure à Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps depuis la gare. Avec la fonction osrmTrip(), c’est possible !
Cette fonction répond au problème du voyageur de commerce, qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.
# Coordonnées des pizzeriapizz <- poi[!is.na(poi$cuisine) & poi$cuisine =="pizza" , "geometry"]# Ajouter le point de départ : la garepizz <-rbind(pt[1, "geometry"], pizz)pizz$id <-1:nrow(pizz)trip <-osrmTrip(loc = pizz,overview ="full",osrm.profile ="foot")trip <- trip[[1]]$tripst_write(trip, "data/mat.gpkg", layer ="trip", delete_layer =TRUE)
trip <-st_read("data/mat.gpkg", layer ="trip", quiet =TRUE)trips <-st_transform(trip, crs ="EPSG:3857")# Pour les labels : récupérer le début de la lignestart <- lwgeom::st_startpoint(trips)start <-do.call(rbind, start)start <-data.frame(id =rownames(trips), X = start[,1], Y = start[,2])start <-st_as_sf(start, coords =c("X", "Y"), crs ="EPSG:3857")om_map(res, theme ="pizza", title ="La route de la pizza")mf_map(trips, col ="black", lwd =3, lty =3, add =TRUE)mf_map(start[start$id !=1, ], pch =20, cex =2, col ="black", add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_label(start, var ="id", cex =1.2, overlap =TRUE, pos =4, col ="black", halo =TRUE)
sum(trip$duration)
[1] 23.86833
sum(trip$distance)
[1] 1.7733
D’après les données OSM, moins de 25 minutes suffisent pour faire la tournée des pizzas à Compiègne, le tout en moins de 2 kilomètres. On est large !!
Cambon, Jesse, Diego Hernangómez, Christopher Belanger, Daniel Possenriede, et Otto Hansen. 2025. tidygeocoder: An intuitive interface for getting data from geocoding services.https://CRAN.R-project.org/package=tidygeocoder.
Padgham, Mark, Bob Rudis, Robin Lovelace, Maëlle Salmon, et Joan Maspons. 2023. osmdata: Import ’OpenStreetMap’ Data as Simple Features or Spatial Objects. https://CRAN.R-project.org/package=osmdata.
Code source
---title: "Utiliser OpenStreetMap avec R"subtitle: "Mise en pratique"author: - name: Louis Laurian orcid: 0009-0002-7275-5536 affiliations: - ref: riate - name: Ronan Ysebaert orcid: 0000-0002-7344-5911 - name: Timothée Giraud orcid: 0000-0002-1932-3323affiliations: - id: riate name: Centre pour l’analyse spatiale et la géovisualisation (UAR RIATE) city: Paris state: France url: https://riate.cnrs.fr/date: 13 05 2025date-format: "D MMM YYYY"lang: frformat: html: embed-resources: true smooth-scroll: true fontsize: 0.9em code-tools: true code-fold: false toc: true toc-title: Utiliser OpenStreetMap avec R - Mise en pratique toc-depth: 3 toc-location: right bibliography: bib.bib linkcolor: "#d31329" css: custom.css include-after-body: "footer.html"editor_options: chunk_output_type: console---## ObjectifsCe document présente un usage possible des packages R qui interfacent OpenStreetMap (OSM): **géocoder une adresse**, **import de tuiles et de couches géographiques**, **import de points d'intérêt**, **routage** et **calcul d'itinéraires**.Il propose aussi des **méthodes spatiales** (carroyages, KDE, segments de lignes) et des **visualisations associées** pour restituer graphiquement la structure spatiale de semis de points ou de lignes. En **utilisant exclusivement des informations issues d'OSM**, l'exercice consiste à **identifier les secteurs accessibles à pied dans le voisinage géographique d'un lieu défini par une adresse**. ## Reproduire l'analysePour reproduire l'analyse, télécharger l'archive ZIP associée au dépôt [https://github.com/riatecom/osm-elementr-2025-application](https://github.com/riatecom/osm-elementr-2025-application), ouvrir le fichier proj.Rproj puis le fichier script.R qui reproduit en tout point le contenu du support de formation. N'oubliez pas d'installer les packages en amont via la commande inscrite dans le README.Cette analyse est reproductible en n'importe quel lieu, avec une adresse (adr) et un nom de ville (lib). Néanmoins les temps de calcul utilisant le serveur de démonstration d'OSRM ont été paramétrés pour ne pas le surcharger. L'exécution de ces blocs de code prennent du temps. En conséquence, 6 jeux de données ont été préparés en amont à partir de l'adresse des gares SNCF de Lille, Toulouse, Grenoble, Nantes, Clermont-Ferrand et Compiègne, ainsi que de la station de métro Front Populaire. Les blocs de code relevant des calculs d'itinéraires / temps de trajets sont placés en commentaires dans le fichier script.R, mais sont fonctionnels.```{r}cs <-read.csv("data/case_studies.csv", encoding ="UTF-8", sep =",")cs```Nous présentons ici les résultats depuis la **gare SNCF de Compiègne**. ```{r}cs <- cs[cs$lib =="Compiegne",]adr <- cs$adrlib <- cs$lib```## 0. Packages```{r, warning = FALSE, message = FALSE}# 1. Packages qui interfacent des données OSM#install.packages('maposm', repos = 'https://riatelab.r-universe.dev')library(tidygeocoder) # géocodagelibrary(maptiles) # import de tuiles OSM (raster)library(osmdata) # import de données OSM (vecteur)library(maposm) # import de données OSM (couches géo)library(osrm) # calcul d'itinéraires# 2. Autres packages utilitaireslibrary(sf) # manipulation de données vectorielleslibrary(terra) # manipulation de données rasterlibrary(mapsf) # cartographie thématiquelibrary(mapiso) # discrétiser des données continueslibrary(maplegend) # légendeslibrary(spatstat) # analyse de semis de pointslibrary(mapview) # cartographie interactivelibrary(stplanr) # segmentiser des lignes```## 1. GéocodageLe package `tidygeocoder`[@R-tidygeocoder] permet d'utiliser des services de géocodage en ligne. Par défaut, le package utilise [Nominatim](https://nominatim.org/) qui se base sur les données OSM.```{r}pt <-geo(address = adr) |>st_as_sf(coords =c("long", "lat"), crs =4326)pt```## 2. Import de tuiles OSMLe package `maptiles`[@R-maptiles] permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au **format raster.** [20 niveaux de zoom](https://wiki.openstreetmap.org/wiki/Zoom_levels) sont disponibles, selon la taille de l'objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.Plusieurs modèles de tuiles sont disponibles [ici](https://wiki.openstreetmap.org/wiki/Raster_tile_providers), certains nécessitent une inscription.```{r}emprise <-st_buffer(pt, 3000) |>st_bbox()tiles <-get_tiles(x = emprise,project =FALSE, # crs en output : 3857crop =TRUE,cachedir ="cache",zoom =14)pt_reproj <-st_transform(pt, crs ="EPSG:3857")mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt =paste0("Gare de ", lib))mf_title("Tuiles OSM - Défaut")# Autre exemple de fournisseurtiles <-get_tiles(x = emprise,project =FALSE,crop =TRUE,provider ="OpenStreetMap.HOT",zoom =14)mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt =paste0("Gare de ", lib))mf_title("Tuiles OSM - OpenStreetmap.HOT")```## 3. Import de couches géographiques OSMPour la cartographie, le package `maposm` (pas sur le CRAN) se base sur les fonctions du package `osmdata`[@R-osmdata] (partie suivante) pour télécharger les données OSM, au **format vectoriel**.La fonction principale, `om_get()`, permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc) autour d'un point, dans un rayon donné. Le résultat prend la forme d'une liste d'objets sf.```{r}#| fig-width: 7#| fig-height: 7res <-om_get(x =c(st_coordinates(pt)[1],st_coordinates(pt)[2]),r =2000)names(res)mf_map(res$zone, col ="#f2efe9", border =NA, add =FALSE)mf_map(res$urban, col ="#e0dfdf", border ="#e0dfdf", lwd = .5, add =TRUE)mf_map(res$green, col ="#c8facc", border ="#c8facc", lwd = .5, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$building, col ="#d9d0c9", border ="#c6bab1", lwd = .5, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Autour de la gare")mf_arrow()mf_scale(size =500, scale_units ="m")mf_credits("OpenStreetMap contributors, 2025", pos ="bottomleft")```## 4. Import de tags avec osmdataLe package `osmdata` permet d'extraire des données OSM à partir de l'API [Overpass](https://overpass-api.de/). Il nécessite dans un premier temps d'établir une *bounding box* au sein de laquelle seront extraites les données, puis de l'exécuter (fonction `add_osm_feature()`). Le résultat est une liste, qu'il faut ensuite convertir en objet sf.Les objets OSM prennent la forme de *nodes* (points), *ways* (lignes ou polygones), ou *relations* (relation entre plusieurs objets). Ces éléments sont décrits par un ensemble de *tags*, recensées sous forme de clé-valeurs (*key-values*). Ici, nous recherchons dans la clé **amenity** les points (*nodes*) correspondant aux restaurants, bars et cafés (*value* = restaurant, bar, cafe). En dehors des clés-valeurs principales, d'autres attributs peuvent être présents, comme en l'occurence le type de cuisine, ou bien les horaires.```{r, warning = FALSE}# Définition d'une bounding box (emprise Compiègne)q <- opq(bbox = emprise, osm_types = "node") # seulement les points# Extraction des restaurants, bars et cafésreq <- add_osm_feature(opq = q, key = 'amenity', value = c("bar", "cafe", "restaurant"))poi <- osmdata_sf(req)poi <- poi$osm_points# Sélectionner les variablespoi <- poi[, c("osm_id", "name", "amenity", "cuisine")]# Intersection et reprojectionpoi <- st_intersection(poi, st_transform(res$zone, crs = 4326))poi_reproj <- st_transform(poi , crs = "EPSG:3857")```La fonction `om_map()` affiche les couches géographiques récupérées.```{r}#| fig-width: 7#| fig-height: 7# Affichageom_map(res, title ="Autour de la gare")mf_map(poi_reproj, type ="typo", var ="amenity", cex = .7, pch =15, leg_pos ="topright", leg_val_cex = .9, leg_title_cex = .9, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, lwd =2, col ="darkblue", add =TRUE)```Il est possible d'utiliser le package `mapview`[@R-mapview] pour afficher les points ainsi que leurs attributs, de manière interactive. Dans le cas de données issues d'OSM, il est intéressant d'observer le nombre d'attributs disponibles par point.```{r}mapview(res$zone, alpha.regions =0, layer.name ="Zone") +mapview(poi, layer.name ="Manger et boire")```::: callout-note### osmextractLe package `osmextract`[@R-osmextract] permet d’extraire des données directement depuis une base de données OSM, sans avoir à passer par une API. Le plus petit niveau géographique auquel il est possible de télécharger les données correspond en France aux anciennes régions. Pour travailler sur Compiègne, il nous faudrait d'abord télécharger l'ensemble des données OSM de la Picardie, disponibles [ici](https://download.geofabrik.de/europe/france/picardie.html), et pesant plus de 240 Mb.Le code ci-dessous effectue les mêmes manipulations que précédemment.```{r, eval = FALSE}library(osmextract)osm_pt <- oe_get(place = "Picardie", layer = "points", extra_tags = "amenity", quiet = TRUE)osm_pt <- osm_pt[osm_pt$amenity %in% c("bar", "cafe", "restaurant"), ]osm_pt <- st_transform(osm_pt, crs = "EPSG:3857")osm_pt <- st_filter(osm_pt, res$zone, .predicates = st_intersects)```:::## 5. CarroyageUne fois les points d'intérêt récupérés (aussi appelés **POI**, pour *Points of Interest*), il est possible de les agréger dans une maille.```{r}#| fig-width: 7#| fig-height: 7# Création d'une grille régulière hexagonale de 300m de résolutiongrid <-st_make_grid(res$zone, cellsize =300, square =FALSE)grid <-st_intersection(grid, res$zone)grid <-st_sf(ID =1:length(grid), geom = grid)# Compter les POIs dans les carreauxinter <-st_intersects(grid, poi_reproj, sparse =TRUE)grid$n_poi <-lengths(inter)# Cartographieom_map(res)mf_map(grid, var ="n_poi", type ="choro", breaks =c(0,1,3,5,10,15),alpha = .6, border =NA, leg_title ="Nombre \nd'aménités",leg_val_rnd =0, leg_pos ="topright", add =TRUE)mf_title("Aménités gustatives - données carroyées")```## 6. Lissage (KDE)Le package `spatstat`[@R-spatstat] regroupe de nombreuses fonctions pour analyser des semis de points. Il génère des objets de classe **ppp** *(planar point pattern)*, dans des fenêtres d'observation, de classe **owin** *(observation window)*.La fonction `density.ppp()` calcule des valeurs à partir d'une estimation par noyau (KDE). Elle prend en entrée une portée (*sigma*), ainsi qu'une résolution au sein de laquelle sera calculée la densité, en pixels (*eps*).```{r}#| fig-cap: "@moraga2023"#| echo: falseknitr::include_graphics(rep('img/kde.png'))``````{r}p <-as.ppp(X =st_coordinates(poi_reproj), W =as.owin(res$zone))ds <-density.ppp(x = p, sigma =100, eps =10, positive =TRUE)```Une fois les densités calculées, nous créons un objet de type *SpatRaster*. La densité étant calculée sur un mètre carré, nous la multiplions par 10.000 afin de récupérer une densité par hectare.```{r, warning = FALSE}#| fig-width: 7#| fig-height: 7# Calcul densité de restaurants par hectarer <- rast(ds) * 100 * 100crs(r) <- st_crs(poi_reproj)$wktrpal <- c( "#418e40", "#a5be6a", "#e2cf6a", "#c19849", "#ab6c32", "#842f1a")mf_map(res$zone, col = "#f2efe9", border = NA, expandBB = c(0,.3,0,0))mf_raster(r, leg_title = "N. restaurants / ha.", leg_pos = "left", pal = pal, add = TRUE)mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(poi_reproj, pch = 20, cex = .3, col = "red", add = TRUE)mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)mf_title("La montagne culinaire - Compiègne")```Le package `mapiso`[@R-mapiso] permet de transformer des grilles régulières (ici notre raster) en contours de polygones (format vectoriel).```{r}#| fig-width: 7#| fig-height: 7# Discrétiser le rasterbks <- bks <-c(0, 0.1, 0.5, 1, minmax(r)[2])isos <-mapiso(x = r, var ="lyr.1 ", breaks = bks)mf_map(res$zone, col ="#f2efe9", border =NA)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5,add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(isos, "isomin", "choro", border =NA, alpha =0.5, breaks = bks,leg_title ="Densité\n (n/ha)", add =TRUE)mf_map(poi_reproj, pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Densités de restaurants - vectorisé")```## 7. Temps d'accèsLe package `osrm`[@R-osrm] interface l'engin de routage [OSRM](https://project-osrm.org/) et propose différentes fonctions pour le calcul de temps de trajet et d'itinéraires. Ces calculs s'effectuent sur la base de **profils .lua** qui pour chaque mode de transport (voiture, vélo, à pied) définissent des algorithmes permettant la sélection de telle ou telle route, une vitesse par défaut selon le type de route, etc. Le profil piéton est disponible [ici](https://github.com/Project-OSRM/osrm-backend/blob/master/profiles/foot.lua).D'après la [documentation de l'API OSRM](https://project-osrm.org/docs/v5.24.0/api/#), le nombre de requêtes **sur le serveur de démo** est limité à **512 par connection à l'API, espacées d'au moins 5 secondes**.Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d'installer une instance OSRM sur son propre serveur (voir [ici](https://rcarto.github.io/posts/build_osrm_server/)).### 7.1. Depuis la garePour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la gare de Compiègne et l'ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction `osrmTable()`, qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d'origine et de destination, et d'un vecteur correspondant à la métrique calculée. **Il est possible de calculer les distances et temps de trajet entre plusieurs couples origine / destination en une seule requête.**Nous construisons d'abord une matrice O/D composée de coordonnées X/Y des points d'origine et de destination (CRS 4326).::: callout-warning#### Utilisation de l'API OSRMPar souci de respect d'utilisation de l'API, une partie du code ci-dessous est passé en commentaire et non joué. Le résultat est enregistré sous format RDS pour être utilisé par la suite. Dans les morceaux de code suivants, le code apparaît dans des chunks sous option *eval = FALSE.*:::```{r}# Création de la matrice O/Do <-data.frame(X =st_coordinates(pt)[, 1],Y =st_coordinates(pt)[, 2])row.names(o) <-1:nrow(o)poi <- poi[poi$amenity =="restaurant",]d <-data.frame(X =st_coordinates(poi)[, 1],Y =st_coordinates(poi)[, 2])# 1 requête : 1 point d'origine, 50 points de destination# mat_dist <- osrmTable(src = o,# dst = d,# measure = "distance",# osrm.profile = 'foot')# saveRDS(mat_dist, "data/mat_dist.rds")# Intervalle de 5 secondes entre les requêtes# Sys.sleep(5)# mat_dur <- osrmTable(src = o,# dst = d,# measure = "duration",# osrm.profile = 'foot')# saveRDS(mat_dur, "data/mat_dur.rds")mat_dist <-readRDS("data/mat_dist.rds")mat_dur <-readRDS("data/mat_dur.rds")dist <-t(mat_dist$distances)dur <-t(mat_dur$durations)poi$dist <- distpoi$dur <- durapply(poi$dur, 2, summary)apply(poi$dist, 2, summary)```Il est possible d'atteindre `r round(nrow(d)/2, -1)` restaurants, cafés ou bars en moins de `r ceiling(median(poi$dur))` minutes à pied depuis la gare de Compiègne, et les trois quarts se situent à moins de `r round(apply(poi$dist, 2, summary)[5],-1)` mètres à pied.Les restaurants à moins de 5 minutes à pied de la gare sont les suivants :```{r}poi$name[poi$dur <5]```La fonction `osrmIsochrone()` du package `osrm` permet de calculer des isochrones depuis un point.```{r, eval = FALSE}isos <- osrmIsochrone(loc = pt, breaks = seq(0,45,5), res = 60, osrm.profile = "foot")st_write(isos, "data/mat.gpkg", layer = "isos", delete_layer = TRUE)``````{r, warning = FALSE}#| fig-width: 8#| fig-height: 7isos <- st_read("data/mat.gpkg", layer = "isos", quiet = TRUE)isos <- st_transform(isos, crs = "EPSG:3857")isos <- st_intersection(isos, res$zone)om_map(res, title = "Isochrones autour de la gare de Compiègne")mf_map(isos, type = "typo", var = "isomax", alpha = .6, pal = "Geyser", border = "white", leg_pos = NA, add = TRUE)mf_map(st_transform(poi, crs = "EPSG:3857"), var = "dur", type = "choro", breaks = sort(unique(c(isos$isomin, isos$isomax))), pal = "Geyser", pch = 21, border = "white", leg_pos = NA, add = TRUE)mf_map(pt_reproj, pch = 24, cex = 1.3, lwd = 2, col = "darkblue", add = TRUE)leg(type = "choro", val = sort(unique(c(isos$isomin, isos$isomax))), pal = "Geyser", title = "Temps\n(minutes)", pos = "left")```### 7.2. Depuis la grilleAfin de construire des indicateurs d'accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville **à pied depuis chaque point de grille.**```{r}nrow(grid) *nrow(d)```Nous avons en tout `r nrow(grid) * nrow(d)` couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos **189** points d'origine et **un** point de destination, sous forme de boucle.Le code pour effectuer les calculs est présenté ci-dessous, mais a été joué en local. Les résultats sont enregistrés dans le dossier *data*.```{r, warning = FALSE}# Coordonnées des points d'origine : centroïdes des cellules de la grilleo <- st_centroid(grid) |> st_transform(crs = 4326)o$X <- st_coordinates(o)[, 1]o$Y <- st_coordinates(o)[, 2]o <- st_set_geometry(o[, c("X", "Y")], NULL)```::: callout-important#### osrmNearest()Pour que l'engin de routage calcule un itinéraire, les points d'origine et de destination doivent être localisés sur une route. Si le point d'entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l'itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d'espace (urbain ou rural).:::```{r, eval = FALSE}matfoot <- list()for(i in 1:nrow(d)){ matfoot[[i]] <- osrmTable(src = o[, c("X", "Y")], dst = d[i , c("X", "Y")], measure = "duration", osrm.profile = 'foot') Sys.sleep(5) matfoot[[i]]$ID_ORIG <- c(1:nrow(o)) matfoot[[i]]$ID_DEST <- rep(i, nrow(o)) matfoot[[i]] <- data.frame("ID_ORIG" = matfoot[[i]]$ID_ORIG, "ID_DEST" = matfoot[[i]]$ID_DEST, "duration" = matfoot[[i]]$durations) colnames(matfoot[[i]]) <- c("ID_ORIG", "ID_DEST", "dur_foot")}matfoot <- do.call(rbind, matfoot)write.csv(matfoot, "data/matfoot.csv", row.names = FALSE)```La carte ci-dessous est un exemple d'indicateur d'accessibilité basé sur le temps de trajet : il s'agit du **temps à parcourir pour atteindre le service le plus proche**. Cet indicateur, que l'on peut décliner pour calculer le deuxième voire le troisième service le plus proche, et ainsi calculer la différence entre le premier et le second service, ouvre la question du **choix de la destination** : un individu est-il prêt à parcourir une plus grande distance pour bénéficier d'un meilleur / différent service ?```{r, warning = FALSE}#| fig-width: 7#| fig-height: 7 mat <- read.csv("data/matfoot.csv")# Restaurant le plus proche pour chaque point de grillepoi_min <- aggregate(dur_foot ~ ID_ORIG, mat, FUN = min, na.rm = TRUE)# Récupérer l'identifiant de destination associé# (Important pour la question du choix de la destination)poi_min <- do.call(rbind, lapply(split(mat, mat$ID_ORIG), function(x) x[which.min(x$dur_foot), ]))# Cartographierpoi_min_sf <- merge(grid, poi_min, by.x = "ID", by.y = "ID_ORIG", all.x = TRUE)om_map(res, title = "Temps d'accès au restaurant le plus proche")mf_map(poi_min_sf, var = "dur_foot", type = "choro", breaks = seq(0,15,3), alpha = .6, border = NA, leg_title = "Temps\n(minutes)", col_na = "#8B000060", leg_no_data = "> 15", leg_val_rnd = 0, pal = "Greens", leg_pos = "topright", add = TRUE)```::: callout-important#### ValhallaValhalla est un engin de routage open source basé sur les données OSM, au même titre que OSRM. Il fonctionne sur un système de **tuiles** contenant différents niveaux de routes, ainsi que les graphes routiers associés (les tuiles de niveau 0 contiennent les autoroutes, celles de niveau 2 les routes départementales). Les routes sont calculées au moment de la requête (tandis qu'elles sont pré-calculées sur OSRM) via différents algorithmes (A* pour les plus courts chemins, *Contraction Hierarchies* pour les longues distances). Valhalla propose de nombreux profils (vélo de location, multimodalité, taxi), ainsi que la possibilité de prendre en compte les pentes.:::Un autre exemple d'indicateur d'accessibilité : les **opportunités cumulées** : combien de restaurants sont disponibles sous un certain seuil de temps ? Parmi eux, combien de restaurants offrant une cuisine spécifique ?```{r}#| fig-width: 7#| fig-height: 7# moins de 5 minutesmat5 <- mat[mat$dur_foot <5 ,]mat5$n <-1mat5 <-aggregate(list(n_resto = mat5$n),list(ID_ORIG = mat5$ID_ORIG),FUN = sum)grid <-merge(grid[, "ID"], mat5,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)grid$n_resto[is.na(grid$n_resto)] <-0om_map(res, title ="Restaurants à moins de 5 minutes à pied")mf_map(grid, var ="n_resto", type ="choro", breaks =c(0,1,5,10,20,max(grid$n_resto)),alpha = .6, border =NA, leg_title ="Nombre de\nrestaurants",leg_val_rnd =0, pal ="Reds", leg_pos ="topright", add =TRUE)```## 8. ItinérairesPour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l'aide de la fonction `osrmRoute()`. Attention, l'API route d'OSRM permet de calculer des itinéraires seulement entre deux points. **Une requête correspond donc à un couple origine / destination.**```{r, eval = FALSE}routes <- list()o <- data.frame(X = st_coordinates(pt)[, 1], Y = st_coordinates(pt)[, 2])row.names(o) <- 1for(i in 1:nrow(d)){ routes[[i]] <- osrmRoute(src = o, dst = d[i,], overview = "full", osrm.profile = 'foot') Sys.sleep(5)}routes <- do.call(rbind, routes)st_write(routes, "data/mat.gpkg", layer = "routes", delete_layer = TRUE)```Le package `stplanr`[@R-stplanr] permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l'utilisons ici à des fins cartographiques, pour découper des lignes en segments.```{r, warning = FALSE}#| fig-width: 7#| fig-height: 7routes <- st_read("data/mat.gpkg", layer = "routes", quiet = TRUE)routes <- st_transform(routes, crs = "EPSG:3857")routes$n <- 1routes$distance <- routes$distance * 1000routes <- routes[order(routes$distance, decreasing = TRUE), ]l_seg <- list()# La fonction line_segment découpe un objet sf LINESTRING en segments réguliersfor(i in 1:nrow(routes)){ l_seg[[i]] <- line_segment( routes[i, ], segment_length = max(routes$distance) / 20)}pal <- mf_get_pal(nrow(l_seg[[1]]), rev = TRUE, palette = "Blues 3")om_map(res, title = "Routes", theme = "dark")for(i in nrow(routes):1){ mf_map(l_seg[[i]], col = pal, lwd = 1.6, add = TRUE)}leg(type = "cont", val = c(min(routes$distance), 1000, max(routes$distance)), pal = pal, pos = "topright", title = "Distance\n(mètres)")# La fonction overline permet de compter le nombre d'occurence des tronçonsseg <- overline(routes, attrib = "n", fun = sum)om_map(res, title = "Routes")mf_map(seg, var = "n", type = "grad", col= "darkred", breaks = "fisher", nbreaks = 5, lwd = c(1,3,5,9,13), leg_val_rnd = 0, leg_pos = "topright", leg_title = 'Fréquence', add = TRUE)```## Bonus : la tournée des pizzeriasJe fais une escale d'une demie-heure à Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps depuis la gare. Avec la fonction `osrmTrip()`, c'est possible !Cette fonction répond au [problème du voyageur de commerce](https://fr.wikipedia.org/wiki/Probl%C3%A8me_du_voyageur_de_commerce), qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.```{r, eval = FALSE}# Coordonnées des pizzeriapizz <- poi[!is.na(poi$cuisine) & poi$cuisine == "pizza" , "geometry"]# Ajouter le point de départ : la garepizz <- rbind(pt[1, "geometry"], pizz)pizz$id <- 1:nrow(pizz)trip <- osrmTrip(loc = pizz, overview = "full", osrm.profile = "foot")trip <- trip[[1]]$tripst_write(trip, "data/mat.gpkg", layer = "trip", delete_layer = TRUE)``````{r}#| fig-width: 7#| fig-height: 7trip <-st_read("data/mat.gpkg", layer ="trip", quiet =TRUE)trips <-st_transform(trip, crs ="EPSG:3857")# Pour les labels : récupérer le début de la lignestart <- lwgeom::st_startpoint(trips)start <-do.call(rbind, start)start <-data.frame(id =rownames(trips), X = start[,1], Y = start[,2])start <-st_as_sf(start, coords =c("X", "Y"), crs ="EPSG:3857")om_map(res, theme ="pizza", title ="La route de la pizza")mf_map(trips, col ="black", lwd =3, lty =3, add =TRUE)mf_map(start[start$id !=1, ], pch =20, cex =2, col ="black", add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_label(start, var ="id", cex =1.2, overlap =TRUE, pos =4, col ="black", halo =TRUE)sum(trip$duration)sum(trip$distance)```D'après les données OSM, moins de 25 minutes suffisent pour faire la tournée des pizzas à Compiègne, le tout en moins de 2 kilomètres. On est large !!## [](img/logo.png)
