Pratiquons les données spatiales avec R !

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=".") 

Exercice 1 : Manipuler des objets sf

1

Importer la carte des iris¹ ‘iris.75.shp’ de Paris.

Indices

Utilisez la fonction sf::st_read().

Solution

library(sf)
iris.75 <- st_read("data/iris_75.shp", quiet = TRUE,
                   stringsAsFactors = F)

2

Afficher la carte de Paris grâce à l’instruction plot(iris.75). Que remarquez-vous ?

Solution

plot(iris.75)

Solution

On remarque que R fait plusieurs graphiques : un par variable contenue dans l’objet sf.

3

A quoi sert la fonction sf::st_geometry() ? Quelle solution au problème précédent proposez-vous ?

Solution

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.

Solution

plot(st_geometry(iris.75))

4

Importez la couche des accidents de la route appelée ‘accidents2019_paris.geojson’ et affichez la carte des accidents dans Paris en utilisant simplement la fonction plot.

Indices

Utilisez 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…

Solution

accidents.2019.paris <- st_read("data/accidents2019_paris.geojson",
                                quiet = TRUE)

Solution

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")

5

Effectuez ces quelques traitement géomatiques :

1. Créez un polygone avec le contour de Paris en agrégeant les iris de Paris (union simple)
2. Faire une zone tampon de 1km autour de celui-ci
3. Extraire les centroïdes des iris
4. Calculez les distances entre les centroïdes des iris
5. Calculez les polygones de Voronoï autour des centroïdes des iris
6. Affichez les différentes couches créées

Indices

Utilisez les fonctions st_union, st_buffer, st_voronoi, st_collection_extract, st_interection, st_join et st_centroid.

Solution

# 1. Carte agrégée Paris
iris.75.u <- st_union(iris.75)
# 2. Zone tampon de 1km
iris.75.b <- st_buffer(x = iris.75.u, dist = 1000)
# 3. Centroïdes
iris.75.c <- st_centroid(iris.75)
# 4. Distance entre centroïdes
mat <- st_distance(x=iris.75.c, y=iris.75.c)
# 5. Calculer les polygones de Voronoïs autour des centroïdes
iris.75.v <- st_collection_extract(st_voronoi(x = st_union(iris.75.c)))
iris.75.v <- st_intersection(iris.75.v,iris.75.u)

# 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)

Exercice 2 : Intersections, Agrégations

Cet exercice vient dans le prolongement de l’exercice 1.

1

Comptez le nombre de personnes accidentées par iris ainsi que le nombre de personnes accidentées mais non blessées. Vous pouvez utiliser les deux méthodes vues dans la partie cours.

Indices

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 :

• 1 : Indemne
• 2 : Tué
• 3 : Blessé hospitalisé
• 4 : Blessé léger

Solution

library(dplyr)
# # Méthode 1 : st_intersect
 inter <- st_intersects(x = iris.75, y = accidents.2019.paris)
inter_nonbless <- st_intersects(x = iris.75, y = accidents.2019.paris %>% filter(grav==1))
iris.75$nbacc <- sapply(inter, length)
iris.75$nbaccnb <- sapply(X = inter_nonbless, FUN = length)
head(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))

2

Utilisez la couche ‘iris.75’, pour créer une nouvelle couche cartographique agrégée appelée ‘com.75’ qui correspond aux ‘arrondissements’ de Paris. Gardez aussi dans cette nouvelle couche l’information sur le nombre de personnes accidentées et le nombre de personnes accidentées non blessées dans chaque arrondissement.

Information

La couche cartographique ‘iris.75’ contient un code de 5 chiffres dans sa variable INSEE_COM qui correspond au code de l’arrondissement.

Indices

Utilisez les fonctions du package classique dplyr : select, group_by et summarize. Ces fonctions fonctionnent également avec les objets sf.

Solution

library(dplyr)
com.75 <- iris.75 %>%
  group_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)

3

Géocodez à partir des variables voie (et éventuellement com) les personnes ayant eu un accident dans le 2ème arrondissement de Paris et en étant sorties indemnes.

Indices

Utilisez les deux méthodes de géocodage :

1. library(banR) pour interroger l’API de la BAN française
2. library(tidygeocoder) pour les adresses internationales

Solution

# Jointure spatiale
accidents.2019.paris.iris <- iris.75 %>% st_join(st_as_sf(accidents.2019.paris))

# Avec banR
# remotes::install_github("joelgombin/banR")
library(banR)
geo_banR <- accidents.2019.paris.iris %>% 
  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)
