Visualisation des incidents de pipelines au Canada avec R

Comment visualiser les incidents de pipelines canadiens avec R ?

Ce billet explore les données ouvertes du gouvernement du Canada sur les événements de pipeline. On géocode les lieux d’incidents, on cartographie les déversements par produit et quantité, on analyse la distribution par province et compagnie, et on construit un arbre de décision pour classifier le type d’événement.

En 2015, je présentais au DataMeetup de Québec une analyse des données ouvertes sur les incidents de pipelines au Canada. Le code R complet est disponible dans le dépôt — je te le présente ici avec les visualisations.

🌘 Source des données

Les données viennent du gouvernement du Canada : Données sur les événements de pipeline de janvier 2004 au présent. Le jeu contient le détail de chaque événement déclaré par les compagnies pipelinières : type d’incident, produit impliqué, lieu, province, quantités déversées, dommages matériels et environnementaux.

Le dictionnaire de données (dd-20150216.csv) décrit chaque colonne et ses libellés. On y réfère dynamiquement dans le code pour extraire les noms de champs.

🌘 Chargement et préparation

On commence par charger les librairies et lire les fichiers. J’utilise data.table pour la performance et dplyr pour la manipulation.

options(unzip = 'internal')
library("data.table")
library("dplyr")
library("maps")
library("ggmap")
library("reshape2")
library("lubridate")
library("rpart")
library("rpart.plot")

Le dictionnaire contient les métadonnées des colonnes ; les données brutes sont dans p2015-07-fr.csv. On convertit l’encodage en latin1 pour éviter les problèmes de caractères français.

datadict <- read.csv(file = "dd-20150216.csv", header = TRUE,
                     stringsAsFactors = TRUE, encoding = "latin1") %>%
  as.data.table()
donnees <- read.csv(file = "p2015-07-fr.csv", header = TRUE,
                    stringsAsFactors = TRUE) %>%
  as.data.table()

class_donnees <- sapply(donnees, class)
character_columns <- names(donnees)[class_donnees == "character"]
num_columns <- names(donnees)[class_donnees != "character"]

donnees_conv <- cbind(
  donnees[, lapply(.SD, iconv, "utf-8", "latin1"),
          .SDcols = character_columns],
  donnees[, .SD, .SDcols = num_columns])

🌘 Géocodage des incidents

Chaque incident a un lieu (LOCATION) et une province (PROVINCE). Pour les cartographier, on géocode ces adresses avec l’API Google via ggmap::geocode(). Les résultats sont sauvegardés dans geocodedLocations.Rdata pour éviter de frapper l’API à chaque exécution.

LocationsUnique <- unique(donnees_conv[, .(LOCATION, PROVINCE)])
catLocationsUnique <- unique(paste(LocationsUnique$LOCATION,
  ifelse(!is.na(LocationsUnique$PROVINCE), LocationsUnique$PROVINCE, ""),
  "CANADA", sep = " "))
geocodedLocations <- cbind(LocationsUnique, geocode(catLocationsUnique))
save(geocodedLocations, file = "geocodedLocations.Rdata")

On joint les coordonnées aux données d’origine et on filtre les lignes sans province ou lieu.

load("geocodedLocations.Rdata")
donnees_conv2 <- as.data.table(merge(donnees_conv, geocodedLocations,
                                     by = c("LOCATION", "PROVINCE")))
donnees_conv3 <- donnees_conv2[!is.na(PROVINCE) & !is.na(LOCATION)]

🌘 Carte de fond OpenStreetMap

Pour le rendu cartographique, on télécharge un fond de carte OpenStreetMap centré sur l’étendue des incidents.

range_lat <- range(donnees_conv3$lat)
range_lon <- range(donnees_conv3$lon)

openstreetmap1 <- get_openstreetmap(
  bbox = c(left = range_lon[1], bottom = range_lat[1],
           right = range_lon[2], top = range_lat[2]),
  scale = 32500000, format = c("png"),
  crop = TRUE, messaging = FALSE, urlonly = FALSE, color = c("color", "bw"))

save(openstreetmap1, file = "openstreetmap1.Rdata")

Le résultat est mis en cache dans openstreetmap1.Rdata pour accélérer les itérations suivantes.

load("openstreetmap1.Rdata")

🌘 Carte des incidents avec déversement

On filtre les incidents où QUANTITY_LOST > 0 et on superpose les points sur la carte. La couleur indique le type de produit, la taille la quantité déversée.

