Feuille d’exercice

Contexte

On va s’intéresser aux données sur les accidents corporels de circulation mises à disposition par le ministère de l’Intérieur, ici en version 2018 : https://www.data.gouv.fr/fr/datasets/base-de-donnees-accidents-corporels-de-la-circulation/

Les objectifs vont être d’importer ses données, de les jointer avec des données sur la population, par département et commune, et de calculer à différents niveaux des indicateurs d’accidentalité.

On représentera ensuite ces indicateurs sur une ou plusieurs cartes qui seront exportées.

Pour cela, il faut télécharger le dossier https://github.com/tvroylandt/r-cartostat qui contient aussi les slides du cours.

Packages

• On charge les packages qui vont nous être utiles.

library(tidyverse)
library(here)
library(sf)
library(readxl)

Import des données

Les données sont disponibles dans le dossier data/dataset/exo.

Elles se composent :

• D’un fichier sur les accidents en 2018 : accidents18_exo.csv ;
• D’un fichier sur la population par département en 2017 : pop_dep17.csv ;
• D’un fichier sur la population par commune en 2017 : pop_com17.xlsx

Données accidents

• Importer le fichier sur les accidents dans un objet appelé df_accidents18.
• Combien de colonnes comporte-t-il ?
• Combien de lignes ?
• Assurez vous que le type des colonnes est cohérent avec leur contenu.
• Sinon, recodez les colonnes qui sont mal importées.
• Visualisez le jeu de données.

# import
df_accidents18 <-
  read_csv2(
    here("data/dataset/exo/accidents18_exo.csv"),
    col_types = cols(
      id_accident = col_character(),
      mois = col_character(),
      code_dep = col_character(),
      code_com = col_character(),
      luminosite = col_character(),
      agglo = col_character(),
      nb_indemnes = col_double(),
      nb_blesses_hospi = col_double(),
      nb_blesses_legers = col_double(),
      nb_tues = col_double(),
      nb_impliques = col_double()
    )
  )

# dimension
dim(df_accidents18)

## [1] 57783    11

# jeu de données
df_accidents18

## # A tibble: 57,783 x 11
##    id_accident mois  code_dep code_com luminosite agglo nb_indemnes nb_blesses_hospi nb_blesses_lege~
##    <chr>       <chr> <chr>    <chr>    <chr>      <chr>       <dbl>            <dbl>            <dbl>
##  1 2018000000~ 1     59       59005    Jour       Hors~           1                1                0
##  2 2018000000~ 2     59       59011    Jour       En a~           1                0                1
##  3 2018000000~ 3     59       59477    Jour       En a~           1                1                0
##  4 2018000000~ 5     59       59052    Jour       En a~           1                1                0
##  5 2018000000~ 6     59       59477    Jour       En a~           1                0                1
##  6 2018000000~ 9     59       59052    Aube/Crép~ En a~           2                1                1
##  7 2018000000~ 9     59       59133    Nuit avec~ En a~           0                1                0
##  8 2018000000~ 11    59       59011    Nuit avec~ En a~           1                0                0
##  9 2018000000~ 2     59       59550    Jour       Hors~           0                0                1
## 10 2018000000~ 3     59       59051    Jour       En a~           1                1                0
## # ... with 57,773 more rows, and 2 more variables: nb_tues <dbl>, nb_impliques <dbl>

Données population

• Importer les données départementales de population dans un objet df_pop_dep17.
• Idem pour les communales dans df_pop_com17.

# départements
df_pop_dep17 <- read_csv2(here("data/dataset/exo/pop_dep17.csv"))

df_pop_dep17

## # A tibble: 100 x 2
##    code_dep pop_tot
##    <chr>      <dbl>
##  1 01        659180
##  2 02        546527
##  3 03        347035
##  4 04        168381
##  5 05        145883
##  6 06       1097496
##  7 07        334688
##  8 08        280032
##  9 09        157210
## 10 10        317118
## # ... with 90 more rows

# communes
df_pop_com17 <- read_xlsx(here("data/dataset/exo/pop_com17.xlsx"))

df_pop_com17