geo_tidygeocoder <- accidents.2019.paris.iris %>% 
  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)

Exercice 3 : Cartes interactives

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.

1

Chargez la base de données ‘accidents.2019.paris’ et affichez les positions des 11 897 personnes victimes d’un accident à Paris en 2019 grâce au package mapview. Essayez d’utiliser différents paramètres pour customiser votre carte. Coloriez notamment les points selon la gravité des accidents.

Information

Par exemple, vous pouvez utiliser les paramètres map.types, col.regions, label, color, legend, layer.name, homebutton, lwd, alpha, zcol … du package mapview.

Solution

library(mapview)
library(sf)

accidents.2019.paris <- st_read("data/accidents2019_paris.geojson", 
                                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)

2

1. Utilisez les polygones de ‘iris.75’ pour extraire la “bounding box” de Paris en projection WGS84.

2. Récupérez ensuite grâce à osmdata, à l’intérieur de cette bounding box, le fond de carte du périphérique parisien (key = "highway", value = "trunk")

3. Faites l’intersection entre les accidents et le périphérique, en prenant soin d’ajouter une zone tampon de 50 mètres autour de ce dernier et appeler ce nouvel ensemble de points accidents.2019.paris.periph.

Indices

Utiliser sf::st_bbox() et sf::st_transform() pour extraire la bounding box. Le code epsg de WGS84 est 4326.

Utiliser :

• osmdata:opq() pour définir la bounding box de la requête osm
• osmdata:add_osm_feature() pour définir la paire key:value recherchée
• osmdata:osmdata_sf() pour récupérer les données osm.

Solution

library(osmdata)

#1. bounding box
bb      <- iris.75 %>% st_transform(4326) %>% st_bbox()
q       <- opq(bbox = bb,timeout = 180)

#2. périphérique
qt      <- add_osm_feature (q, key = 'highway',value = 'trunk', value_exact = FALSE)
roads    <- c(osmdata_sf(qt))$osm_lines %>% st_transform(st_crs(iris.75))
Encoding(roads$name) <- "UTF-8" #si problème d'encodage

#3. zone tampon et intersection

# Méthode 1
accidents.2019.paris.periph <- st_intersection(
  accidents.2019.paris,
  st_intersection(st_geometry(roads),iris.75)  %>% st_buffer(dist = 50) %>% st_union()
)  

# Méthode 2
periph <- st_geometry(roads %>% filter(name=="Boulevard Périphérique Extérieur"))

periph_count <- st_sf(periph %>% st_buffer(50),
                      id=1:length(periph)) %>%
  st_join(accidents.2019.paris) %>% count(id)

3

Affichez les positions des 2 073 personnes victimes d’un accident à Paris SUR LE PERIPHERIQUE en 2019 grâce au package mapview. Essayez à nouveau d’utiliser différents paramètres pour customiser votre carte.

Information

Par exemple, vous pouvez utiliser les paramètres map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … du package mapview.

Solution

library(mapview)
library(sf)

mapview(accidents.2019.paris.periph, map.types = "OpenStreetMap",
        col.regions = "#940000",
        color = "white", legend = TRUE, layer.name = "Accidents",
        homebutton = FALSE, lwd = 0.2)

Bonus : compter les points dans une grille

4

Utiliser la fonction pt_in_grid ci-dessous pour compter le nombre de personnes accidentées dans des cellules de 500m de côté.

Utilisez mapview pour afficher la grille choroplèthe.

pt_in_grid <- function(feat, adm, cellsize = 1000){
  grid <- st_make_grid(x = adm, cellsize = cellsize, what = "polygons")
  . <- st_intersects(grid, adm)
  grid <- grid[sapply(X = ., FUN = length)>0]
  . <- st_intersects(grid, feat)
  grid <- st_sf(n = sapply(X = ., FUN = length), grid)
  return(grid)
}

Solution

library(RColorBrewer)
gr <- pt_in_grid(accidents.2019.paris,iris.75,500)
bks = quantile(gr$n)
cols <- brewer.pal(length(bks), "Reds")

mapview(st_as_sf(gr)  %>% st_transform(4326),
        map.types = "Stamen.TonerLite",
        color = "white",
        col.regions = cols,
        alpha = 0.9,
        at = bks, 
        legend = TRUE,
        layer.name = "Nombre d'accidents par</br>carreau de 500 metres",
        homebutton = FALSE, lwd = 0.2)

