#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ build_recensement.py — Construit data/staging-recensement.db pour Lou-Ka. Tables produites (schéma contractuel) : 1. da_poly(dauid, lat_min, lat_max, lng_min, lng_max, poly) - polygones des aires de diffusion (AD) 2021 du Québec, en WGS84, restreints à la zone de couverture des annonces (lat 44.9–47.7 / lng -74.6–-70.2 : Québec, Lévis, Grand Montréal + périphérie). 2. da_stats(dauid, population, densite, age_median, revenu_median, pct_locataires, loyer_moyen, pct_francais, pct_univ) - statistiques du Recensement de 2021 (Profil du recensement). 3. meta(cle, valeur) — provenance et attribution. Sources (Statistique Canada, licence ouverte) : - Polygones : Fichiers des limites cartographiques 2021 (92-160-X / index 92-169-X), shapefile « aires de diffusion », projection Lambert conforme de Statistique Canada (EPSG:3347) : https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/boundary-limites/files-fichiers/lda_000b21a_e.zip (page index : https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/boundary-limites/index2021-fra.cfm) - Statistiques : Profil du recensement 2021, 98-401-X2021006, fichier « Québec » au niveau des aires de diffusion (CSV, ~2 Go décompressé) : https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/details/download-telecharger/comp/GetFile.cfm?Lang=F&FILETYPE=CSV&GEONO=006_Quebec (page : https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/details/download-telecharger.cfm?Lang=F) Usage : .venv/bin/python scripts/build_recensement.py \ --shp /tmp/louka-recensement/lda/lda_000b21a_e.shp \ --csv /tmp/louka-recensement/profil/.csv \ --db data/staging-recensement.db \ [--etape poly|stats|valider|tout] Le CSV du profil est STREAMÉ (jamais chargé en mémoire). Les valeurs supprimées par StatCan (confidentialité : symboles x, F, ..., etc.) restent NULL — rien n'est inventé. """ import argparse import csv import datetime import json import os import sqlite3 import sys import pyproj import shapefile # pyshp # --------------------------------------------------------------------------- # Zone de couverture (bbox des annonces Lou-Ka, WGS84) # 2026-08-08 : élargie à toute la province habitée (expansion provinciale — # Gatineau, Estrie, Mauricie, Saguenay, Abitibi, Bas-Saint-Laurent, Côte-Nord, # Gaspésie). Toutes les AD 24* passent le filtre. # --------------------------------------------------------------------------- LAT_MIN, LAT_MAX = 44.5, 63.0 LNG_MIN, LNG_MAX = -80.0, -56.0 # Projection source du shapefile StatCan : Lambert conforme conique # « NAD83_Statistics_Canada_Lambert » = EPSG:3347. Cible : WGS84 (EPSG:4326). TRANSFORMER = pyproj.Transformer.from_crs("EPSG:3347", "EPSG:4326", always_xy=True) # --------------------------------------------------------------------------- # Variables retenues dans le Profil du recensement (98-401-X2021006, français). # Chaque entrée : nom exact de la caractéristique (colonne NOM_CARACTÉRISTIQUE, # sans l'indentation) -> clé interne. Les ID (ID_CARACTÉRISTIQUE) sont # découverts dynamiquement en scannant la liste des caractéristiques de la # première géographie du fichier, ce qui rend le script robuste aux # renumérotations éventuelles. # --------------------------------------------------------------------------- CARACTERISTIQUES = { # population et densité (ID attendus : 1 et 6) "Population, 2021": "population", "Densité de la population au kilomètre carré": "densite", # âge médian de la population (ID 40) "Âge médian de la population": "age_median", # revenu total médian des ménages en 2020 (ID 243) "Revenu total médian des ménages en 2020 ($)": "revenu_median", # mode d'occupation : total des ménages (ID 1414) ; « Locataire » # (ID 1416) est résolu par position (voir ENFANTS) "Total - Ménages privés selon le mode d'occupation - Données-échantillon (25 %)": "menages_total", # loyer mensuel moyen des logements loués = frais de logement mensuels # moyens des ménages locataires (ID 1495) "Frais de logement mensuels moyens pour les logements occupés par un ménage locataire ($)": "loyer_moyen", # langue parlée le plus souvent à la maison : total (ID 735) ; le # « Français » (ID 739) est résolu par position (voir ENFANTS) "Total - Langue parlée le plus souvent à la maison pour la population totale à l'exclusion des résidents d'un établissement institutionnel - Données intégrales (100 %)": "langue_total", # scolarité 25-64 ans : total (ID 2014) ; « Baccalauréat ou grade # supérieur » (ID 2024) est résolu par position (voir ENFANTS) "Total - Plus haut certificat, diplôme ou grade pour la population âgée de 25 à 64 ans dans les ménages privés - Données-échantillon (25 %)": "scol_total", } # Sous-caractéristiques ambiguës (leur nom apparaît plusieurs fois dans le # profil, p. ex. « Français » dans chaque section linguistique) : on prend la # PREMIÈRE occurrence qui SUIT la ligne « Total » de leur section. # clé du total parent -> (nom normalisé de l'enfant, clé interne de l'enfant) ENFANTS = { "menages_total": ("Locataire", "locataires"), # ID 1416 "langue_total": ("Français", "francais"), # ID 739 "scol_total": ("Baccalauréat ou grade supérieur", "pct_univ_taux"), # ID 2024 } def _normaliser_nom(nom): """Nettoie un NOM_CARACTÉRISTIQUE : enlève l'indentation et uniformise les apostrophes (le fichier mélange 0x27 et 0x92/'’').""" return nom.strip().replace("\x92", "'").replace("’", "'") # clés normalisées une fois pour toutes CARACTERISTIQUES = {_normaliser_nom(k): v for k, v in CARACTERISTIQUES.items()} def log(msg): print(msg, flush=True) # =========================================================================== # ÉTAPE 1 — POLYGONES (da_poly) # =========================================================================== def construire_poly(chemin_shp, conn): """Lit le shapefile des AD, filtre Québec (PRUID 24) + bbox de couverture, reprojette EPSG:3347 -> WGS84 et écrit da_poly.""" conn.execute("DROP TABLE IF EXISTS da_poly") conn.execute( """CREATE TABLE da_poly( dauid TEXT PRIMARY KEY, lat_min REAL, lat_max REAL, lng_min REAL, lng_max REAL, poly TEXT)""" ) lecteur = shapefile.Reader(chemin_shp) champs = [f[0] for f in lecteur.fields[1:]] i_dauid = champs.index("DAUID") i_pruid = champs.index("PRUID") n_qc = n_gardees = 0 lot = [] for sr in lecteur.iterShapeRecords(): rec = sr.record if rec[i_pruid] != "24": continue n_qc += 1 shp = sr.shape # Pré-filtre grossier : on reprojette les 4 coins de la bbox projetée ; # si même la bbox élargie ne touche pas la zone, on saute la géométrie # complète (économise ~40 % du temps de reprojection). x0, y0, x1, y1 = shp.bbox cx, cy = TRANSFORMER.transform([x0, x1, x0, x1], [y0, y0, y1, y1]) marge = 0.05 # la bbox projetée n'est pas alignée sur les méridiens if (max(cy) + marge < LAT_MIN or min(cy) - marge > LAT_MAX or max(cx) + marge < LNG_MIN or min(cx) - marge > LNG_MAX): continue # Reprojection complète de la géométrie (tous les anneaux). pts = shp.points xs = [p[0] for p in pts] ys = [p[1] for p in pts] lngs, lats = TRANSFORMER.transform(xs, ys) lat_min, lat_max = min(lats), max(lats) lng_min, lng_max = min(lngs), max(lngs) # Filtre exact : intersection de la bbox WGS84 avec la zone visée. if (lat_max < LAT_MIN or lat_min > LAT_MAX or lng_max < LNG_MIN or lng_min > LNG_MAX): continue # Reconstitution des anneaux (parts) : le 1er anneau d'un polygone # shapefile est l'anneau extérieur ; les trous / autres polygones # suivent. Coordonnées [lng, lat], arrondies à 6 décimales (~10 cm). bornes = list(shp.parts) + [len(pts)] anneaux = [] for a in range(len(shp.parts)): d, f = bornes[a], bornes[a + 1] anneaux.append( [[round(lngs[i], 6), round(lats[i], 6)] for i in range(d, f)] ) lot.append(( rec[i_dauid], round(lat_min, 6), round(lat_max, 6), round(lng_min, 6), round(lng_max, 6), json.dumps(anneaux, separators=(",", ":")), )) n_gardees += 1 if len(lot) >= 500: conn.executemany("INSERT INTO da_poly VALUES (?,?,?,?,?,?)", lot) lot = [] if lot: conn.executemany("INSERT INTO da_poly VALUES (?,?,?,?,?,?)", lot) conn.commit() log(f"da_poly : {n_gardees} AD retenues sur {n_qc} AD au Québec") return n_gardees # =========================================================================== # ÉTAPE 2 — STATISTIQUES (da_stats) # =========================================================================== def _detecter_colonnes(entete): """Repère les colonnes utiles du CSV (le fichier FR utilise des noms français ; on tolère les variantes EN par prudence).""" def trouver(*candidats): for c in candidats: for i, nom in enumerate(entete): if nom.strip().lstrip("").upper() == c.upper(): return i raise KeyError(f"colonne introuvable : {candidats} dans {entete}") return { "geo_level": trouver("NIVEAU_GÉO", "GEO_LEVEL"), "alt_geo": trouver("CODE_GÉO_ALT", "ALT_GEO_CODE"), "car_id": trouver("ID_CARACTÉRISTIQUE", "CHARACTERISTIC_ID"), "car_nom": trouver("NOM_CARACTÉRISTIQUE", "CHARACTERISTIC_NAME"), "total": trouver("C1_CHIFFRE_TOTAL", "C1_COUNT_TOTAL"), } def decouvrir_ids(chemin_csv): """1re passe (rapide) : lit les caractéristiques de la première géographie du fichier pour associer chaque variable à son ID_CARACTÉRISTIQUE.""" ids = {} en_attente = [] # [(nom enfant recherché, clé interne)] with open(chemin_csv, newline="", encoding="latin-1") as f: lecteur = csv.reader(f) entete = next(lecteur) col = _detecter_colonnes(entete) premier_geo = None for ligne in lecteur: geo = ligne[col["alt_geo"]] if premier_geo is None: premier_geo = geo elif geo != premier_geo: break # une géographie = la liste complète des caractéristiques cid = int(ligne[col["car_id"]]) nom = _normaliser_nom(ligne[col["car_nom"]]) # sous-caractéristique attendue après son « Total » parent ? for i, (nom_enfant, cle_enfant) in enumerate(en_attente): if nom == nom_enfant and cle_enfant not in ids: ids[cle_enfant] = cid en_attente.pop(i) break cle = CARACTERISTIQUES.get(nom) if cle and cle not in ids: ids[cle] = cid if cle in ENFANTS: en_attente.append(ENFANTS[cle]) return ids, col def _nombre(txt): """Convertit une cellule du profil en float, ou None si valeur supprimée (x, F, .., ..., vide) — on n'invente rien.""" txt = (txt or "").strip().replace(",", ".") if not txt or txt in {"x", "F", "..", "...", "r", "t"}: return None try: return float(txt) except ValueError: return None def construire_stats(chemin_csv, conn): """2e passe : streame tout le CSV, ne garde que les lignes des AD retenues dans da_poly et les 12 caractéristiques utiles, puis calcule les 8 variables finales.""" dauids = {r[0] for r in conn.execute("SELECT dauid FROM da_poly")} if not dauids: raise SystemExit("da_poly est vide — lancer l'étape poly d'abord") ids, col = decouvrir_ids(chemin_csv) attendues = set(CARACTERISTIQUES.values()) | {c for _, c in ENFANTS.values()} manquants = attendues - set(ids) if manquants: raise SystemExit(f"IDs de caractéristiques introuvables : {manquants}") log(f"IDs des caractéristiques : {ids}") ids_voulus = {v: k for k, v in ids.items()} # cid -> clé interne # accumulation : dauid -> {clé interne: valeur brute} donnees = {} niveaux_ad = {"Aire de diffusion", "Dissemination area"} n_lignes = 0 with open(chemin_csv, newline="", encoding="latin-1") as f: lecteur = csv.reader(f) next(lecteur) # entête i_niv, i_geo = col["geo_level"], col["alt_geo"] i_cid, i_tot = col["car_id"], col["total"] for ligne in lecteur: n_lignes += 1 if ligne[i_niv] not in niveaux_ad: continue geo = ligne[i_geo] if geo not in dauids: continue cle = ids_voulus.get(int(ligne[i_cid])) if cle is None: continue donnees.setdefault(geo, {})[cle] = _nombre(ligne[i_tot]) log(f"CSV : {n_lignes} lignes lues, {len(donnees)} AD avec données") conn.execute("DROP TABLE IF EXISTS da_stats") conn.execute( """CREATE TABLE da_stats( dauid TEXT PRIMARY KEY, population INTEGER, densite REAL, age_median REAL, revenu_median REAL, pct_locataires REAL, loyer_moyen REAL, pct_francais REAL, pct_univ REAL)""" ) def pct(part, total): """Pourcentage part/total, NULL si l'un des deux est supprimé/nul.""" if part is None or not total: return None return round(100.0 * part / total, 1) lot = [] for dauid, d in donnees.items(): lot.append(( dauid, int(d["population"]) if d.get("population") is not None else None, d.get("densite"), d.get("age_median"), d.get("revenu_median"), pct(d.get("locataires"), d.get("menages_total")), d.get("loyer_moyen"), pct(d.get("francais"), d.get("langue_total")), pct(d.get("pct_univ_taux"), d.get("scol_total")), )) conn.executemany("INSERT INTO da_stats VALUES (?,?,?,?,?,?,?,?,?)", lot) conn.commit() log(f"da_stats : {len(lot)} AD insérées") return len(lot) # =========================================================================== # ÉTAPE 3 — MÉTADONNÉES (meta) # =========================================================================== def ecrire_meta(conn): conn.execute("DROP TABLE IF EXISTS meta") conn.execute("CREATE TABLE meta(cle TEXT PRIMARY KEY, valeur TEXT)") meta = { "source_polygones": ( "Statistique Canada, Fichiers des limites cartographiques du " "Recensement de 2021 (92-160-X), aires de diffusion, " "lda_000b21a_e.zip — https://www12.statcan.gc.ca/census-recensement/" "2021/geo/sip-pis/boundary-limites/files-fichiers/lda_000b21a_e.zip" ), "source_stats": ( "Statistique Canada, Profil du recensement, Recensement de la " "population de 2021, no 98-401-X2021006 au catalogue, fichier " "Québec au niveau des aires de diffusion — " "https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/" "details/download-telecharger.cfm?Lang=F" ), "date_construction": datetime.date.today().isoformat(), "attribution": ( "Statistique Canada, Recensement de la population de 2021 " "(reproduit et diffusé « tel quel » avec la permission de " "Statistique Canada — Licence ouverte de Statistique Canada)" ), "projection_source": "EPSG:3347 (NAD83 Statistics Canada Lambert) -> EPSG:4326 (WGS84)", "zone_couverture": f"lat {LAT_MIN}–{LAT_MAX}, lng {LNG_MIN}–{LNG_MAX} (Québec/Lévis/Grand Montréal + périphérie)", } conn.executemany("INSERT INTO meta VALUES (?,?)", meta.items()) conn.commit() # =========================================================================== # VALIDATIONS # =========================================================================== def point_dans_poly(lat, lng, anneaux): """Ray casting pair/impair sur l'ensemble des anneaux (gère les trous).""" dedans = False for anneau in anneaux: n = len(anneau) j = n - 1 for i in range(n): xi, yi = anneau[i] xj, yj = anneau[j] if (yi > lat) != (yj > lat) and \ lng < (xj - xi) * (lat - yi) / (yj - yi) + xi: dedans = not dedans j = i return dedans def trouver_dauid(conn, lat, lng): """Retourne le DAUID contenant le point (candidats par bbox, puis PIP).""" for dauid, poly in conn.execute( "SELECT dauid, poly FROM da_poly " "WHERE ? BETWEEN lat_min AND lat_max AND ? BETWEEN lng_min AND lng_max", (lat, lng), ): if point_dans_poly(lat, lng, json.loads(poly)): return dauid return None def valider(conn): n_poly = conn.execute("SELECT COUNT(*) FROM da_poly").fetchone()[0] n_stats = conn.execute( "SELECT COUNT(*) FROM da_stats WHERE dauid IN (SELECT dauid FROM da_poly)" ).fetchone()[0] log(f"Validation : {n_poly} AD dans da_poly ; couverture da_stats = " f"{n_stats}/{n_poly} ({100.0 * n_stats / max(n_poly, 1):.1f} %)") # point de contrôle : colline Parlementaire / Vieux-Québec dauid = trouver_dauid(conn, 46.8139, -71.2329) ok = dauid == "24231035" log(f"Point (46.8139, -71.2329) -> DAUID {dauid} " f"({'OK' if ok else 'ÉCHEC, attendu 24231035'})") # échantillon lisible for nom, lat, lng in [ ("Saint-Roch (Québec)", 46.8163, -71.2258), ("Sillery (Québec)", 46.7702, -71.2601), ("Plateau Mont-Royal (Mtl)", 45.5230, -73.5817), ]: d = trouver_dauid(conn, lat, lng) row = conn.execute( "SELECT population, revenu_median, pct_locataires, loyer_moyen, " "pct_francais, pct_univ FROM da_stats WHERE dauid=?", (d,) ).fetchone() if d else None log(f" {nom}: DAUID={d} stats={row}") return ok and n_poly > 0 and n_stats >= 0.95 * n_poly # =========================================================================== def main(): ap = argparse.ArgumentParser(description=__doc__) ap.add_argument("--shp", help="chemin du shapefile lda_000b21a_e.shp") ap.add_argument("--csv", help="chemin du CSV du profil (98-401-X2021006, Québec)") ap.add_argument("--db", default="data/staging-recensement.db") ap.add_argument("--etape", default="tout", choices=["poly", "stats", "valider", "tout"]) args = ap.parse_args() conn = sqlite3.connect(args.db) try: if args.etape in ("poly", "tout"): if not args.shp: ap.error("--shp requis pour l'étape poly") construire_poly(args.shp, conn) if args.etape in ("stats", "tout"): if not args.csv: ap.error("--csv requis pour l'étape stats") construire_stats(args.csv, conn) ecrire_meta(conn) if args.etape in ("valider", "tout"): ok = valider(conn) sys.exit(0 if ok else 1) finally: conn.close() if __name__ == "__main__": main()