## # A tibble: 34,995 x 2
##    code_com pop_tot
##    <chr>      <dbl>
##  1 01001        794
##  2 01002        249
##  3 01004      14428
##  4 01005       1723
##  5 01006        117
##  6 01007       2841
##  7 01008        767
##  8 01009        339
##  9 01010       1132
## 10 01011        387
## # ... with 34,985 more rows

Bonus - offres

On a construit un jeu de données sur les offres par département : offres_stmt19.sas7bdat.

• Importer ce jeu de données.

library(haven)

df_off_dep19 <-
  read_sas(here("data/dataset/exo/offres_stmt19.sas7bdat"))

df_off_dep19 

## # A tibble: 100 x 2
##    DPTETA nb_off
##    <chr>   <dbl>
##  1 73      41954
##  2 28      12814
##  3 60      26753
##  4 52       6105
##  5 89      12723
##  6 01      26490
##  7 77      42081
##  8 84      33083
##  9 27      18903
## 10 45      34033
## # ... with 90 more rows

Calculs départementaux

On va regarder ce que ça donne au niveau des départements.

Regroupement

• Regrouper les données par département.
• Pour cela on somme les indicateurs.
• Calculer aussi un nombre d’accidents.

df_accidents_dep18 <- df_accidents18 %>% 
  group_by(code_dep) %>% 
  summarise(nb_accidents = n(),
            nb_indemnes = sum(nb_indemnes),
            nb_blesses_legers = sum(nb_blesses_legers),
            nb_blesses_hospi = sum(nb_blesses_hospi),
            nb_tues = sum(nb_tues),
            nb_impliques = sum(nb_impliques)) %>% 
  ungroup()

df_accidents_dep18

## # A tibble: 101 x 7
##    code_dep nb_accidents nb_indemnes nb_blesses_legers nb_blesses_hospi nb_tues nb_impliques
##    <chr>           <int>       <dbl>             <dbl>            <dbl>   <dbl>        <dbl>
##  1 01                454         394               373              275      36         1078
##  2 02                261         242               173              187      31          633
##  3 03                253         213               158              154      25          550
##  4 04                204         160               155              131      19          465
##  5 05                241         242               247               89      19          597
##  6 06               1156        1066               957              432      57         2512
##  7 07                219         168               145              131      12          456
##  8 08                124         101                72               87      12          272
##  9 09                 90          82                54               68      12          216
## 10 10                327         282               301              115      18          716
## # ... with 91 more rows

Jointure

• Jointer sur les données de population.

df_accidents_dep_pop <- df_accidents_dep18 %>% 
  left_join(df_pop_dep17, by = "code_dep")

df_accidents_dep_pop

## # A tibble: 101 x 8
##    code_dep nb_accidents nb_indemnes nb_blesses_legers nb_blesses_hospi nb_tues nb_impliques pop_tot
##    <chr>           <int>       <dbl>             <dbl>            <dbl>   <dbl>        <dbl>   <dbl>
##  1 01                454         394               373              275      36         1078  659180
##  2 02                261         242               173              187      31          633  546527
##  3 03                253         213               158              154      25          550  347035
##  4 04                204         160               155              131      19          465  168381
##  5 05                241         242               247               89      19          597  145883
##  6 06               1156        1066               957              432      57         2512 1097496
##  7 07                219         168               145              131      12          456  334688
##  8 08                124         101                72               87      12          272  280032
##  9 09                 90          82                54               68      12          216  157210
## 10 10                327         282               301              115      18          716  317118
## # ... with 91 more rows

Calculer des indicateurs

• Calculer le nombre d’accidents pour 10 000 habitants pour chaque département.
• Calculer le nombre de tués pour 10 000 habitants pour chaque département.
• Découper ces deux parts en classes (à choisir manuellement).
• Ne conserver que les variables d’intérêt.
• Enlever les départements à blanc (20 et 976).