Exercice 4 : ggplot2 (graphique et carte)

Nous aimerions créer avec le package ggplot2 une carte des arrondissements de Paris qui combine le nombre de personnes accidentées et la part de celles qui n’ont pas été blessées.

Avant de passer aux cartes, nous allons nous entraîner à faire des graphiques avec ce même package ggplot2.

Graphique

1

En utilisant le package ggplot2, essayez de reproduire le graphique suivant.

Indices

Utilisez la fonction ggplot2::geom_barplot() pour faire un diagramme en barres. Et mettez en place la grammaire des graphiques dans ce cas précis.

1. Ne garder que les véhicules catv pour lesquels il y plus de 100 personnes accidentées ;
2. Créez une base de données avec une colonne prop qui correspond à la répartition en pourcentages de la gravité grav de l’accident (Tué, Blessé hospitalisé, Blessé léger et Indemne) pour chaque type de véhicule catv impliqué sélectionné dans le 1.
3. Veillez à transformer les deux variables en factor pour bien les ordonner

• grav dans l’ordre suivant : c("Tué", "Blessé hospitalisé","Blessé léger","Indemne")
• catv dans l’ordre du graphique qui correspond à l’ordre croissant de la part de personnes indemnes par catégorie de véhicule

4. Réalisez le graphique (geom_barplot, scale_fill_manual, theme_bw, geom_text, labs)

Solution

# 1. Ne garder que les véhicules +100 personnes accidentées
vehic_plus100_accidents <- accidents.2019.paris %>% st_drop_geometry %>% 
  count(catv) %>% filter(n>100) %>% pull(catv)

# 2. Créer bases de données avec pourcentages
gg <- accidents.2019.paris %>%
  st_drop_geometry() %>% 
  filter(catv %in% vehic_plus100_accidents) %>% 
  mutate(catv_f = factor(catv)) %>% 
  mutate(grav_f = factor(grav, #3. facteur bien ordonné grav_f
                         levels = c(2,3,4,1),
                         labels = c("Tué", "Blessé hospitalisé", 
                                    "Blessé léger","Indemne")
                         )) %>% 
  count(catv_f,grav_f) %>%  
  group_by(catv_f) %>% 
  mutate(prop = 100*n/sum(n))  

#  facteur bien ordonné catv_f
ordre_catv_f <- gg %>% 
  filter(grav_f=="Indemne") %>% 
  arrange(prop) %>% 
  select(catv_f) %>% 
  pull()

# 4. Graphique
library(ggplot2)
ggplot(data = gg)+
  geom_bar(aes(x = prop, y= factor(catv_f,levels=ordre_catv_f),fill = grav_f), stat="identity") + 
  scale_fill_manual("Gravité (%)",
                    values =  c("darkred","red","orange","darkgreen"))+
  theme_bw()+
  labs(title = "Répartition de la gravité des accidents par type de véhicule",
       subtitle = "à Paris en 2019",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",
       x = "", y = "Type de véhicule")

Carte

2

Préparation des données pour la carte :

• Chargez le fond de carte ‘com75_shp’ (qui contient le nombre de personnes accidentées, en tout et non blessées, dans chaque arrondissement) et créez une variable appelée part_non_blesses qui correspond à la part des personnes non blessées parmi les accidentés dans chaque arrondissement.
  
• Créez un vecteur des quartiles de la variable part_non_blesses.
  
• Créez le vecteur de couleur qui correspond au nombre de classes définies plus tôt.
• Ajouter une variable appelée typo à ‘com.75’ qui indique la classe de l’arrondissement selon la discrétisation contenue dans bks pour la variable part_non_blesses.
  

Information

Pour la création de ‘bks’ et de ‘cols’, utilisez les fonctions quantile et RColorBrewer::brewer.pal. Pour la création de la variable typo, vous pouvez utiliser la fonction cut avec les paramètres digit.lab = 2 et include.lowest = TRUE.

Solution

library(sf)
# Importer les données
com.75 <- st_read("data/com_75.shp", quiet = TRUE)
# Créer la variable
com.75$part_non_blesses <- 100 * com.75$nbaccnb / com.75$nbacc
# Définir les bks par quantile
bks <- quantile(com.75$part_non_blesses, na.rm = TRUE)
# Définir une palette de couleurs
library(RColorBrewer)
cols <- brewer.pal(length(bks)-1,"Greens")

