SPB Git

spb/lou-ka Public

Lou·Ka — tous les logements à louer du Québec, un seul endroit.

HTML 99.7%
19.4 KB · 465 lines python
Raw Blame History
1#!/usr/bin/env python32# -*- coding: utf-8 -*-3"""4build_recensement.py — Construit data/staging-recensement.db pour Lou-Ka.56Tables produites (schéma contractuel) :7  1. da_poly(dauid, lat_min, lat_max, lng_min, lng_max, poly)8     - polygones des aires de diffusion (AD) 2021 du Québec, en WGS84,9       restreints à la zone de couverture des annonces10       (lat 44.9–47.7 / lng -74.6–-70.2 : Québec, Lévis, Grand Montréal + périphérie).11  2. da_stats(dauid, population, densite, age_median, revenu_median,12              pct_locataires, loyer_moyen, pct_francais, pct_univ)13     - statistiques du Recensement de 2021 (Profil du recensement).14  3. meta(cle, valeur) — provenance et attribution.1516Sources (Statistique Canada, licence ouverte) :17  - Polygones : Fichiers des limites cartographiques 2021 (92-160-X / index 92-169-X),18    shapefile « aires de diffusion », projection Lambert conforme de19    Statistique Canada (EPSG:3347) :20    https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/boundary-limites/files-fichiers/lda_000b21a_e.zip21    (page index : https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/boundary-limites/index2021-fra.cfm)22  - Statistiques : Profil du recensement 2021, 98-401-X2021006, fichier23    « Québec » au niveau des aires de diffusion (CSV, ~2 Go décompressé) :24    https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/details/download-telecharger/comp/GetFile.cfm?Lang=F&FILETYPE=CSV&GEONO=006_Quebec25    (page : https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/details/download-telecharger.cfm?Lang=F)2627Usage :28  .venv/bin/python scripts/build_recensement.py \29      --shp /tmp/louka-recensement/lda/lda_000b21a_e.shp \30      --csv /tmp/louka-recensement/profil/<fichier>.csv \31      --db  data/staging-recensement.db \32      [--etape poly|stats|valider|tout]3334Le CSV du profil est STREAMÉ (jamais chargé en mémoire). Les valeurs35supprimées par StatCan (confidentialité : symboles x, F, ..., etc.)36restent NULL — rien n'est inventé.37"""3839import argparse40import csv41import datetime42import json43import os44import sqlite345import sys4647import pyproj48import shapefile  # pyshp4950# ---------------------------------------------------------------------------51# Zone de couverture (bbox des annonces Lou-Ka, WGS84)52# 2026-08-08 : élargie à toute la province habitée (expansion provinciale —53# Gatineau, Estrie, Mauricie, Saguenay, Abitibi, Bas-Saint-Laurent, Côte-Nord,54# Gaspésie). Toutes les AD 24* passent le filtre.55# ---------------------------------------------------------------------------56LAT_MIN, LAT_MAX = 44.5, 63.057LNG_MIN, LNG_MAX = -80.0, -56.05859# Projection source du shapefile StatCan : Lambert conforme conique60# « NAD83_Statistics_Canada_Lambert » = EPSG:3347. Cible : WGS84 (EPSG:4326).61TRANSFORMER = pyproj.Transformer.from_crs("EPSG:3347", "EPSG:4326", always_xy=True)6263# ---------------------------------------------------------------------------64# Variables retenues dans le Profil du recensement (98-401-X2021006, français).65# Chaque entrée : nom exact de la caractéristique (colonne NOM_CARACTÉRISTIQUE,66# sans l'indentation) -> clé interne. Les ID (ID_CARACTÉRISTIQUE) sont67# découverts dynamiquement en scannant la liste des caractéristiques de la68# première géographie du fichier, ce qui rend le script robuste aux69# renumérotations éventuelles.70# ---------------------------------------------------------------------------71CARACTERISTIQUES = {72    # population et densité (ID attendus : 1 et 6)73    "Population, 2021": "population",74    "Densité de la population au kilomètre carré": "densite",75    # âge médian de la population (ID 40)76    "Âge médian de la population": "age_median",77    # revenu total médian des ménages en 2020 (ID 243)78    "Revenu total médian des ménages en 2020 ($)": "revenu_median",79    # mode d'occupation : total des ménages (ID 1414) ; « Locataire »80    # (ID 1416) est résolu par position (voir ENFANTS)81    "Total - Ménages privés selon le mode d'occupation - Données-échantillon (25 %)": "menages_total",82    # loyer mensuel moyen des logements loués = frais de logement mensuels83    # moyens des ménages locataires (ID 1495)84    "Frais de logement mensuels moyens pour les logements occupés par un ménage locataire ($)": "loyer_moyen",85    # langue parlée le plus souvent à la maison : total (ID 735) ; le86    # « Français » (ID 739) est résolu par position (voir ENFANTS)87    "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",88    # scolarité 25-64 ans : total (ID 2014) ; « Baccalauréat ou grade89    # supérieur » (ID 2024) est résolu par position (voir ENFANTS)90    "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",91}9293# Sous-caractéristiques ambiguës (leur nom apparaît plusieurs fois dans le94# profil, p. ex. « Français » dans chaque section linguistique) : on prend la95# PREMIÈRE occurrence qui SUIT la ligne « Total » de leur section.96#   clé du total parent -> (nom normalisé de l'enfant, clé interne de l'enfant)97ENFANTS = {98    "menages_total": ("Locataire", "locataires"),               # ID 141699    "langue_total": ("Français", "francais"),                   # ID 739100    "scol_total": ("Baccalauréat ou grade supérieur", "pct_univ_taux"),  # ID 2024101}102103104def _normaliser_nom(nom):105    """Nettoie un NOM_CARACTÉRISTIQUE : enlève l'indentation et uniformise106    les apostrophes (le fichier mélange 0x27 et 0x92/'’')."""107    return nom.strip().replace("\x92", "'").replace("’", "'")108109110# clés normalisées une fois pour toutes111CARACTERISTIQUES = {_normaliser_nom(k): v for k, v in CARACTERISTIQUES.items()}112113114def log(msg):115    print(msg, flush=True)116117118# ===========================================================================119# ÉTAPE 1 — POLYGONES (da_poly)120# ===========================================================================121def construire_poly(chemin_shp, conn):122    """Lit le shapefile des AD, filtre Québec (PRUID 24) + bbox de couverture,123    reprojette EPSG:3347 -> WGS84 et écrit da_poly."""124    conn.execute("DROP TABLE IF EXISTS da_poly")125    conn.execute(126        """CREATE TABLE da_poly(127               dauid   TEXT PRIMARY KEY,128               lat_min REAL, lat_max REAL,129               lng_min REAL, lng_max REAL,130               poly    TEXT)"""131    )132133    lecteur = shapefile.Reader(chemin_shp)134    champs = [f[0] for f in lecteur.fields[1:]]135    i_dauid = champs.index("DAUID")136    i_pruid = champs.index("PRUID")137138    n_qc = n_gardees = 0139    lot = []140    for sr in lecteur.iterShapeRecords():141        rec = sr.record142        if rec[i_pruid] != "24":143            continue144        n_qc += 1145        shp = sr.shape146147        # Pré-filtre grossier : on reprojette les 4 coins de la bbox projetée ;148        # si même la bbox élargie ne touche pas la zone, on saute la géométrie149        # complète (économise ~40 % du temps de reprojection).150        x0, y0, x1, y1 = shp.bbox151        cx, cy = TRANSFORMER.transform([x0, x1, x0, x1], [y0, y0, y1, y1])152        marge = 0.05  # la bbox projetée n'est pas alignée sur les méridiens153        if (max(cy) + marge < LAT_MIN or min(cy) - marge > LAT_MAX154                or max(cx) + marge < LNG_MIN or min(cx) - marge > LNG_MAX):155            continue156157        # Reprojection complète de la géométrie (tous les anneaux).158        pts = shp.points159        xs = [p[0] for p in pts]160        ys = [p[1] for p in pts]161        lngs, lats = TRANSFORMER.transform(xs, ys)162163        lat_min, lat_max = min(lats), max(lats)164        lng_min, lng_max = min(lngs), max(lngs)165        # Filtre exact : intersection de la bbox WGS84 avec la zone visée.166        if (lat_max < LAT_MIN or lat_min > LAT_MAX167                or lng_max < LNG_MIN or lng_min > LNG_MAX):168            continue169170        # Reconstitution des anneaux (parts) : le 1er anneau d'un polygone171        # shapefile est l'anneau extérieur ; les trous / autres polygones172        # suivent. Coordonnées [lng, lat], arrondies à 6 décimales (~10 cm).173        bornes = list(shp.parts) + [len(pts)]174        anneaux = []175        for a in range(len(shp.parts)):176            d, f = bornes[a], bornes[a + 1]177            anneaux.append(178                [[round(lngs[i], 6), round(lats[i], 6)] for i in range(d, f)]179            )180181        lot.append((182            rec[i_dauid],183            round(lat_min, 6), round(lat_max, 6),184            round(lng_min, 6), round(lng_max, 6),185            json.dumps(anneaux, separators=(",", ":")),186        ))187        n_gardees += 1188        if len(lot) >= 500:189            conn.executemany("INSERT INTO da_poly VALUES (?,?,?,?,?,?)", lot)190            lot = []191192    if lot:193        conn.executemany("INSERT INTO da_poly VALUES (?,?,?,?,?,?)", lot)194    conn.commit()195    log(f"da_poly : {n_gardees} AD retenues sur {n_qc} AD au Québec")196    return n_gardees197198199# ===========================================================================200# ÉTAPE 2 — STATISTIQUES (da_stats)201# ===========================================================================202def _detecter_colonnes(entete):203    """Repère les colonnes utiles du CSV (le fichier FR utilise des noms204    français ; on tolère les variantes EN par prudence)."""205    def trouver(*candidats):206        for c in candidats:207            for i, nom in enumerate(entete):208                if nom.strip().lstrip("").upper() == c.upper():209                    return i210        raise KeyError(f"colonne introuvable : {candidats} dans {entete}")211212    return {213        "geo_level": trouver("NIVEAU_GÉO", "GEO_LEVEL"),214        "alt_geo": trouver("CODE_GÉO_ALT", "ALT_GEO_CODE"),215        "car_id": trouver("ID_CARACTÉRISTIQUE", "CHARACTERISTIC_ID"),216        "car_nom": trouver("NOM_CARACTÉRISTIQUE", "CHARACTERISTIC_NAME"),217        "total": trouver("C1_CHIFFRE_TOTAL", "C1_COUNT_TOTAL"),218    }219220221def decouvrir_ids(chemin_csv):222    """1re passe (rapide) : lit les caractéristiques de la première géographie223    du fichier pour associer chaque variable à son ID_CARACTÉRISTIQUE."""224    ids = {}225    en_attente = []  # [(nom enfant recherché, clé interne)]226    with open(chemin_csv, newline="", encoding="latin-1") as f:227        lecteur = csv.reader(f)228        entete = next(lecteur)229        col = _detecter_colonnes(entete)230        premier_geo = None231        for ligne in lecteur:232            geo = ligne[col["alt_geo"]]233            if premier_geo is None:234                premier_geo = geo235            elif geo != premier_geo:236                break  # une géographie = la liste complète des caractéristiques237            cid = int(ligne[col["car_id"]])238            nom = _normaliser_nom(ligne[col["car_nom"]])239            # sous-caractéristique attendue après son « Total » parent ?240            for i, (nom_enfant, cle_enfant) in enumerate(en_attente):241                if nom == nom_enfant and cle_enfant not in ids:242                    ids[cle_enfant] = cid243                    en_attente.pop(i)244                    break245            cle = CARACTERISTIQUES.get(nom)246            if cle and cle not in ids:247                ids[cle] = cid248                if cle in ENFANTS:249                    en_attente.append(ENFANTS[cle])250    return ids, col251252253def _nombre(txt):254    """Convertit une cellule du profil en float, ou None si valeur supprimée255    (x, F, .., ..., vide) — on n'invente rien."""256    txt = (txt or "").strip().replace(",", ".")257    if not txt or txt in {"x", "F", "..", "...", "r", "t"}:258        return None259    try:260        return float(txt)261    except ValueError:262        return None263264265def construire_stats(chemin_csv, conn):266    """2e passe : streame tout le CSV, ne garde que les lignes des AD retenues267    dans da_poly et les 12 caractéristiques utiles, puis calcule les 8268    variables finales."""269    dauids = {r[0] for r in conn.execute("SELECT dauid FROM da_poly")}270    if not dauids:271        raise SystemExit("da_poly est vide — lancer l'étape poly d'abord")272273    ids, col = decouvrir_ids(chemin_csv)274    attendues = set(CARACTERISTIQUES.values()) | {c for _, c in ENFANTS.values()}275    manquants = attendues - set(ids)276    if manquants:277        raise SystemExit(f"IDs de caractéristiques introuvables : {manquants}")278    log(f"IDs des caractéristiques : {ids}")279    ids_voulus = {v: k for k, v in ids.items()}  # cid -> clé interne280281    # accumulation : dauid -> {clé interne: valeur brute}282    donnees = {}283    niveaux_ad = {"Aire de diffusion", "Dissemination area"}284    n_lignes = 0285    with open(chemin_csv, newline="", encoding="latin-1") as f:286        lecteur = csv.reader(f)287        next(lecteur)  # entête288        i_niv, i_geo = col["geo_level"], col["alt_geo"]289        i_cid, i_tot = col["car_id"], col["total"]290        for ligne in lecteur:291            n_lignes += 1292            if ligne[i_niv] not in niveaux_ad:293                continue294            geo = ligne[i_geo]295            if geo not in dauids:296                continue297            cle = ids_voulus.get(int(ligne[i_cid]))298            if cle is None:299                continue300            donnees.setdefault(geo, {})[cle] = _nombre(ligne[i_tot])301    log(f"CSV : {n_lignes} lignes lues, {len(donnees)} AD avec données")302303    conn.execute("DROP TABLE IF EXISTS da_stats")304    conn.execute(305        """CREATE TABLE da_stats(306               dauid          TEXT PRIMARY KEY,307               population     INTEGER,308               densite        REAL,309               age_median     REAL,310               revenu_median  REAL,311               pct_locataires REAL,312               loyer_moyen    REAL,313               pct_francais   REAL,314               pct_univ       REAL)"""315    )316317    def pct(part, total):318        """Pourcentage part/total, NULL si l'un des deux est supprimé/nul."""319        if part is None or not total:320            return None321        return round(100.0 * part / total, 1)322323    lot = []324    for dauid, d in donnees.items():325        lot.append((326            dauid,327            int(d["population"]) if d.get("population") is not None else None,328            d.get("densite"),329            d.get("age_median"),330            d.get("revenu_median"),331            pct(d.get("locataires"), d.get("menages_total")),332            d.get("loyer_moyen"),333            pct(d.get("francais"), d.get("langue_total")),334            pct(d.get("pct_univ_taux"), d.get("scol_total")),335        ))336    conn.executemany("INSERT INTO da_stats VALUES (?,?,?,?,?,?,?,?,?)", lot)337    conn.commit()338    log(f"da_stats : {len(lot)} AD insérées")339    return len(lot)340341342# ===========================================================================343# ÉTAPE 3 — MÉTADONNÉES (meta)344# ===========================================================================345def ecrire_meta(conn):346    conn.execute("DROP TABLE IF EXISTS meta")347    conn.execute("CREATE TABLE meta(cle TEXT PRIMARY KEY, valeur TEXT)")348    meta = {349        "source_polygones": (350            "Statistique Canada, Fichiers des limites cartographiques du "351            "Recensement de 2021 (92-160-X), aires de diffusion, "352            "lda_000b21a_e.zip — https://www12.statcan.gc.ca/census-recensement/"353            "2021/geo/sip-pis/boundary-limites/files-fichiers/lda_000b21a_e.zip"354        ),355        "source_stats": (356            "Statistique Canada, Profil du recensement, Recensement de la "357            "population de 2021, no 98-401-X2021006 au catalogue, fichier "358            "Québec au niveau des aires de diffusion — "359            "https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/"360            "details/download-telecharger.cfm?Lang=F"361        ),362        "date_construction": datetime.date.today().isoformat(),363        "attribution": (364            "Statistique Canada, Recensement de la population de 2021 "365            "(reproduit et diffusé « tel quel » avec la permission de "366            "Statistique Canada — Licence ouverte de Statistique Canada)"367        ),368        "projection_source": "EPSG:3347 (NAD83 Statistics Canada Lambert) -> EPSG:4326 (WGS84)",369        "zone_couverture": f"lat {LAT_MIN}{LAT_MAX}, lng {LNG_MIN}{LNG_MAX} (Québec/Lévis/Grand Montréal + périphérie)",370    }371    conn.executemany("INSERT INTO meta VALUES (?,?)", meta.items())372    conn.commit()373374375# ===========================================================================376# VALIDATIONS377# ===========================================================================378def point_dans_poly(lat, lng, anneaux):379    """Ray casting pair/impair sur l'ensemble des anneaux (gère les trous)."""380    dedans = False381    for anneau in anneaux:382        n = len(anneau)383        j = n - 1384        for i in range(n):385            xi, yi = anneau[i]386            xj, yj = anneau[j]387            if (yi > lat) != (yj > lat) and \388               lng < (xj - xi) * (lat - yi) / (yj - yi) + xi:389                dedans = not dedans390            j = i391    return dedans392393394def trouver_dauid(conn, lat, lng):395    """Retourne le DAUID contenant le point (candidats par bbox, puis PIP)."""396    for dauid, poly in conn.execute(397        "SELECT dauid, poly FROM da_poly "398        "WHERE ? BETWEEN lat_min AND lat_max AND ? BETWEEN lng_min AND lng_max",399        (lat, lng),400    ):401        if point_dans_poly(lat, lng, json.loads(poly)):402            return dauid403    return None404405406def valider(conn):407    n_poly = conn.execute("SELECT COUNT(*) FROM da_poly").fetchone()[0]408    n_stats = conn.execute(409        "SELECT COUNT(*) FROM da_stats WHERE dauid IN (SELECT dauid FROM da_poly)"410    ).fetchone()[0]411    log(f"Validation : {n_poly} AD dans da_poly ; couverture da_stats = "412        f"{n_stats}/{n_poly} ({100.0 * n_stats / max(n_poly, 1):.1f} %)")413414    # point de contrôle : colline Parlementaire / Vieux-Québec415    dauid = trouver_dauid(conn, 46.8139, -71.2329)416    ok = dauid == "24231035"417    log(f"Point (46.8139, -71.2329) -> DAUID {dauid} "418        f"({'OK' if ok else 'ÉCHEC, attendu 24231035'})")419420    # échantillon lisible421    for nom, lat, lng in [422        ("Saint-Roch (Québec)", 46.8163, -71.2258),423        ("Sillery (Québec)", 46.7702, -71.2601),424        ("Plateau Mont-Royal (Mtl)", 45.5230, -73.5817),425    ]:426        d = trouver_dauid(conn, lat, lng)427        row = conn.execute(428            "SELECT population, revenu_median, pct_locataires, loyer_moyen, "429            "pct_francais, pct_univ FROM da_stats WHERE dauid=?", (d,)430        ).fetchone() if d else None431        log(f"  {nom}: DAUID={d} stats={row}")432    return ok and n_poly > 0 and n_stats >= 0.95 * n_poly433434435# ===========================================================================436def main():437    ap = argparse.ArgumentParser(description=__doc__)438    ap.add_argument("--shp", help="chemin du shapefile lda_000b21a_e.shp")439    ap.add_argument("--csv", help="chemin du CSV du profil (98-401-X2021006, Québec)")440    ap.add_argument("--db", default="data/staging-recensement.db")441    ap.add_argument("--etape", default="tout",442                    choices=["poly", "stats", "valider", "tout"])443    args = ap.parse_args()444445    conn = sqlite3.connect(args.db)446    try:447        if args.etape in ("poly", "tout"):448            if not args.shp:449                ap.error("--shp requis pour l'étape poly")450            construire_poly(args.shp, conn)451        if args.etape in ("stats", "tout"):452            if not args.csv:453                ap.error("--csv requis pour l'étape stats")454            construire_stats(args.csv, conn)455            ecrire_meta(conn)456        if args.etape in ("valider", "tout"):457            ok = valider(conn)458            sys.exit(0 if ok else 1)459    finally:460        conn.close()461462463if __name__ == "__main__":464    main()465