df_accidents_dep_pop_r <- df_accidents_dep_pop %>%
  mutate(
    accidents_pop = nb_accidents / pop_tot * 10000,
    tues_pop = nb_tues / pop_tot * 10000,
    accidents_pop_cut = cut(
      accidents_pop,
      include.lowest = TRUE,
      breaks = c(0, 5, 10, 15, 10000)
    ),
    tues_pop_cut = cut(
      tues_pop,
      include.lowest = TRUE,
      breaks = c(0, 0.3, 0.6, 1, 10000),
      labels = c("0,3 et moins", "De 0,3 à 0,6", "De 0,6 à 1", "Plus de 1")
    )
  ) %>%
  select(code_dep,
         nb_accidents,
         accidents_pop,
         accidents_pop_cut,
         tues_pop,
         tues_pop_cut) %>%
  filter(!is.na(tues_pop))

df_accidents_dep_pop_r

## # A tibble: 100 x 6
##    code_dep nb_accidents accidents_pop accidents_pop_cut tues_pop tues_pop_cut
##    <chr>           <int>         <dbl> <fct>                <dbl> <fct>       
##  1 01                454          6.89 (5,10]               0.546 De 0,3 à 0,6
##  2 02                261          4.78 [0,5]                0.567 De 0,3 à 0,6
##  3 03                253          7.29 (5,10]               0.720 De 0,6 à 1  
##  4 04                204         12.1  (10,15]              1.13  Plus de 1   
##  5 05                241         16.5  (15,1e+04]           1.30  Plus de 1   
##  6 06               1156         10.5  (10,15]              0.519 De 0,3 à 0,6
##  7 07                219          6.54 (5,10]               0.359 De 0,3 à 0,6
##  8 08                124          4.43 [0,5]                0.429 De 0,3 à 0,6
##  9 09                 90          5.72 (5,10]               0.763 De 0,6 à 1  
## 10 10                327         10.3  (10,15]              0.568 De 0,3 à 0,6
## # ... with 90 more rows

# pour choisir les classes
ggplot(df_accidents_dep_pop_r,
       aes(x = accidents_pop)) +
  geom_histogram()

ggplot(df_accidents_dep_pop_r,
       aes(x = tues_pop)) +
  geom_histogram()

• Refaire les opérations en conservant la variable mois.

df_accidents_dep_pop_mois_r <- df_accidents18 %>%
  group_by(code_dep, mois) %>%
  summarise(
    nb_accidents = n(),
    nb_indemnes = sum(nb_indemnes),
    nb_blesses_legers = sum(nb_blesses_legers),
    nb_blesses_hospi = sum(nb_blesses_hospi),
    nb_tues = sum(nb_tues),
    nb_impliques = sum(nb_impliques)
  ) %>%
  ungroup() %>%
  left_join(df_pop_dep17, by = "code_dep") %>%
  mutate(
    accidents_pop = nb_accidents / pop_tot * 10000,
    tues_pop = nb_tues / pop_tot * 10000,
    accidents_pop_cut = cut(
      accidents_pop,
      include.lowest = TRUE,
      breaks = c(0, 5, 10, 15, 10000)
    ),
    tues_pop_cut = cut(
      tues_pop,
      include.lowest = TRUE,
      breaks = c(0, 0.3, 0.6, 1, 10000),
      labels = c("0,3 et moins", "De 0,3 à 0,6", "De 0,6 à 1", "Plus de 1")
    )
  ) %>%
  select(code_dep,
         mois,
         nb_accidents,
         accidents_pop,
         accidents_pop_cut,
         tues_pop,
         tues_pop_cut) %>%
  filter(!is.na(tues_pop))

df_accidents_dep_pop_mois_r

## # A tibble: 1,199 x 7
##    code_dep mois  nb_accidents accidents_pop accidents_pop_cut tues_pop tues_pop_cut
##    <chr>    <chr>        <int>         <dbl> <fct>                <dbl> <fct>       
##  1 01       1               32         0.485 [0,5]               0.0455 0,3 et moins
##  2 01       10              46         0.698 [0,5]               0      0,3 et moins
##  3 01       11              37         0.561 [0,5]               0.0455 0,3 et moins
##  4 01       12              36         0.546 [0,5]               0.0303 0,3 et moins
##  5 01       2               22         0.334 [0,5]               0.0759 0,3 et moins
##  6 01       3               41         0.622 [0,5]               0.0455 0,3 et moins
##  7 01       4               33         0.501 [0,5]               0.0152 0,3 et moins
##  8 01       5               42         0.637 [0,5]               0.0152 0,3 et moins
##  9 01       6               49         0.743 [0,5]               0.137  0,3 et moins
## 10 01       7               37         0.561 [0,5]               0.0759 0,3 et moins
## # ... with 1,189 more rows

