J’ai envie de savoir qui a eu le plus souvent besoin de souper aux chandelles à cause d’une panne de courant.
library(dplyr)
library(sf)
library(lubridate)
library(purrr)
library(furrr)
library(DBI)
library(pool)
library(mapview)
library(ggplot2)
library(ggmap)
library(gganimate)
con <- dbPool(
RPostgres::Postgres(),
host = "192.168.2.15",
dbname = Sys.getenv("POSTGIS_DBNAME"),
port = 5432,
user = Sys.getenv("POSTGIS_USER"),
password = Sys.getenv("POSTGIS_PASSWORD"),
bigint = "numeric"
)Hydro-Québec rend les polygones derrières les données d’info-panne disponible sur son site.
Les fichiers KMZ d’info-panne sont disponibles à l’adresse http://pannes.hydroquebec.com/pannes/donnees/v3_0/open_data/.
Si vous le visitez aujourd’hui, vous allez voir les polygones des pannes des 24 dernières heures:
J’imagine que c’est pratique pour répondre à des questions existentielles comme “est-ce que j’ai de l’électricité présentement dans ma maison?” quand on n’a aucun appareil électrique pour vérifier.
Ça ne permet pas de reculer dans le temps pour faire une bonne analyse.
C’est décevant, mais le Simon du passé a une bonne surprise: j’ai commencé à archiver ces fichiers il y a 3 ans.
Aujourd’hui c’est enfin le moment de trouver qui est le pauvre hère qui a manqué le plus longtemps de courant durant ces 3 années!
Voici à quoi ressemble le contenu d’un d’entre eux, celui du 6 avril 2023 à 10h40 du matin (20230406104016.kmz):
url <- "https://raw.githubusercontent.com/simoncoulombe/hydroquebecoutage/main/data/20230406104016.kmz"
# Define local destination
destfile <- "20230406104016.kmz"
# Download the file
#download.file(url, destfile, mode = "wb")
# Function to read KMZ files with better error handling
read_kmz_manual <- function(path, add_metadata = FALSE) {
tryCatch({
# Create unique temporary directory for each file
tmpdir <- file.path(tempdir(), basename(path))
dir.create(tmpdir, showWarnings = FALSE, recursive = TRUE)
# Use utils::unzip instead of zip package - more robust
utils::unzip(path, exdir = tmpdir, junkpaths = TRUE)
# Find KML file
kml_files <- list.files(tmpdir, pattern = "\\.kml$", full.names = TRUE)
if (length(kml_files) == 0) {
stop("No KML file found in KMZ archive")
}
# Read the first KML file
sf_obj <- sf::st_read(kml_files[1], quiet = TRUE)
# Clean up temporary directory
unlink(tmpdir, recursive = TRUE)
# Remove Z/M dimensions unconditionally
sf_obj <- sf::st_zm(sf_obj, drop = TRUE, what = "ZM")
# Check if the file is empty (no features)
if (nrow(sf_obj) == 0) {
cat("KMZ file", basename(path), "is empty (0 features). Skipping.\n")
return(NULL)
}
# Add metadata
if(add_metadata){
sf_obj$source_file <- basename(path)
sf_obj$file_path <- path
sf_obj$row_number <- 1:nrow(sf_obj)
}
return(sf_obj)
}, error = function(e) {
cat("Failed to read", basename(path), ":", conditionMessage(e), "\n")
# Clean up on error
if (exists("tmpdir") && dir.exists(tmpdir)) {
unlink(tmpdir, recursive = TRUE)
}
return(NULL)
})
}
kmz20230406104016 <-read_kmz_manual("20230406104016.kmz")Failed to read 20230406104016.kmz : No KML file found in KMZ archive kmz20230406104016 %>% knitr::kable()Ce sont seulement des polygônes. On n’a pas d’information comme la durée de la panne ou la durée de la panne.
Voici à quoi il peut ressembler sur une carte. Il y a quand même pas mal de pannes. C’est fait exprès, j’ai choisi le moment où le plus de bâtiments étaient affectés, soit 742 955 bâtiments (sachant que les bâtiments collés se font fusionner..)
mapview(kmz20230406104016)NULLMon plan initial était de trouver qui avait manqué le “plus longtemps” de courant. Mais on n’a pas ça dans le polygone. Je pensais initialement faire une hypothèse: les pannes indiquées dans un fichier dureraient du timestamp du fichier jusqu’au timestamp du fichier suivant.
C’est quand même une hypothèse forte: ça implique que le site d’hydro-québec est mis à jour dès qu’une panne se termine et aussi que mon système broche-à-foin récupère tous les KMZ qu’hydro rend disponible.
J’ai donc décidé de compter le nombre de jours avec au moins une panne.
J’ai quand même créé les colonnes “outage_start” , “outage_end” et “outage_duration_sec”, mais je ne les ai pas utilisées.
Mon processus d’ingestion de données a simplement été de faire une boucle sur les 39 000 KMZ, d’ajouter un peu de metadata aux polygones ( source_file, source_path et row_number) et d’append ça à la table “outages” sur mon serveur postgis.
Ensuite j’ai lu chacun des KMZ, j’ai ajouté un peu de metadata aux pannes (outage_start, outage_end, outage_duration_sec, row_num, file_name) et j’ai ajouté ça à mon serveur PostGIS.
Au final, ma table de pannes compte 2 641 483 polygones de pannes, entre le 31 mai 2022 et le 27 août 2025 (1184 jours?):
#dbGetQuery(con, "SELECT count(*) as count , min(source_file) as min_date, max(source_file) as max_date FROM outages ;") %>% gt::gt()st_read(con,query = " select * FROM outages limit 5 ;") %>% knitr::kable(caption = "exemple de lignes dans outages")| Name | Description | source_file | file_path | row_number | outage_start | outage_end | outage_duration_sec | geom |
|---|---|---|---|---|---|---|---|---|
| 20230402220003.kmz | data/20230402220003.kmz | 4 | 2023-04-02 22:00:03 | 2023-04-02 22:10:05 | 602 | POLYGON ((-72.84491 45.2916… | ||
| 20230405114008.kmz | data/20230405114008.kmz | 123 | 2023-04-05 11:40:08 | 2023-04-05 11:50:10 | 602 | POLYGON ((-75.88507 45.4538… | ||
| 20230405114008.kmz | data/20230405114008.kmz | 124 | 2023-04-05 11:40:08 | 2023-04-05 11:50:10 | 602 | POLYGON ((-73.70173 45.4722… | ||
| 20230405114008.kmz | data/20230405114008.kmz | 125 | 2023-04-05 11:40:08 | 2023-04-05 11:50:10 | 602 | POLYGON ((-74.64625 45.7582… | ||
| 20230405114008.kmz | data/20230405114008.kmz | 126 | 2023-04-05 11:40:08 | 2023-04-05 11:50:10 | 602 | POLYGON ((-76.34218 45.5549… |
J’aurais voulu des données ouvertes avec toutes les adresses du Québec. Il me semble que ça existait, mais je n’ai pas trouvé dans mon google rapide.
Je me suis donc rabattu sur le référentiel québécois sur les bâtiments , qui contient les polygônes de tous les bâtiments du Québec. Ça semble avoir été créé à partir d’orthophotographie et d’une bone dose d’intelligence artificielle.
Ça a un gros défaut: c’est plein de beau polygones, mais on n’a pas d’adresse et surtout un jumelé (ou plusieurs maison de villes sur les rues de Montréal) va apparaître comme un seul “bâtiment”.
La table du référentiel bâti compte 3 289 166 bâtiments.
#dbGetQuery(con, "SELECT count(*) FROM batiments_4326 ;") %>% gt::gt()st_read(con,query =" select * FROM batiments_4326 limit 5 ;") %>% knitr::kable(caption = "exemple de lignes dans batiments_4326")| id | idbati | noseq | version | datecreation | methoprod | producteur | source | datesource | altmin | altmax | altmoy | niveaucompl | perimetre | superficie | geom | geom_inside |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 208665 | NA | 953619 | BATI20250401 | 2023-03-31 | Extraction automatique | Ressources naturelles Canada en partenariat avec le gouvernement du Québec | Lidar | 2015-05-31 | NA | NA | NA | NC-0 | 57.70 | 172.63 | POLYGON ((-75.83341 45.3817… | POINT (-75.83329 45.38168) |
| 208666 | NA | 953630 | BATI20250401 | 2023-03-31 | Extraction automatique | Ressources naturelles Canada en partenariat avec le gouvernement du Québec | Lidar | 2015-05-31 | NA | NA | NA | NC-0 | 60.06 | 169.71 | POLYGON ((-75.40774 45.5455… | POINT (-75.40766 45.54556) |
| 208667 | NA | 953636 | BATI20250401 | 2023-03-31 | Extraction automatique | Ressources naturelles Canada en partenariat avec le gouvernement du Québec | Lidar | 2014-10-22 | NA | NA | NA | NC-0 | 36.79 | 84.37 | POLYGON ((-77.27872 45.9947… | POINT (-77.27865 45.99476) |
| 208668 | NA | 953647 | BATI20250401 | 2023-03-31 | Extraction automatique | Ressources naturelles Canada en partenariat avec le gouvernement du Québec | Lidar | 2014-10-22 | NA | NA | NA | NC-0 | 28.53 | 50.20 | POLYGON ((-77.07098 45.9107… | POINT (-77.07092 45.91071) |
| 208669 | NA | 953651 | BATI20250401 | 2023-03-31 | Extraction automatique | Ressources naturelles Canada en partenariat avec le gouvernement du Québec | Lidar | 2019-10-16 | NA | NA | NA | NC-0 | 21.15 | 30.07 | POLYGON ((-68.46981 48.4998… | POINT (-68.46978 48.49983) |
J’ai commencé par rouler cette query qui regarde quel bâtiment intersecte avec quelle panne. Ça devait me redonner une ligne par bâtiment-panne:
buildings_sf <- st_read(con, query = "
+ SELECT DISTINCT b.*, o.source_file, o.outage_start, o.outage_end, o.outage_duration_sec
+ FROM batiments_4326 AS b
+ JOIN outages AS o
+ ON ST_Intersects(b.geom, o.geom)Mais après un très long moment j’ai dû me rentre à l’évidence que je devais simplifier car trouver l’intersection de 2.6 millions de pannes avec 3.2 millions de bâtiments ça va être long.
J’ai donc créé une nouvelle colonne de points pour les centroïdes de chaque bâtiments et je regarde ensuite si le centroide est dans une panne. C’est beaucoup moins lourd que de regarder si “une partie du building” est dans la panne
Aussi, j’ai pas besoin de retourner la geometry des batiments à chaque ligne, alors j’ai flushé ça.
dbExecute(con, "
ALTER TABLE batiments_4326
ADD COLUMN geom_inside geometry(Point, 4326);
")
dbExecute(con, "
UPDATE batiments_4326
SET geom_inside = ST_PointOnSurface(geom);
")
dbExecute(con, "
CREATE INDEX IF NOT EXISTS batiments_inside_idx
ON batiments_4326 USING GIST(geom_inside);
")
dbExecute(con, "
CREATE TABLE points_in_outages AS
SELECT
b.id AS building_id,
o.source_file,
o.row_number,
o.outage_duration_sec,
CAST(LEFT(regexp_replace(split_part(o.source_file, '/', -1), '\\.kmz$', ''), 8) AS date) AS outage_date
FROM batiments_4326 AS b
JOIN outages AS o
ON ST_Contains(o.geom, b.geom_inside);
")La table points_in_outage qui indique simplement chaque combinaison de panne-bâtiments compte un débile 284 841 906 lignes et a pris une vingtaine de minutes minutes à calculer sur mon ryzen 5600g.
#dbGetQuery(con, "select count(*) as count from points_in_outages")La table ressemble à ceci :
dbGetQuery(con, "select * from points_in_outages limit 5") building_id source_file row_number outage_duration_sec outage_date
1 269732 20230406215005.kmz 2283 602 2023-04-06
2 260865 20230406215005.kmz 2283 602 2023-04-06
3 77607 20230406215005.kmz 2283 602 2023-04-06
4 442280 20230406215005.kmz 2283 602 2023-04-06
5 71756 20230406215005.kmz 2283 602 2023-04-06C’est normal que le nombre de ligne soit élevé: une grosse panne qui affecterait le québec en entier affecterait 3 millions de bâtiments et contriburait autant de lignes à la table.
Tiens, c’est quoi le moment avec le plus de “bâtiments” affectés ? voici la distribution du nombre de bâtiments affectés par fichier de panne
# tictoc::tic()
# dbExecute(con, "
# CREATE TABLE buildings_per_file AS
# SELECT source_file, COUNT(DISTINCT building_id) AS building_count
# FROM points_in_outages
# GROUP BY source_file;"
# )
# tictoc::toc() # 295s elapsed # 38944 rows
buildings_per_file <- dbGetQuery(con, "select * from buildings_per_file")C’était le 6 avril juillet 2023 au matin, avec 742955 “buildings” affectés.
buildings_per_file %>% filter(building_count == max(building_count)) %>% head(1) source_file building_count
1 20230406104016.kmz 742955Ça ressemblait à ça sur une carte:
pire_panne_poly <- st_read(con, query = "select * from outages where source_file = '20230406104016.kmz' ")
mapview(pire_panne_poly)
mettons que je veux les polygones des 745 000 buildings, je ferais qqch comme
# dbExecute(con, "CREATE TEMP TABLE peak_ids AS
# SELECT DISTINCT building_id
# FROM points_in_outages
# WHERE source_file = '20230406104016.kmz';")
#
# peak_buildings <- st_read(con, query= "SELECT b.id, b.geom
# FROM batiments_4326 AS b
# JOIN peak_ids AS p
# ON b.id = p.building_id;")On va commencer par créer une table
# one row per building per day, even if a building appears in multiple polygons that day.
dbExecute(con, "
CREATE TABLE buildings_in_day AS
SELECT DISTINCT
building_id,
outage_date
FROM points_in_outages;
") # 26758624
# Number of buildings affected per day (outage_date × total buildings). To find the worst day
dbExecute(con, "
CREATE TABLE buildings_count_per_day AS
SELECT
outage_date,
COUNT(*) AS building_count
FROM buildings_in_day
GROUP BY outage_date
ORDER BY outage_date;
") # 1168
# Number of outage days per building per month for each building (building_id × month × days_with_outage)
dbExecute(con, "
CREATE TABLE building_outage_days_per_month AS
SELECT
building_id,
DATE_TRUNC('month', outage_date)::date AS month,
COUNT(*) AS days_with_outage
FROM buildings_in_day
GROUP BY building_id, month
ORDER BY building_id, month;
")
dbExecute(con, "
CREATE TABLE building_outage_ever AS
SELECT
building_id,
COUNT(*) AS days_with_outage
FROM buildings_in_day
GROUP BY building_id
ORDER BY building_id;
") # 3100871Bon, je pensais faire un beau geocoding pour sortir les adresses, mais on va se contenter d’une carte:
top10_pauvres_heres <- st_read(con, query = "
SELECT b.id, b.geom_inside, top10.days_with_outage
FROM batiments_4326 AS b
JOIN (
SELECT building_id, days_with_outage
FROM building_outage_ever
ORDER BY days_with_outage DESC
LIMIT 100
) AS top10
ON b.id = top10.building_id;
")top10_pauvres_heres %>% st_drop_geometry() %>% head(10) %>% knitr::kable()| id | days_with_outage |
|---|---|
| 429738 | 88 |
| 39955 | 88 |
| 250624 | 86 |
| 1565816 | 68 |
| 1578947 | 68 |
| 2151690 | 67 |
| 2137128 | 67 |
| 2131254 | 67 |
| 1592203 | 67 |
| 2154164 | 67 |
mapview(top10_pauvres_heres)
Avec qq screenshots zoomés.
Voici le champion toute catégorie, le nu,éro 429738 (et sa grappe), sur l’avenue Boomfield à Montréal, qui a eu des pannes pendant 88 journées distinctes.
Ceci-dit, j’ai envie de l’accuser de dopage, parce que lui et sa grappe n’aurait pas eu de courant pendant tout l’été 2023 alors je pense que c’était peut-être fait exprès?
La 2e plus grosse grappe c’est un paquet de maisons au bord du Lac des Seize-Iles.
Ce pauvre type a perdu le courant 68 jours au cours des 3 dernières années. J’avoue que j’ai envie de le déclarer champion car il a accumulé un peu de pannes à tous les mois.
Voici la distribution de leur nombre de pannes mensuel:
top2 <- dbGetQuery(con, "select * from building_outage_days_per_month where building_id in ( 429738, 1565816) ")top2 <- top2 %>%
mutate(building_name = case_when(
building_id == 429738 ~ "Montreal",
building_id == 1565816 ~ "Lac des Seize Îles",
TRUE ~ as.character(building_id)
))
myplot <- top2 %>%
ggplot(aes(x = month, y = days_with_outage, fill = building_name)) +
geom_col(show.legend = FALSE) +
geom_text(aes(label = days_with_outage), vjust = -0.5, size = 3) +
facet_wrap(~building_name, ncol = 1) +
scale_x_date(
date_labels = "%b %Y",
date_breaks = "1 month",
limits = as.Date(c("2022-04-01", "2025-09-01"))
) +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
labs(
title = "Distribution du nombre mensuel de pannes pour le top 2",
x = NULL,
y = "Nombre de jours avec au moins une panne"
) +
theme_minimal(base_size = 12) +
theme(
strip.text = element_text(face = "bold", size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid.major.x = element_blank()
)
ggsave("thumbnail.png", plot = myplot, width = 8, height = 5)
myplot
Pour rire, une petite animation rapide des tourments de notre maison aux Seize Iles. L’idée c’est de montrer que je ne suis pas buggé et qu’il a bien eu plein de pannes différentes…
center_point <- st_read(
con,
query = "
SELECT *
FROM batiments_4326
WHERE id = 1565816
"
)
outages <- st_read(
con,
query = "
SELECT o.source_file, o.row_number, o.geom,
CAST(LEFT(regexp_replace(split_part(o.source_file, '/', -1), '\\.kmz$', ''), 8) AS date) AS outage_date
FROM outages o
JOIN (
SELECT pio.source_file, pio.row_number
FROM points_in_outages pio
WHERE pio.building_id = 1565816
) pio
ON o.source_file = pio.source_file
AND o.row_number = pio.row_number
"
)
outages_1_per_day <- outages %>%
group_by(outage_date) %>%
slice(1) %>%
ungroup()
mutate(frame_id = row_number())
outages_1_per_day$date_label <- paste("Date de la panne:", outages_1_per_day$outage_date)
# Bounding box with padding
bbox <- st_bbox(center_point)
xrange <- bbox["xmax"] - bbox["xmin"]
yrange <- bbox["ymax"] - bbox["ymin"]
bbox_expanded <- c(
left = as.numeric(bbox["xmin"] - 400 * xrange),
bottom = as.numeric(bbox["ymin"] - 400 * xrange),
right = as.numeric(bbox["xmin"] + 400 * xrange),
top = as.numeric(bbox["ymin"] + 400 * xrange)
)
ggmap::register_stadiamaps(key = Sys.getenv("stadiamaps"))
# Fetch free Stamen tile map via Stadia (static map)
stadia_map <- get_stadiamap(
bbox = bbox_expanded,
maptype = "stamen_terrain", # try "stamen_terrain", "stamen_toner", "stamen_watercolor"
zoom = 12 # lower zoom avoids tile errors
)
# Get the map bounds to position text in corner
map_bounds <- attr(stadia_map, "bb")
corner_x <- map_bounds$ur.lon
corner_y <- map_bounds$ur.lat
p <- ggmap(stadia_map) +
geom_sf(data = outages_1_per_day, fill = "orange", color = "orange", alpha = 0.3, size = 1.2, inherit.aes = FALSE) +
geom_sf(data = center_point, aes(geometry = geom_inside), fill = "black", size = 4, alpha = 1, color = "black", inherit.aes = FALSE) +
geom_text(data = outages_1_per_day, aes(label = date_label),
x = corner_x, y = corner_y, hjust = 1, vjust = 1,
size = 4, color = "white", fontface = "bold", inherit.aes = FALSE) +
theme_void() +
theme(
legend.position = "none",
plot.margin = margin(0, 0, 0, 0, "cm")
)
anim <- p + transition_manual(frames = outages_1_per_day$outage_date)
animate(anim, nframes = nrow(outages_1_per_day), fps = 2, width = 600, height = 600,
renderer = gifski_renderer("house_outages.gif"))
library(magick)
img <- image_read("house_outages.gif")
img_trimmed <- image_trim(img)
image_write(img_trimmed, "house_outages_trimmed.gif")
Bon, alors j’imagine que la réponse à ma question c’est la série de blocs à Montréal, avec 88 jours avec une panne de courant durant les 1 184 jours observés (7.4% des jours ). Une moyenne de 2.2 pannes par mois. Mais je ne trouve pas ça naturel, parce que le gros des pannes provient d’une longue coupure à l’été 2023.
Suite à cette disqualification, je décerne donc le championnat au pauvre hère qui habite au 202 Côté Estate du Lac à Lac des Seize avec ses 68 jours avec une panne.
On le salue bien bas. Bravo.
Quelques pistes à regarder avec ce data dans l’avenir:
# # TODO: Quelle est la distribution du % de jour sans panne ?
Je me demande c’est quoi la médiane du nombre de jour de panne par année en 2023 et 2024 (années complètes). mais faut pas oublier que si t’es pas dans mes outages alors t’as eu 0 pannes..
Si j’étais une personne importante dans ma région, j’aimerais savoir quelle est la fréquence où des grosses pannes se produisent. On va regarder la distribution du nombre quotidiens de buildings affectés pour chaque région administrative.
Faut commencer par aller chercher les polygones des régions administratives.
… mais ça c’est pour une autre fois.
Pour ceux que ça intéresse (absolument personne) , voici comment j’ai fait quelques trucs. C’est surtout pour m’aider à ne pas l’oublier.
pour organiser toutes les backups.
Dans le repo hydroquebecoutage, il y a :
Le script hydroquebecoutage.R qui utilise le package rvestpour déterminer la liste des fichiers qui sont disponibles et qui les télécharge.
La github action download_and_save_to_git.yaml qui s’exécute plusieurs fois par jour. Sa job est d’exécuter hydroquebecoutage.R, puis de git committous les fichiers télécharges dans le répertoire data/du repo.
Le répertoire data/ contient aujourd’hui plus de 39 000 fichiers.
Ça ressemble un peu à ça:
st_write(
kmz_data,
con,
"outage",
append = TRUE,
row.names = FALSE,
layer_options = "GEOMETRY_NAME=geom"
)J’ai aussi créé un index géospatial pour que ça aille un peu plus vite quand on va vouloir déterminer quelle bâtiment est dans quelle panne..
dbExecute(con, "
CREATE INDEX IF NOT EXISTS outages_start_idx
ON outages(outage_start DESC);
")Je l’ai téléchargé en format GPKG et je l’ai aussi ajouté à mon serveur PostGIS avec ogr2ogr. Tant qu’à être là, on reprojette en CRS4326 pour être compatible avec la projection des pannes.
ogr2ogr -f "PostgreSQL" PG:"host=192.168.2.15 port=5432 dbname=dbname user=user password=password" "C:\Users\Simon\Downloads\RefBati_GPKG\Referentiel_batiment.gpkg" -nln batiments_4326 -overwrite -t_srs EPSG:4326 -lco GEOMETRY_NAME=geom -lco FID=idOn fait unindex sur lui aussi:
CREATE INDEX idx_batiments_geom ON batiments_4326 USING GIST(geom);
@online{coulombe2025,
author = {Coulombe, Simon},
title = {Qui a manqué le plus de courant au Québec depuis 2022\,?},
date = {2025-08-28},
url = {https://www.simoncoulombe.com/posts/2025-08-28-qui-a-manqué-le-plus-de-courant-au-québec-depuis-2022-/},
langid = {fr}
}