# For ggplot2 maps - Create a "typo" variable
library(dplyr)
com.75 <- com.75 %>%
  mutate(typo = cut(part_non_blesses,
                    breaks = bks,
                    labels = paste0(
                      round(bks[1:(length(bks)-1)]),
                      " à ",round(bks[2:length(bks)])
                      ),
                    include.lowest = TRUE))

3

En utilisant le package ggplot2, créez une carte qui contient en choroplèthe la variable part_non_blesses et en cercles proportionnels la variable nbacc.

Solution

library(ggplot2)

map_ggplot <- ggplot() +
  geom_sf(data = com.75, aes(fill = typo), colour = "grey80") +
  scale_fill_manual(name = "Part des non-blessés parmi les\naccidentés de la route (en %)",
                    values = cols) +
  geom_sf(data = com.75 %>%  st_centroid(),
          aes(size = nbacc), fill = "#f5f5f5", color = "grey20", shape = 21, 
          stroke = 1, alpha = 0.8, show.legend = "point") +
  scale_size_area(max_size = 12, name = "Nombre de personnes\n accidentées") +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(com.75)[c(1,3)],
           ylim = st_bbox(com.75)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "cornsilk", color = NA), 
        legend.position = "bottom", plot.background = element_rect(fill = "cornsilk",color=NA)) +
  labs(title = "Accidents de la route à Paris en 2019",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021") +
  guides(size = guide_legend(label.position = "bottom", title.position = "top",
                             override.aes = list(alpha = 1, color = "#ffffff")),
         fill = guide_legend(label.position = "bottom", title.position = "top"))

plot(map_ggplot)

Exercice 5 : Données Airbnb

Ce dernier exercice s’appuie sur les données des listings Airbnb à Paris.

Elle sont disponibles sur le site Inside Airbnb.

1

Préparation des données Airbnb :

• Chargez le fichier listings.csv et convertissez le en data.frame sf.
• Explorez un peu ces données avec mapview
• Projetez le en Lambert 93 (code epsg 2154)

Solution

library(readr)
library(dplyr)
library(sf)
listings <- read_csv("data/listings.csv") %>% 
  st_as_sf(coords = c("longitude", "latitude"),
           crs = 4326,
           agr = "constant")

listings.l93 <- st_transform(listings,2154)
plot(listings.l93 %>% st_geometry(), col="#88888822")

2

Préparation des données OSM :

• Récupérez les localisations des stations de métros de Paris via le package osm data (key = 'station', value = 'subway').
• Faites en un data.frame spatiale propre
  • filtrez-la pour ne conserver que les lignes dont le champ name est renseigné
  • supprimez les doublons éventuels (sur le champ name).
• Projetez-la en Lambert 93 (code epsg 2154). Appelez cet objet subway.l93.

Solution

library(osmdata)
bb      <- listings %>% st_bbox()
q       <- opq(bbox = bb,timeout = 180)
qm      <- add_osm_feature (q, key = 'station',
                            value = 'subway',
                            value_exact = FALSE)
subway  <- osmdata_sf(qm)
subway.l93  <- st_transform(subway$osm_points,2154) %>%
  filter(!is.na(name),!duplicated(name))
Encoding(subway.l93$name) <- "UTF-8" #pour les problèmes d'encodage
plot(subway.l93 %>% st_geometry(), pch = 20, col = "red",  cex = 1)

3

Préparation des polygones de Voronoï :

1. Créez les polygones de Voronoï associés aux stations de métro ;
2. Faites l’intersection entre ces voronoïs et la bounding box des listings. Cet objet s’appellera subway.voronoi.

Solution

# 1. Création des polygones de Voronoïs
subway.voronoi <- st_collection_extract(
  st_voronoi(x = st_union(subway.l93 %>% st_geometry()))
  )
# 2. Bounding box et intersection
bb <- st_as_sfc(st_bbox(listings.l93),crs=2154)
subway.voronoi <- subway.voronoi %>% st_intersection(bb)    

# Carte pour visualiser                                 
plot(subway.voronoi) # Contour des voronoi
plot(listings.l93 %>% st_geometry(),
     pch = 20, cex = 1,add=TRUE,col="#55555511") #Listings airbnb
plot(subway.l93 %>% st_geometry(),
     pch = 20, col = "red", add=TRUE, cex = 1) #Stations de métro

4

Résumé statistique :

1. Créez l’objet subway.poly, un data.frame sf créé à partir des polygones subway.voronoi, qui continent l’identifiant osm des stations de métro, ainsi que leurs noms (contenus dans subway.l93).
2. Agrégez les listings sur les voronoïs et calculez le nombre d’annonces, le prix moyen et médian, ainsi que la variance des prix.

Solution

# 1. Création de subway.poly
subway.poly <- subway.l93 %>% select(osm_id,name)
subway.poly$geometry <- subway.voronoi

# 2. Agrégation des listings
## Au niveau des polygones
listings.voronoi <- subway.poly %>% st_join(listings.l93) %>% 
  group_by(osm_id,name.x) %>% 
  summarize(n=n(),
            price_med=median(price,na.rm=TRUE),
            price_mean=mean(price,na.rm=TRUE),
            price_sd=sd(price,na.rm=TRUE))
head(listings.voronoi)

## Au niveau des centroïdes
subway.point <- subway.l93 %>%
  select(osm_id,name) %>%
  left_join(listings.voronoi %>% st_drop_geometry())

Simple feature collection with 6 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 653378 ymin: 6857596 xmax: 654637.2 ymax: 6865268
Projected CRS: RGF93 / Lambert-93
# A tibble: 6 x 7
# Groups:   osm_id [6]
  osm_id  name.x        n price_med price_mean price_sd                           geometry
  <chr>   <chr>     <int>     <dbl>      <dbl>    <dbl>                      <POLYGON [m]>
1 122478~ Avron       677        80      105.      90.6 ((653868.7 6861305, 654029.4 6861~
2 123268~ Boulevar~   157        80       94.6     67.1 ((653709.6 6862344, 653795.5 6862~
3 123268~ Boulevar~   321        85       99.2     67.2 ((653609.2 6861866, 653813.9 6862~
4 123268~ Boulevar~    65        60       78.2     66.4 ((654005.2 6858581, 654500.5 6858~
5 123268~ Rue Chau~   167        65      102.     383.  ((653635.9 6864921, 653906.9 6865~
6 123268~ Rue Tait~   468        70       82.1     55.7 ((653924.7 6863926, 654271.7 6864~

5

Essayez de reproduire cette figure :

Information

Vous pouvez vous servir de ce vecteur de sélection de quelques stations de métro.

sel_metro <- c('Filles du Calvaire', 'Colonel Fabien','Saint-Placide','Rennes', 'Les Halles','Tolbiac', 'Denfert-Rochereau','Oberkampf','Montparnasse-Bienvenüe','La Tour Maubourg','Mairie des Lilas', 'École Militaire','Saint-Germain-des-Prés','Boulogne - Pont de Saint-Cloud',"Gare de l'Est", "Place d'Italie", 'Richard Lenoir','Saint-Lazare','Porte de la Villette','Palais Royal - Musée du Louvre','Château Rouge', 'Gare du Nord (Métro)','Strasbourg - Saint-Denis', 'Poissonnière','Balard','Gare de Lyon',"Porte d'Italie", 'Gare d'Austerlitz')

Indices

Utilisez la fonction ggplot2::geom_boxplot() pour faire une boîte à moustache. Et mettez en place la grammaire des graphiques dans ce cas précis.

Solution

sel_metro <- c('Filles du Calvaire',
               'Colonel Fabien','Saint-Placide','Rennes',
               'Les Halles','Tolbiac','Denfert-Rochereau',
               'Oberkampf','Montparnasse-Bienvenüe',
               'La Tour Maubourg','Mairie des Lilas',
               'École Militaire','Saint-Germain-des-Prés',
               'Boulogne - Pont de Saint-Cloud',"Gare de l'Est",
               "Place d'Italie",'Richard Lenoir','Saint-Lazare',
               'Porte de la Villette','Palais Royal - Musée du Louvre',
               'Château Rouge','Gare du Nord (Métro)',
               'Strasbourg - Saint-Denis','Poissonnière',
               'Balard','Gare de Lyon',"Porte d'Italie",
               "Gare d'Austerlitz")


gg <- subway.poly %>% st_join(listings.l93) %>% 
  st_drop_geometry() %>% 
  filter(name.x %in% sel_metro, price < 750) %>%  
  mutate(name.x=factor(name.x,levels=sel_metro))

library(ggplot2)
ggplot(gg) +
  geom_boxplot(aes(y = name.x, x = price)) + 
  theme_bw()+ 
  labs(title = "Répartition des prix d'une nuitée à Paris",
       subtitle = "à proximité de certaines stations de métro",
       caption = "fichier inside Airbnb,\nantuki & comeetie, 2021",
       x = "Prix d'une nuitée (€)",y = "")

6

Essayez de reproduire cette carte :

Indices

Les étiquettes s’affichent avec geom_sf_text et correspondent aux 10 lieux avec le plus d’hôtes.

Solution

bks <- round(quantile(subway.point$price_med,
                      na.rm=TRUE,
                      probs=seq(0,1,0.2)))
subway.point <- subway.point %>%
  mutate(price_cat= cut(price_med,bks))

roads.geom <- st_read(dsn = "data/road.shp", quiet = TRUE)
river.geom <- st_read(dsn = "data/river.shp", quiet = TRUE)

library(ggplot2)
library(ggspatial)
ggplot(subway.point)+
  geom_sf(aes(size=n,color=price_cat)) + 
  geom_sf(data = river.geom, colour = "#87cdde",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  geom_sf_text(data=subway.point%>% arrange(desc(n)) %>% head(),aes(label=name))+
  scale_color_brewer("Prix médian",palette="Greens") + 
  scale_size(name = "Nombre d'hôtes",
             breaks = c(1,100,250,500,750),
             range = c(0,5)) +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(listings.l93)[c(1,3)],
           ylim = st_bbox(listings.l93)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "ivory",color=NA),
        plot.background = element_rect(fill = "ivory",color=NA)) +
  labs(title = "Prix médian d'une nuitée Airbnb",
       subtitle = "suivant la station de métro la plus proche.",
       caption = "fichier inside Airbnb, \nantuki & comeetie, 2021",
       x = "",y = "")

---

---

Reproducibilité

sessionInfo()

R version 4.0.5 (2021-03-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)

Matrix products: default

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
[5] LC_TIME=French_France.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggspatial_1.1.5    readr_1.4.0        ggplot2_3.3.3      RColorBrewer_1.1-2
 [5] osmdata_0.1.5      mapview_2.9.9      tidygeocoder_1.0.3 banR_0.2.2        
 [9] dplyr_1.0.5        sf_0.9-8           knitr_1.33        

loaded via a namespace (and not attached):
 [1] httr_1.4.2              sass_0.3.1              jsonlite_1.7.2         
 [4] bslib_0.2.4             assertthat_0.2.1        sp_1.4-5               
 [7] highr_0.9               stats4_4.0.5            yaml_2.2.1             
[10] pillar_1.6.0            lattice_0.20-41         glue_1.4.2             
[13] uuid_0.1-4              digest_0.6.27           rvest_1.0.0            
[16] colorspace_2.0-0        leaflet.providers_1.9.0 htmltools_0.5.1.1      
[19] pkgconfig_2.0.3         raster_3.4-10           purrr_0.3.4            
[22] scales_1.1.1            webshot_0.5.2           brew_1.0-6             
[25] svglite_2.0.0           satellite_1.0.2         tibble_3.1.1           
[28] proxy_0.4-25            generics_0.1.0          farver_2.1.0           
[31] ellipsis_0.3.2          withr_2.4.2             cli_2.5.0              
[34] magrittr_2.0.1          crayon_1.4.1            mime_0.10              
[37] evaluate_0.14           fansi_0.4.2             xml2_1.3.2             
[40] class_7.3-18            tools_4.0.5             hms_1.0.0              
[43] lifecycle_1.0.0         stringr_1.4.0           munsell_0.5.0          
[46] compiler_4.0.5          jquerylib_0.1.4         e1071_1.7-6            
[49] systemfonts_1.0.1       rlang_0.4.11            classInt_0.4-3         
[52] units_0.7-1             grid_4.0.5              unilur_0.4.0.9000      
[55] leafpop_0.1.0           rstudioapi_0.13         htmlwidgets_1.5.3      
[58] crosstalk_1.1.1         labeling_0.4.2          leafem_0.1.6           
[61] base64enc_0.1-3         rmarkdown_2.9.1         gtable_0.3.0           
[64] codetools_0.2-18        DBI_1.1.1               curl_4.3.1             
[67] R6_2.5.0                lubridate_1.7.10        utf8_1.2.1             
[70] KernSmooth_2.23-18      stringi_1.6.2           Rcpp_1.0.6             
[73] vctrs_0.3.8             png_0.1-7               leaflet_2.0.4.1        
[76] tidyselect_1.1.1        xfun_0.24              

---

1. Les iris sont un zonage statistique de l’Insee dont l’acronyme signifie « Ilots Regroupés pour l’Information Statistique ». Leur taille est de 2000 habitants par unité.↩︎