Luminosité

• Calculer une part d’accident qui ont lieu le jour, par département.
• Découper cette part en classes selon les quantiles.
• Découper cette part en classes définies manuellement.
• Garder code_dep et les variables de pourcentages.

df_accidents_dep_jour18 <- df_accidents18 %>%
  group_by(luminosite, code_dep) %>%
  count() %>%
  group_by(code_dep) %>%
  mutate(perc_jour = n / sum(n) * 100,
         perc_jour = round(perc_jour, 1)) %>%
  filter(luminosite == "Jour") %>%
  ungroup() %>%
  mutate(
    perc_jour_quantile = cut(
      perc_jour,
      breaks = quantile(.$perc_jour),
      include.lowest = TRUE
    ),
    perc_jour_cut = cut(
      perc_jour,
      include.lowest = TRUE,
      breaks = c(0, 70, 75, 80, 100)
    )
  ) %>% 
  select(code_dep, starts_with("perc_jour"))

df_accidents_dep_jour18

## # A tibble: 101 x 4
##    code_dep perc_jour perc_jour_quantile perc_jour_cut
##    <chr>        <dbl> <fct>              <fct>        
##  1 01            67.4 (65,67.8]          [0,70]       
##  2 02            64.4 [56.3,65]          [0,70]       
##  3 03            71.1 (70.6,76.5]        (70,75]      
##  4 04            76.5 (70.6,76.5]        (75,80]      
##  5 05            76.3 (70.6,76.5]        (75,80]      
##  6 06            74.2 (70.6,76.5]        (70,75]      
##  7 07            74.4 (70.6,76.5]        (70,75]      
##  8 08            61.3 [56.3,65]          [0,70]       
##  9 09            72.2 (70.6,76.5]        (70,75]      
## 10 10            68.8 (67.8,70.6]        [0,70]       
## # ... with 91 more rows

Bonus - Agglomération

• Construire un score de gravité pour chaque accident :

  • Indemne = 0.1 ;
  • Blessés légers = 1 ;
  • Blessés hospitalisés = 2 ;
  • Tués = 5.

• Calculer le score moyen dans et hors agglomération, par département.

df_accidents_score <- df_accidents18 %>% 
  mutate(score_acc = 0.1 * nb_indemnes + 1 * nb_blesses_legers + 2 * nb_blesses_hospi + 5 * nb_tues) %>% 
  group_by(agglo, code_dep) %>% 
  summarise(score_acc_mean = mean(score_acc)) %>% 
  ungroup() 

df_accidents_score

## # A tibble: 202 x 3
##    agglo            code_dep score_acc_mean
##    <chr>            <chr>             <dbl>
##  1 En agglomération 01                 2.15
##  2 En agglomération 02                 2.08
##  3 En agglomération 03                 2.13
##  4 En agglomération 04                 2.30
##  5 En agglomération 05                 1.76
##  6 En agglomération 06                 1.82
##  7 En agglomération 07                 1.89
##  8 En agglomération 08                 2.02
##  9 En agglomération 09                 2.48
## 10 En agglomération 10                 1.70
## # ... with 192 more rows

• Calculer la différence de score entre dans et hors agglomération, par département (utiliser pivot_wider puis pivot_longer).

df_accidents_score_diff <- df_accidents_score %>% 
  pivot_wider(names_from = agglo,
              values_from = score_acc_mean) %>% 
  mutate(diff_score = `Hors agglomération` - `En agglomération`)

df_accidents_score_diff

