La base de données utilisée dans ce TP concerne la géolocalisation des accidents de la route à Paris. Il s’agit plus précisément des bases de données annuelles des accidents corporels de la circulation routière, en particulier le millésime 2019.
« Pour chaque accident corporel (soit un accident survenu sur une voie ouverte à la circulation publique, impliquant au moins un véhicule et ayant fait au moins une victime ayant nécessité des soins), des saisies d’information décrivant l’accident sont effectuées par l’unité des forces de l’ordre (police, gendarmerie, etc.) qui est intervenue sur le lieu de l’accident. Ces saisies sont rassemblées dans une fiche intitulée bulletin d’analyse des accidents corporels. L’ensemble de ces fiches constitue le fichier national des accidents corporels de la circulation dit « Fichier BAAC » administré par l’Observatoire national interministériel de la sécurité routière “ONISR”.
Un certain nombre d’indicateurs issus de cette base font l’objet d’une labellisation par l’autorité de la statistique publique (arrêté du 27 novembre 2019). »
Vous pouvez télécharger les bases de données brutes ici ou utiliser R pour les télécharger dans votre dossier actuel :
# télécharger le dataset
download.file("https://github.com/comeetie/quantilille/raw/main/exercises/data.zip",
destfile = "data.zip")
# dézipper
unzip("data.zip",exdir=".")
sf
sf::st_read()
.
library(sf)
.75 <- st_read("data/iris_75.shp", quiet = TRUE,
irisstringsAsFactors = F)
plot(iris.75)
. Que remarquez-vous ?
plot(iris.75)
sf
.
sf::st_geometry()
? Quelle solution au problème précédent proposez-vous ?
sf::st_geometry()
permet d’isoler l’information contenue dans la colonne geometry
de l’objet sf
. Cela permet de mettre de côté les autres variables et de n’en afficher qu’une.
plot(st_geometry(iris.75))
plot
.
sf::st_read()
et sf::st_geometry()
. Vous pouvez aussi customiser la carte en utilisant différents paramètres de la fonction plot
: bg, col, lwd, border, pch, cex…
2019.paris <- st_read("data/accidents2019_paris.geojson",
accidents.quiet = TRUE)
plot(st_geometry(iris.75), bg = "cornsilk", col = "lightblue",
border = "white", lwd = .5)
plot(st_geometry(accidents.2019.paris), col = "red", pch = 20, cex = .2, add=TRUE)
title("Accidents à Paris")
Effectuez ces quelques traitement géomatiques :
st_union
, st_buffer
, st_voronoi
, st_collection_extract
, st_interection
, st_join
et st_centroid
.
# 1. Carte agrégée Paris
75.u <- st_union(iris.75)
iris.# 2. Zone tampon de 1km
75.b <- st_buffer(x = iris.75.u, dist = 1000)
iris.# 3. Centroïdes
75.c <- st_centroid(iris.75)
iris.# 4. Distance entre centroïdes
<- st_distance(x=iris.75.c, y=iris.75.c)
mat # 5. Calculer les polygones de Voronoïs autour des centroïdes
75.v <- st_collection_extract(st_voronoi(x = st_union(iris.75.c)))
iris.75.v <- st_intersection(iris.75.v,iris.75.u)
iris.
# 6. Affichage des cartes
plot(st_geometry(iris.75.b), lwd=2, border ="red",col=NA)
plot(st_geometry(iris.75), ltw=5, col="#999999", add = TRUE)
plot(st_geometry(iris.75.u), border="blue", ltw=5, col=NA, add = TRUE)
plot(st_geometry(iris.75.c), pch = 20, cex = .2,col="red", add = TRUE)
plot(st_geometry(iris.75.v), ltw=5, col=NA,border="blue", add = TRUE)
Cet exercice vient dans le prolongement de l’exercice 1.
Méthode 1 : Utilisez sf::st_intersects()
et sapply()
.
Méthode 2 (approche tidyverse) : Utilisez sf::st_join()
(jointure spatiale), dplyr::count
et dplyr::left_join
.
Documentation de la variable grav
: Gravité de blessure de l’usager, les usagers accidentés sont classés en trois catégories de victimes plus les indemnes :
library(dplyr)
# # Méthode 1 : st_intersect
<- st_intersects(x = iris.75, y = accidents.2019.paris)
inter <- st_intersects(x = iris.75, y = accidents.2019.paris %>% filter(grav==1))
inter_nonbless .75$nbacc <- sapply(inter, length)
iris.75$nbaccnb <- sapply(X = inter_nonbless, FUN = length)
irishead(iris.75)
Simple feature collection with 6 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 648979.6 ymin: 6861839 xmax: 654353.1 ymax: 6866159
Projected CRS: RGF93 / Lambert-93
CODE_IRIS INSEE_COM geometry nbacc nbaccnb
1 751197316 75119 MULTIPOLYGON (((653970.6 68... 4 2
2 751176716 75117 MULTIPOLYGON (((649189.3 68... 7 3
3 751103703 75110 MULTIPOLYGON (((652767.6 68... 14 5
4 751187104 75118 MULTIPOLYGON (((652827.6 68... 4 2
5 751114314 75111 MULTIPOLYGON (((654272.9 68... 11 4
6 751103707 75110 MULTIPOLYGON (((652960.8 68... 11 5
# # Méthode 2 : jointure et agrégation spatiale
# accidents.2019.paris.iris <- iris.75 %>% st_join(st_as_sf(accidents.2019.paris))
# iris.75 <- iris.75 %>%
# # sans le st_drop_geometry, le traitement est beaucoup plus long
# left_join(accidents.2019.paris.iris %>% st_drop_geometry() %>%
# count(CODE_IRIS) %>% setNames(c("CODE_IRIS","nbacc"))) %>%
# left_join(accidents.2019.paris.iris %>% st_drop_geometry() %>% filter(grav==1) %>%
# count(CODE_IRIS) %>% setNames(c("CODE_IRIS","nbaccnb"))) %>%
# mutate(nbaccnb = ifelse(is.na(nbaccnb),0,nbaccnb))
dplyr
: select
, group_by
et summarize
. Ces fonctions fonctionnent également avec les objets sf
.
library(dplyr)
.75 <- iris.75 %>%
comgroup_by(INSEE_COM) %>%
summarize(nbacc = sum(nbacc),
nbaccnb = sum(nbaccnb))
plot(st_geometry(iris.75), col = "ivory3", border = "ivory1")
plot(st_geometry(com.75), col = NA, border = "ivory4", lwd = 2, add = TRUE)
voie
(et éventuellement com
) les personnes ayant eu un accident dans le 2ème arrondissement de Paris et en étant sorties indemnes.
Utilisez les deux méthodes de géocodage :
library(banR)
pour interroger l’API de la BAN françaiselibrary(tidygeocoder)
pour les adresses internationales
# Jointure spatiale
2019.paris.iris <- iris.75 %>% st_join(st_as_sf(accidents.2019.paris))
accidents.
# Avec banR
# remotes::install_github("joelgombin/banR")
library(banR)
<- accidents.2019.paris.iris %>%
geo_banR filter(grav == 1 & INSEE_COM == "75102") %>%
geocode_tbl(adresse = voie,code_insee = com) %>%
select(latitude,longitude) %>%
st_as_sf(coords = c("longitude", "latitude"),
crs = 4326, agr = "constant") %>%
st_transform(2154)
# Avec tidygeocoder
# install.packages("tidygeocoder")
library(tidygeocoder)
<- accidents.2019.paris.iris %>%
geo_tidygeocoder filter(grav == 1 & INSEE_COM == "75102") %>%
mutate(addr = paste(voie, ", Paris, France")) %>%
geocode(addr,method = "osm") %>% select(lat, long) %>%
st_as_sf(coords = c("long", "lat"),crs = 4326, agr = "constant") %>%
st_transform(2154)
Dans cet exercice, nous allons utiliser mapview
pour explorer les accidents de la route ayant eu lieu à Paris en 2019.
En complément, nous allons utiliser les données d’OSM permettant de croiser le lieu des accidents avec les routes empruntées.
mapview
. Essayez d’utiliser différents paramètres pour customiser votre carte. Coloriez notamment les points selon la gravité des accidents.
map.types
, col.regions
, label
, color
, legend
, layer.name
, homebutton
, lwd
, alpha
, zcol
… du package mapview
.
library(mapview)
library(sf)
2019.paris <- st_read("data/accidents2019_paris.geojson",
accidents.quiet = TRUE)
mapview(accidents.2019.paris %>%
mutate(grav_f = factor(grav,
levels = c(2,3,4,1),
labels = c("Tué", "Blessé hospitalisé",
"Blessé léger","Indemne")
)),col.regions = c("darkred","red","orange","darkgreen"),
label = accidents.2019.paris$Num_Acc, #label
color = "white", legend = TRUE,
zcol="grav_f", alpha=0.9,
map.types = "Stamen.TonerLite",
layer.name = "Gravite",
homebutton = FALSE, lwd = 0.2)