donnees_incidents <- donnees_conv3[QUANTITY_LOST > 0]
ggmap(openstreetmap1) +
  geom_point(aes(x = lon, y = lat, colour = PRODUCT_TYPE,
                 size = QUANTITY_LOST),
             data = donnees_incidents) +
  scale_size(range = c(2, 8),
    name = datadict[datadict$Nom.de.colonne.PODS == "QUANTITY_LOST"][2]) +
  scale_colour_discrete(
    name = datadict[datadict$Nom.de.colonne.PODS == "PRODUCT_TYPE"][2]) +
  ggtitle(label = "Incidents avec déversement") +
  xlab("Longitude") + ylab("Latitude")

Carte des incidents avec déversement — la couleur indique le type de produit, la taille la quantité déversée

🌘 Distribution des incidents par province et type

Une petite fonction utilitaire pour extraire les initiales des provinces (ex.: « QC » pour Québec), suivie d’une agrégation par province et type d’événement.

initials <- function(str) {
  regex <- '(?<=^|\\s|-)[[:alpha:]]'
  ini <- regmatches(str, gregexpr(regex, str, perl = TRUE))
  toupper(sapply(ini, paste0, collapse = ''))
}

dist_incidents <- donnees_conv3 %>%
  select(PROVINCE, OCC_TYPE, PRODUCT_TYPE) %>%
  group_by(PROVINCE = factor(initials(PROVINCE)),
           OCC_TYPE = factor(OCC_TYPE)) %>%
  summarise(N = n())

Le graphique à barres montre la répartition des types d’incidents par province.

ggplot(dist_incidents, aes(x = PROVINCE, fill = OCC_TYPE)) +
  geom_bar(aes(y = N), stat = "identity") +
  ggtitle(paste0(
    datadict[datadict$Nom.de.colonne.PODS == "OCC_TYPE"][2],
    ", par ",
    datadict[datadict$Nom.de.colonne.PODS == "PROVINCE"][2])) +
  scale_fill_discrete(
    name = datadict[datadict$Nom.de.colonne.PODS == "OCC_TYPE"][2]) +
  xlab(datadict[datadict$Nom.de.colonne.PODS == "PROVINCE"][2]) +
  ylab("Nombre")

Distribution des incidents par province et type d'événement

🌘 Série chronologique : incidents cumulés par compagnie

On prépare une série temporelle des événements cumulés par compagnie pipelinière. La date est nettoyée avec lubridate::parse_date_time().

data_serie <- donnees_conv3 %>%
  select(PIPELINE_CO, OCC_DATE, QUANTITY_LOST) %>%
  mutate(N = 1, date_clean = parse_date_time(OCC_DATE, orders = "%Y-%m-%d")) %>%
  arrange(PIPELINE_CO, date_clean) %>%
  group_by(PIPELINE_CO, date_clean) %>%
  summarise(Nsum = sum(N)) %>%
  mutate(Ncumsum = cumsum(Nsum))

Le graphique linéaire montre l’évolution du nombre d’incidents cumulés pour chaque compagnie.

ggplot(data_serie, aes(x = date_clean, y = Ncumsum, color = PIPELINE_CO)) +
  geom_line(stat = "identity") +
  xlab("Date") +
  ylab("Nombre d'incidents cumulés") +
  ggtitle("Nombre d'incidents cumulés, par compagnie")

Nombre d'incidents cumulés par compagnie

🌘 Arbre de décision

Pour terminer, on construit un arbre de classification qui tente de prédire le type d’événement (OCC_TYPE) à partir de caractéristiques comme le type d’installation, le produit, les dommages environnementaux et matériels, ou la présence de feu ou d’explosion.

tree1 <- rpart(OCC_TYPE ~ FACILITY_TYPE + PROD_RELEASED_IND +
                 PIPELINE_TYPE_CO + ENVIRO_DAMAGE_IND +
                 PROPERTY_DAMAGE_IND + FIRE_IND + EXPLOSION_IND,
               data = donnees_conv3, method = "class")
summary(tree1)

L’affichage de l’arbre donne une vision claire des règles de classification.

  rpart.plot::rpart.plot(tree1,
    main = paste0("Arbre de classification pour ",
      datadict[datadict$Nom.de.colonne.PODS == "OCC_TYPE"]$Nom.du.champ))

Arbre de classification pour le type d'événement

🌘 Version Quarto

Ce document a été converti au format Quarto :

🌘 Pour aller plus loin

Tous les fichiers sources (données, code R, résultats géocodés et carte) sont disponibles dans le dossier static/images/blogue/2019/2019-08-12-visualisation-incidents-pipelines/ de ce site. Tu peux réexécuter l’analyse au complet en ouvrant pipelineViz.Rmd dans RStudio.

Le jeu de données du gouvernement du Canada est mis à jour régulièrement — n’hésite pas à re-télécharger la version la plus récente pour voir l’évolution des incidents.