## # A tibble: 101 x 4
##    code_dep `En agglomération` `Hors agglomération` diff_score
##    <chr>                 <dbl>                <dbl>      <dbl>
##  1 01                     2.15                 2.77      0.622
##  2 02                     2.08                 3.32      1.23 
##  3 03                     2.13                 2.66      0.526
##  4 04                     2.30                 2.69      0.384
##  5 05                     1.76                 2.87      1.11 
##  6 06                     1.82                 2.33      0.516
##  7 07                     1.89                 2.38      0.481
##  8 08                     2.02                 3.28      1.26 
##  9 09                     2.48                 3.07      0.592
## 10 10                     1.70                 2.68      0.979
## # ... with 91 more rows

Recoupement avec les données spatiales

Chargement

• Charger le fond de carte départemental qui se trouve dans data/shp.

shp_dep <- st_read(here("data/shp/FR_DEP_DOM_IDF.shp"), quiet = TRUE)

Jointure

• On jointe les différentes données avec les fichiers spatiaux :

  • par département ;
  • par département/mois ;
  • pour les accidents de jour ;
  • pour le score agglomération/hors agglomération.

df_accidents_dep_shp <- shp_dep %>% 
  inner_join(df_accidents_dep_pop_r, by = "code_dep")

df_accidents_dep_mois_shp <- shp_dep %>% 
  inner_join(df_accidents_dep_pop_mois_r, by = "code_dep")

df_accidents_dep_jour_shp <- shp_dep %>% 
  inner_join(df_accidents_dep_jour18, by = "code_dep")

df_accidents_dep_agglo_score_shp <- shp_dep %>% 
  inner_join(df_accidents_score, by = "code_dep")

Représentation cartographique

Aplats de couleur

• Représenter la part des accidents qui arrivent le jour en continu.
• Mettre un thème blanc et changer la légende.

ggplot(data = df_accidents_dep_jour_shp) +
  geom_sf(aes(fill = perc_jour)) +
  scale_fill_viridis_c(name = "Part d'accidents de jour", 
                       direction = -1) +
  theme_void()

• Faire de même en zoomant sur les Hauts-de-France
• Qu’observe-t-on sur la légende ?

ggplot(data = df_accidents_dep_jour_shp %>%  filter(code_reg == "32")) +
  geom_sf(aes(fill = perc_jour)) +
  scale_fill_viridis_c(name = "Part d'accidents de jour",
                       direction = -1) +
  theme_void()

• Faire de même avec la variable découpée en quantile.

ggplot(data = df_accidents_dep_jour_shp) +
  geom_sf(aes(fill = perc_jour_quantile)) +
  scale_fill_viridis_d(name = "Part d'accidents de jour", 
                       direction = -1) +
  theme_void()

Ajouter le texte

• Ajouter la valeur, arrondie à l’unité, sur les départements.
• Si possibles, pour l’Île-de-France, ne le faire que sur le zoom.

ggplot() +
  geom_sf(data = df_accidents_dep_jour_shp, aes(fill = perc_jour_quantile)) +
  geom_sf_text(
    data = df_accidents_dep_jour_shp %>% filter(zoom_idf == "1" |
                                                  code_reg != "11"),
    aes(label = round(perc_jour, 0)),
    fontface = "bold",
    color = "black"
  ) +
  scale_fill_viridis_d(name = "Part d'accidents de jour",
                       direction = -1) +
  theme_void()

Région

• Charger le fond de carte régions.
• Faire apparaître en plus gros les contours des régions.

shp_reg <- st_read(here("data/shp/FR_REG_DOM.shp"))

## Reading layer `FR_REG_DOM' from data source `D:\r-cartostat\data\shp\FR_REG_DOM.shp' using driver `ESRI Shapefile'
## Simple feature collection with 18 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 99787.91 ymin: 5990209 xmax: 1242111 ymax: 7109552
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

ggplot(data = df_accidents_dep_jour_shp) +
  geom_sf(aes(fill = perc_jour_quantile)) +
  geom_sf(
    data = shp_reg,
    alpha = 0,
    color = "black",
    size = 1.2
  ) +
  scale_fill_viridis_d(name = "Part d'accidents de jour",
                       direction = -1) +
  theme_void()

Cercle proportionnels

• Représenter le nombre d’accidents par des cercles proportionnels.
• En colorer l’intérieur par la part d’accidents par habitant, en continu.

ggplot(df_accidents_dep_shp) +
  geom_sf() +
  stat_sf_coordinates(aes(size = nb_accidents, fill = accidents_pop), shape = 21) +
  scale_fill_viridis_c(name = "Ratio d'accidents \npar habitant", direction = -1) +
  scale_size_continuous(name = "Nombre d'accidents", range = c(1, 10)) +
  theme_void()

• En colorer l’intérieur par la part de tués par habitant, en classes.

ggplot(df_accidents_dep_shp) +
  geom_sf() +
  stat_sf_coordinates(aes(size = nb_accidents, fill = tues_pop_cut), shape = 21) +
  scale_fill_viridis_d(name = "Ratio de tués \npar habitant", guide = guide_legend(override.aes = list(size =
                                                                                              5))) +
  scale_size_continuous(name = "Nombre d'accidents", range = c(1, 10)) +
  theme_void()

Effet de la saisonnalité

• Faire une carte du nombre d’accidents par habitatn par mois et département.
• On peut facetter avec facet_wrap ou facet_grid.

df_accidents_dep_mois_shp %>%
  filter(!is.na(mois)) %>%
  mutate(mois = fct_relevel(mois,
                            "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12")) %>% 
  ggplot() +
  geom_sf(aes(fill = accidents_pop)) +
  scale_fill_viridis_c(name = "Nombre d'accidents \npour 10 000 habitants",
                       direction = -1) +
  facet_wrap(vars(mois)) +
  theme_void()

Bonus - Agglomération

• Réprésenter le score moyen en et hors agglomération, avec une facette.

ggplot(df_accidents_dep_agglo_score_shp) +
  geom_sf(aes(fill = score_acc_mean)) +
  scale_fill_viridis_c(direction = -1) +
  facet_wrap(vars(agglo)) +
  theme_void()

Niveau bassin d’emploi

On se penche maintenant sur le niveau bassin d’emploi.

Regroupement en bassins d’emploi

• Importer la table de passage commune -> bassin d’emploi.

passage_com_bassin <- read_xlsx(here("data/passage/passage_commune_bassin_bmo.xlsx"))

• Jointer sur les accidents pour construire le niveau bassin.
• Construire le nombre d’accidents par bassin.

df_accidents_bassin <- df_accidents18 %>% 
  inner_join(passage_com_bassin, by = c("code_com" = "code_commune")) %>% 
  group_by(code_bassin_bmo) %>% 
  count(name = "nb_accidents")

• A-t-on des trous ?
• Sur quelles communes ? Regarder d’où ça vient sur le site de l’Insee.

trous_bassin <- df_accidents18 %>% 
  anti_join(passage_com_bassin, by = c("code_com" = "code_commune"))

trous_bassin %>% 
  group_by(code_com) %>% 
  count()

## # A tibble: 7 x 2
## # Groups:   code_com [7]
##   code_com     n
##   <chr>    <int>
## 1 50439        1
## 2 59248        1
## 3 59540        5
## 4 61377        2
## 5 71260        3
## 6 71511        1
## 7 72318        1

• Importer le fond de carte au niveau bassin d’emploi.

shp_bassin <- st_read(here("data/shp/FR_BASSIN_BMO_DOM_IDF_2019.shp"), quiet = TRUE)

Cartes

• Jointer sur le fond de carte.
• Réaliser une carte du nombre d’accidents par bassin d’emploi, en cercles proportionnels.
• Avec en surimpression les régions.

shp_bassin %>%
  inner_join(df_accidents_bassin, by = c("cd_bss_" = "code_bassin_bmo")) %>%
  ggplot() +
  geom_sf() +
  stat_sf_coordinates(aes(size = nb_accidents, fill = nb_accidents), shape = 21) +
  geom_sf(
    data = shp_reg,
    alpha = 0,
    color = "black",
    size = 1.1
  ) +
  scale_fill_viridis_c(name = "Nombre d'accidents",
                       direction = -1) +
  scale_size_continuous(name = "", range = c(1, 10)) +
  theme_void()

Export

• Exporter la dernière carte.

ggsave(
  "accidents_bassin.png",
  width = 10,
  height = 8,
  units = "cm"
)