#!/usr/bin/env python """ WP38 - External regional covariates for the 17 Quebec administrative regions. Re-downloads (idempotently, with SHA-256 checksums) every authoritative public source listed in SOURCES, derives the geography products, and rebuilds data/external/regions_external.csv (one row per region, tidy) data/external/manifest.json (one entry per file: url, sha256, vintage, licence, notes) data/external/geo/*.gpkg (MERN admin regions, StatCan 2021 ER / FSA for Quebec) data/external/fsa_to_region_2021.csv (FSA -> region assignment, centroid + dominant-area) data/external/cd_to_region_concordance_2021.csv data/external/ra_er_overlap_2021.csv (does StatCan ER == Quebec RA ? measured, not assumed) data/external/access_log.json (HTTP status of blocked / alternative endpoints) Rules (CLAUDE.md rule 1): nothing is typed in by hand. Every number in the CSV is parsed from a downloaded file whose checksum is in the manifest. Cells that could not be obtained stay EMPTY (never 0). Usage ----- .venv/bin/python src/ingest/external.py # download what is missing, rebuild everything .venv/bin/python src/ingest/external.py --refresh # force re-download of every source .venv/bin/python src/ingest/external.py --no-geo # skip the (slow) StatCan REST geometry pulls .venv/bin/python src/ingest/external.py --strict # fail if a re-downloaded file's sha256 changed Only official hosts are contacted: statistique.quebec.ca (ISQ), www150.statcan.gc.ca (StatCan WDS), geo.statcan.gc.ca (StatCan boundary REST service), www12.statcan.gc.ca (probe only - returns 403 to non-browser clients), diffusion.mern.gouv.qc.ca (MRNF, via donneesquebec.ca catalogue), transformation-numerique.ulaval.ca (ATN / NETendances). """ from __future__ import annotations import argparse import csv import hashlib import io import json import re import sys import time import unicodedata import zipfile from datetime import datetime, timezone from pathlib import Path import pandas as pd import requests PROJECT = Path(__file__).resolve().parents[2] EXT = PROJECT / "data" / "external" DL = EXT / "dl" GEO = EXT / "geo" MANIFEST = EXT / "manifest.json" ACCESS_LOG = EXT / "access_log.json" OUT_CSV = EXT / "regions_external.csv" UA = ("Mozilla/5.0 (Macintosh; Intel Mac OS X 14_0) AppleWebKit/537.36 " "(KHTML, like Gecko) Chrome/124.0.0.0 Safari/537.36 WP38-UQO-research/1.0") SESSION = requests.Session() SESSION.headers.update({"User-Agent": UA}) # -------------------------------------------------------------------------------------- # The 17 administrative regions (official codes and spellings; en-dash in 02 and 11) # -------------------------------------------------------------------------------------- REGIONS = [ ("01", "Bas-Saint-Laurent"), ("02", "Saguenay–Lac-Saint-Jean"), ("03", "Capitale-Nationale"), ("04", "Mauricie"), ("05", "Estrie"), ("06", "Montréal"), ("07", "Outaouais"), ("08", "Abitibi-Témiscamingue"), ("09", "Côte-Nord"), ("10", "Nord-du-Québec"), ("11", "Gaspésie–Îles-de-la-Madeleine"), ("12", "Chaudière-Appalaches"), ("13", "Laval"), ("14", "Lanaudière"), ("15", "Laurentides"), ("16", "Montérégie"), ("17", "Centre-du-Québec"), ] CODE_BY_NAME = {} def norm(s: str) -> str: """Accent/dash/space-insensitive key used to match region names across sources.""" s = str(s) s = s.replace("\n", " ") s = re.sub(r"[\u2013\u2014\-]+", "-", s) # en/em dash and double hyphen -> '-' (before ASCII fold) s = unicodedata.normalize("NFKD", s).encode("ascii", "ignore").decode() s = re.sub(r"\s+", " ", s).strip().lower() s = re.sub(r"\s*\(\d+\)\s*$", "", s) # trailing "(01)" s = re.sub(r"[\d,]+$", "", s) # ISQ footnote digits "Estrie2,3" s = s.replace(", quebec", "") s = s.replace(" ", "") return s for _c, _n in REGIONS: CODE_BY_NAME[norm(_n)] = _c # StatCan spellings (double hyphen) and ISQ footnoted variants map through norm() already. # -------------------------------------------------------------------------------------- # Sources # -------------------------------------------------------------------------------------- LIC_ISQ = ("ISQ - Licence d'utilisation des produits de l'Institut (reproduction permitted with " "attribution; see https://statistique.quebec.ca/fr/institut/politique-de-confidentialite " "and the 'Conditions d'utilisation' page)") LIC_STATCAN = "Statistics Canada Open Licence (https://www.statcan.gc.ca/en/reference/licence)" LIC_MERN = "Creative Commons Attribution 4.0 (CC-BY 4.0), as stated on donneesquebec.ca dataset page" LIC_ATN = "Copyright Académie de la transformation numérique (Université Laval); report PDF kept for citation only" SOURCES = { # ---------------- ISQ ---------------- "isq_pop": dict( file="isq_population_ra_1986_2025_fichier_complet_3595.xlsx", url="https://statistique.quebec.ca/docs-ken/multimedia/Fichier_complet_3595.xlsx", source_name="ISQ - Estimations de la population et indicateurs d'accroissement des régions administratives, " "Québec, 1er juillet 1986 à 2025 (tableau no 3595, 'Fichier complet')", page="https://statistique.quebec.ca/fr/produit/tableau/estimations-population-regions-administratives", licence=LIC_ISQ, vintage="1 July 2025 (provisional, 'p'); 2022-2024 revised; ISQ update 14 Jan 2026; RA boundaries as of 1 July 2025", notes="Source of the ISQ table: Statistics Canada, Annual Demographic Estimates (subprovincial areas, January 2026), " "adapted by ISQ. Estrie/Montérégie reflect the 28 July 2021 boundary change (MRC La Haute-Yamaska and " "Brome-Missisquoi moved from Montérégie to Estrie); Estrie/Chaudière-Appalaches reflect the 1 Jan 2024 " "Courcelles merger.", ), "isq_biz_ind": dict( file="isq_entreprises_actives_ra_industrie_dec2025.xlsx", url="https://statistique.quebec.ca/fr/fichier/entreprises-actives-regions-adminsitratives-industrie-decembre-2025-donnees-semstrielles.xlsx", source_name="ISQ - Nombre d'entreprises actives (avec employés) dans les régions administratives du Québec selon " "l'industrie (SCIAN 2-4 chiffres), décembre 2025 (données semestrielles)", page="https://statistique.quebec.ca/fr/document/nombre-entreprises-actives-quebec", licence=LIC_ISQ, vintage="December 2025 (Business Register snapshot); ISQ compilation of StatCan table 33-10-1097-01, data captured 2 March 2026", notes="Unit = business 'statistical locations' WITH EMPLOYEES only (employer locations). Locations without " "employees are published by ISQ at the province level only. ISQ aggregates StatCan census-subdivision " "counts (33-10-1097-01) to administrative regions; 703 employer locations are 'non classifiées par " "municipalité' and are excluded from the regional totals.", ), "isq_biz_size": dict( file="isq_entreprises_actives_ra_municipalites_dec2025.xlsx", url="https://statistique.quebec.ca/fr/fichier/entreprises-actives-regions-administratives-municipalites-decembre-2025-donnees-semestrielles.xlsx", source_name="ISQ - Nombre d'entreprises avec employés par région administrative et par tranche d'effectif, " "ensemble des industries, décembre 2025", page="https://statistique.quebec.ca/fr/document/nombre-entreprises-actives-quebec", licence=LIC_ISQ, vintage="December 2025", notes="Same universe as isq_biz_ind (employer locations); used as cross-check of regional totals and for the " "employment-size distribution.", ), "isq_income": dict( file="isq_revenu_disponible_ra_tableau_3962.csv", url="https://statistique.quebec.ca/pls/ken/ken411_data_explt_v2.p_retrn_data", params={"p_id_tabl": "3962", "p_champs": "co_tertr,tertr,de_codfc_refrn,mesr,titre_23,perd_23,titre_22,perd_22,titre_21,perd_21," "titre_20,perd_20,titre_19,perd_19,titre_18,perd_18,titre_17,perd_17,titre_16,perd_16," "perd_17_sign,titre_15,perd_15,titre_14,perd_14,titre_13,perd_13,titre_12,perd_12,titre_11," "perd_11,titre_10,perd_10,titre_9,perd_9,titre_8,perd_8,titre_7,perd_7,titre_6,perd_6," "titre_5,perd_5,titre_4,perd_4,titre_3,perd_3,titre_2,perd_2,titre_1,perd_1", "p_tri": "tri_tertr asc,tri_codfc asc"}, source_name="ISQ - Revenu disponible et ses composantes, régions administratives et ensemble du Québec, " "2002-2024 (tableau dynamique no 3962; data endpoint used by the ISQ web table)", page="https://statistique.quebec.ca/fr/document/revenu-disponible/tableau/revenu-disponible-et-composantes-regions-administratives-et-ensemble-du-quebec", licence=LIC_ISQ, vintage="2024 provisional ('p'); 2022-2023 revised; ISQ release 1 June 2026", notes="Disposable income = household-sector aggregate of the System of National Accounts (ISQ regional " "economic accounts), per capita in current dollars ('$/hab'). Not comparable with MFQ tax statistics.", ), "isq_gdp": dict( file="isq_pib_regional_ensemble_ra.xlsx", url="https://statistique.quebec.ca/docs-ken/multimedia/pib-regional-ensemble-ra.xlsx", source_name="ISQ - Produit intérieur brut aux prix de base par industrie, régions administratives, Québec " "(tableau no 2299, 'Fichier complet RA')", page="https://statistique.quebec.ca/fr/produit/tableau/produit-interieur-brut-aux-prix-de-base-par-industrie-regions-administratives-du-quebec", licence=LIC_ISQ, vintage="2024 estimated ('e'); 2023 estimated-revised ('er'); 2011-2022 revised; ISQ update 15 Sep 2026", notes="GDP at basic prices, thousands of current dollars, by industry (NAICS). Provincial totals from " "StatCan 36-10-0487-01 (Nov 2025) for 2011-2022; 'X' = confidential (suppressed) cells left empty.", ), "isq_eqami2016": dict( file="isq_eqami_2016_ensemble_resultats.xlsx", url="https://statistique.quebec.ca/fr/fichier/xls-ensemble-des-resultats-de-lenquete-quebecoise-sur-lacces-des-menages-a-internet-2016.xlsx", source_name="ISQ - Enquête québécoise sur l'accès des ménages à Internet (EQAMI) 2016, ensemble des résultats", page="https://statistique.quebec.ca/fr/document/acces-a-internet", licence=LIC_ISQ, vintage="2016 (and 2012 comparisons)", notes="Checked: 2016 tables are broken down by income quartile, children, education and RMR / hors RMR only; " "NO administrative-region breakdown in 2016. Kept as evidence of that limitation.", ), "isq_eqami2012_ra": dict( file="isq_eqami_menages_branches_situation_geographique_tableau_3532.csv", url="https://statistique.quebec.ca/pls/ken/ken411_data_explt_v2.p_retrn_data", params={"p_id_tabl": "3532", "p_champs": "co_tertr,de_codfc_refrn,an_2012,an_2012_sign,an_2016,an_2016_sign"}, source_name="ISQ - Proportion des ménages branchés à Internet, selon la situation géographique " "(régions administratives et MRC), Québec, 2012 et 2016 (tableau dynamique no 3532)", page="https://statistique.quebec.ca/fr/document/acces-a-internet/tableau/proportion-des-menages-branches-a-internet-selon-la-situation-geographique-quebec-2012-et-2016", licence=LIC_ISQ, vintage="2012 (regional values); 2016 column is '..' (not available) at RA/MRC level", notes="Only regional (RA-level) household Internet-access source found. Côte-Nord and Nord-du-Québec are " "published combined ('09-10'). Quality code 'A' = excellent.", ), # ---------------- Statistics Canada (WDS full-table CSV) ---------------- "sc_98100009": dict( file="statcan_98100009-eng.zip", url="https://www150.statcan.gc.ca/n1/tbl/csv/98100009-eng.zip", source_name="Statistics Canada, Table 98-10-0009-01 Population and dwelling counts: Canada, provinces and " "territories, and economic regions (Census 2021)", page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=9810000901", licence=LIC_STATCAN, vintage="Census 2021 (released 9 Feb 2022); 2021 Census economic-region boundaries", notes="Land area (km2) and population density by economic region. CAUTION: 2021 ER boundaries for Estrie " "(2430) and Montérégie (2435) predate the 28 July 2021 RA boundary change (see README).", ), "sc_14100464": dict( file="statcan_14100464-eng.zip", url="https://www150.statcan.gc.ca/n1/tbl/csv/14100464-eng.zip", source_name="Statistics Canada, Table 14-10-0464-01 Labour force characteristics by province, territory and " "economic region, annual", page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1410046401", licence=LIC_STATCAN, vintage="2025 annual averages (released 9 Jan 2026)", notes="Labour Force Survey; persons in thousands; population 15+. Côte-Nord and Nord-du-Québec are a single " "LFS economic region ('Côte-Nord and Nord-du-Québec').", ), "sc_14100466": dict( file="statcan_14100466-eng.zip", url="https://www150.statcan.gc.ca/n1/tbl/csv/14100466-eng.zip", source_name="Statistics Canada, Table 14-10-0466-01 Employment characteristics by economic region, annual", page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1410046601", licence=LIC_STATCAN, vintage="2025 annual averages (released 9 Jan 2026)", notes="LFS employment by industry (NAICS sectors); manufacturing and 'wholesale and retail trade' extracted. " "Retail alone (44-45) is not separated from wholesale (41) in this table.", ), "sc_22100134": dict( file="statcan_22100134-eng.zip", url="https://www150.statcan.gc.ca/n1/tbl/csv/22100134-eng.zip", source_name="Statistics Canada, Table 22-10-0134-01 Access to the Internet at home by geography (CIUS 2022)", page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=2210013401", licence=LIC_STATCAN, vintage="2022", notes="Evidence that the Canadian Internet Use Survey is published only for province x (CMA/CA vs outside " "CMA/CA); no Quebec administrative-region breakdown. Not used for regional cells.", ), # ---------------- MRNF / Données Québec ---------------- "mern_bdat": dict( file="mern_bdat_adm_100k_shp.zip", url="https://diffusion.mern.gouv.qc.ca/diffusion/RGQ/Vectoriel/Theme/Regional/SDA_100k/SHP/BDAT(adm)_SHP.zip", source_name="MRNF (Québec) - Découpages administratifs 1/100 000 (BDAT adm), Shapefile; catalogued on Données Québec", page="https://www.donneesquebec.ca/recherche/dataset/decoupages-administratifs", licence=LIC_MERN, vintage="V2026-10 (file last modified 27 Oct 2026 on diffusion.mern.gouv.qc.ca; catalogue modified 2026-10-27)", notes="Layers regio_s (régions administratives, 21 polygons for 17 regions) and mrc_s (MRC) extracted to " "geo/qc_regions_admin_mern_bdat100k.gpkg; the 87 MB archive itself is deleted after extraction unless " "--keep-archives (its sha256 is recorded here). CRS EPSG:4269 (NAD83).", delete_after_extract=True, ), # ---------------- ATN (ex-CEFRIO) ---------------- "atn_netendances": dict( file="atn_netendances_2023_commerce_electronique.pdf", url="https://transformation-numerique.ulaval.ca/wp-content/uploads/2024/04/netendances-2023-commerce-electronique.pdf", source_name="Académie de la transformation numérique (Université Laval), NETendances 2023 - Le commerce " "électronique au Québec", page="https://transformation-numerique.ulaval.ca/", licence=LIC_ATN, vintage="Collection January 2024 (edition 2023)", notes="Regional results are published for 5 grouped zones (RMR Montréal, RMR Québec, Est, Centre, Ouest du " "Québec), not for the 17 administrative regions. Not used for regional cells; kept for citation.", ), } # StatCan WDS coordinate pulls (population estimates by ER, 2021 boundaries) WDS_POP = dict( key="sc_17100150", file="statcan_17100150_quebec_er_total_population.json", product_id=17100150, members=list(range(21, 39)), # 21 = Quebec, 22..38 = the 17 Quebec ERs source_name="Statistics Canada, Table 17-10-0150-01 Population estimates, July 1, by economic region, 2021 boundaries", page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1710015001", licence=LIC_STATCAN, vintage="1 July 2025 (released 14 Jan 2026); 2021 Census ER boundaries", notes="Pulled through the WDS endpoint getDataFromCubePidCoordAndLatestNPeriods (total gender, all ages, " "latest 5 periods). Same underlying estimates as the ISQ table, but on 2021 ER boundaries (Estrie / " "Montérégie differ from the ISQ RA figures for that reason).", ) # StatCan 2021 cartographic boundary files served by the official Esri REST service REST_BASE = "https://geo.statcan.gc.ca/geo_wa/rest/services/2021/Cartographic_boundary_files/MapServer" REST_LAYERS = { "er": dict(layer=2, fields="ERUID,DGUID,ERNAME,LANDAREA,PRUID", id_field="ERUID", gpkg="qc_er_2021_statcan_cbf.gpkg", source_name="Statistics Canada, 2021 Census Economic Region cartographic boundary file (ler_000b21a_e), " "Quebec subset via geo.statcan.gc.ca REST layer 2"), "cd": dict(layer=4, fields="CDUID,DGUID,CDNAME,CDTYPE,LANDAREA,PRUID", id_field="CDUID", gpkg=None, # polygons not kept (size); only representative points + concordance source_name="Statistics Canada, 2021 Census Division cartographic boundary file (lcd_000b21a_e), " "Quebec subset via geo.statcan.gc.ca REST layer 4"), "fsa": dict(layer=14, fields="CFSAUID,DGUID,PRUID,PRNAME,LANDAREA", id_field="CFSAUID", gpkg="qc_fsa_2021_statcan_cbf.gpkg", source_name="Statistics Canada, 2021 Census Forward Sortation Area cartographic boundary file " "(lfsa000b21a_e), Quebec subset via geo.statcan.gc.ca REST layer 14"), } WWW12_FSA_ZIP = ("https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/boundary-limites/" "files-fichiers/lfsa000b21a_e.zip") WWW12_ER_DEF = "https://www12.statcan.gc.ca/census-recensement/2021/ref/dict/az/Definition-eng.cfm?ID=geo022" SGC2016_INTRO = "https://www.statcan.gc.ca/en/subjects/standard/sgc/2016/introduction" # -------------------------------------------------------------------------------------- # helpers # -------------------------------------------------------------------------------------- def utcnow() -> str: return datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ") def sha256_path(p: Path) -> str: h = hashlib.sha256() with open(p, "rb") as f: for chunk in iter(lambda: f.read(1 << 20), b""): h.update(chunk) return h.hexdigest() def sha256_bytes(b: bytes) -> str: return hashlib.sha256(b).hexdigest() def log(msg: str) -> None: print(f"[external] {msg}", flush=True) def load_manifest() -> dict: if MANIFEST.exists(): return json.loads(MANIFEST.read_text(encoding="utf-8")) return {"generated_utc": None, "files": []} def http_get(url: str, params=None, timeout=300, retries=3) -> requests.Response: last = None for i in range(retries): try: r = SESSION.get(url, params=params, timeout=timeout) if r.status_code == 200: return r last = RuntimeError(f"HTTP {r.status_code} for {r.url}") except requests.RequestException as e: # noqa: PERF203 last = e time.sleep(2 * (i + 1)) raise last class Manifest: def __init__(self, old: dict, strict: bool): self.old = {e["file"]: e for e in old.get("files", [])} self.entries: list[dict] = [] self.strict = strict self.warnings: list[str] = [] def add(self, file: str, **kw) -> None: kw["file"] = file prev = self.old.get(file) if prev and "sha256" in kw and prev.get("sha256") and prev["sha256"] != kw["sha256"] and file.startswith("dl/"): # only raw downloads are checked; derived files legitimately change when the script changes msg = f"sha256 changed for {file}: {prev['sha256'][:12]}... -> {kw['sha256'][:12]}..." self.warnings.append(msg) log("WARNING " + msg) kw["previous_sha256"] = prev["sha256"] kw["previous_retrieved_utc"] = prev.get("retrieved_utc") if self.strict: raise SystemExit("--strict: checksum mismatch, aborting") self.entries.append(kw) def write(self, access_log: dict) -> None: out = {"generated_utc": utcnow(), "project": "UQO WP38", "note": "Every entry is a file in data/external (relative path). sha256 is of the file as stored. " "For GeoPackages the file hash changes on each rebuild (gpkg_contents.last_change timestamp), " "so content_sha256 (sorted feature ids + WKB) is the stable checksum.", "warnings": self.warnings, "files": self.entries} MANIFEST.write_text(json.dumps(out, ensure_ascii=False, indent=2), encoding="utf-8") ACCESS_LOG.write_text(json.dumps(access_log, ensure_ascii=False, indent=2), encoding="utf-8") def fetch_source(key: str, spec: dict, mf: Manifest, refresh: bool) -> Path: p = DL / spec["file"] retrieved = None if refresh or not p.exists(): log(f"downloading {key}: {spec['url']}") r = http_get(spec["url"], params=spec.get("params")) p.write_bytes(r.content) retrieved = utcnow() else: prev = mf.old.get(f"dl/{spec['file']}", {}) retrieved = prev.get("retrieved_utc") or datetime.fromtimestamp(p.stat().st_mtime, timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ") log(f"cached {key}: {p.name}") full_url = spec["url"] + ("?" + "&".join(f"{k}={requests.utils.quote(str(v))}" for k, v in spec["params"].items()) if spec.get("params") else "") mf.add(f"dl/{spec['file']}", source_name=spec["source_name"], url=full_url, landing_page=spec.get("page"), retrieved_utc=retrieved, licence=spec["licence"], vintage=spec["vintage"], notes=spec["notes"], sha256=sha256_path(p), bytes=p.stat().st_size) return p # -------------------------------------------------------------------------------------- # parsers (each returns {region_code: value}) -- all values parsed, never typed # -------------------------------------------------------------------------------------- def to_num(x): """ISQ/StatCan cell -> float or None. '..', 'X', 'x', 'F', '' -> None.""" if x is None or (isinstance(x, float) and pd.isna(x)): return None if isinstance(x, (int, float)): return float(x) s = str(x).strip().replace(" ", "").replace(" ", "") if s in {"", "..", "...", "X", "x", "F", "—", "–", "-", "nan", "NaN"}: return None s = s.replace(",", ".") s = re.sub(r"[^\d.\-]", "", s) # strip 'p','r','ᵖ','ʳ' flags try: return float(s) except ValueError: return None def parse_isq_population(p: Path) -> tuple[dict, dict, str]: df = pd.read_excel(p, header=None) hdr = df.iloc[1].tolist() # latest column = last header that starts with 4 digits cols = [(i, str(h)) for i, h in enumerate(hdr) if re.match(r"^\d{4}", str(h))] i_last, label = cols[-1] i_2021 = next(i for i, h in cols if str(h).startswith("2021")) out, out21 = {}, {} for _, row in df.iloc[3:].iterrows(): code = row[0] if isinstance(code, (int, float)) and not pd.isna(code) and 1 <= int(code) <= 17: c = f"{int(code):02d}" assert CODE_BY_NAME[norm(row[1])] == c, (row[1], c) out[c] = to_num(row[i_last]); out21[c] = to_num(row[i_2021]) return out, out21, label # label like '2025p' def parse_isq_business(p: Path) -> dict: df = pd.read_excel(p, sheet_name="Qc RA_industrie", header=None) hdr = df.iloc[4].tolist() colmap = {} for i, h in enumerate(hdr): m = re.search(r"\((\d{2})\)\s*$", str(h)) if m: colmap[m.group(1)] = i assert len(colmap) == 17, colmap rows = {"all": "Ensemble des industries", "mfg": "Fabrication [31-33]", "retail": "Commerce de détail [44-45]"} out = {k: {} for k in rows} seen = set() for _, row in df.iterrows(): lab = str(row[0]).strip() for k, want in rows.items(): if lab == want and k not in seen: # first occurrence = 2-digit NAICS block seen.add(k) for c, i in colmap.items(): out[k][c] = to_num(row[i]) assert seen == set(rows), seen return out def parse_isq_business_size(p: Path) -> dict: df = pd.read_excel(p, sheet_name="Qc par RA_taille", header=None) out = {} for _, row in df.iterrows(): m = re.search(r"\((\d{2})\)\s*$", str(row[0])) if m: out[m.group(1)] = to_num(row[1]) return out def parse_isq_income(p: Path) -> tuple[dict, dict]: rows = list(csv.reader(open(p, encoding="utf-8"), delimiter=";")) hdr = rows[0] pairs = [(i, i + 1) for i, h in enumerate(hdr) if h.startswith("titre_") and i + 1 < len(hdr) and hdr[i + 1].startswith("perd_")] best: dict[str, tuple[int, str, float]] = {} for r in rows[1:]: if r[2].strip() != "Revenu disponible" or r[3].strip() != "$/hab": continue code = r[0].strip() if code not in {c for c, _ in REGIONS}: continue for it, iv in pairs: lab = r[it].strip() m = re.match(r"^(\d{4})", lab) if not m: continue v = to_num(r[iv]) if v is None: continue y = int(m.group(1)) if code not in best or y > best[code][0]: best[code] = (y, lab, v) vals = {c: best[c][2] for c in best} labels = {c: best[c][1] for c in best} return vals, labels def parse_isq_gdp(p: Path) -> tuple[dict, str]: x = pd.ExcelFile(p) sheet = x.sheet_names[-1] # most recent year (e.g. '2024e') df = x.parse(sheet, header=None) hdr = df.iloc[2].tolist() colmap = {} for i, h in enumerate(hdr[1:], start=1): k = norm(h) if k in CODE_BY_NAME: colmap[CODE_BY_NAME[k]] = i assert len(colmap) == 17, colmap rows = {"all": "Ensemble des industries (T001)", "mfg": "Fabrication (31-33)", "retail": "Commerce de détail (44-45)", "goods": "Secteur de production de biens (T002)", "services": "Secteur des services (T003)"} out = {k: {} for k in rows} for _, row in df.iterrows(): lab = str(row[0]).strip() for k, want in rows.items(): if lab == want: for c, i in colmap.items(): out[k][c] = to_num(row[i]) for k in rows: assert out[k], k return out, sheet def parse_isq_eqami_ra(p: Path) -> tuple[dict, dict]: rows = list(csv.reader(open(p, encoding="utf-8"), delimiter=";")) out, combined = {}, {} for r in rows[1:]: code, name = r[0].strip(), r[1].strip() k = norm(name) if code == "09-10": combined["09"] = to_num(r[2]); combined["10"] = to_num(r[2]) elif len(code) == 2 and k in CODE_BY_NAME and CODE_BY_NAME[k] == code and code not in out: out[code] = to_num(r[2]) # an_2012 return out, combined def read_statcan_zip(p: Path) -> pd.DataFrame: with zipfile.ZipFile(p) as z: name = [n for n in z.namelist() if n.endswith(".csv") and "MetaData" not in n][0] return pd.read_csv(io.BytesIO(z.read(name)), low_memory=False) def parse_census_er(p: Path) -> dict: d = read_statcan_zip(p) d.columns = [c.split(": ")[-1].split(" [")[0] if c.startswith("Population and") else c for c in d.columns] q = d[d["DGUID"].astype(str).str.startswith("2021S050024")] out = {"pop": {}, "area": {}, "dens": {}, "eruid": {}} for _, r in q.iterrows(): c = CODE_BY_NAME[norm(r["GEO"])] out["pop"][c] = to_num(r["Population, 2021"]) out["area"][c] = to_num(r["Land area in square kilometres, 2021"]) out["dens"][c] = to_num(r["Population density per square kilometre, 2021"]) out["eruid"][c] = str(r["DGUID"])[-4:] assert len(out["pop"]) == 17 return out def parse_lfs(p464: Path, p466: Path) -> tuple[dict, int]: d = read_statcan_zip(p464) year = int(d["REF_DATE"].max()) q = d[(d["REF_DATE"] == year) & d["GEO"].str.endswith(", Quebec")] out = {"emp": {}, "urate": {}, "erate": {}, "pop15": {}, "combined": {}} for _, r in q.iterrows(): geo = r["GEO"].replace(", Quebec", "") ch = r["Labour force characteristics"] key = {"Employment": "emp", "Unemployment rate": "urate", "Employment rate": "erate", "Population": "pop15"}.get(ch) if key is None: continue if geo == "Côte-Nord and Nord-du-Québec": out["combined"].setdefault(key, to_num(r["VALUE"])) continue out[key][CODE_BY_NAME[norm(geo)]] = to_num(r["VALUE"]) d2 = read_statcan_zip(p466) y2 = int(d2["REF_DATE"].max()); assert y2 == year q2 = d2[(d2["REF_DATE"] == year) & d2["GEO"].str.endswith(", Quebec")] out["mfg"], out["wrt"] = {}, {} for _, r in q2.iterrows(): geo = r["GEO"].replace(", Quebec", "") ch = r["Employment characteristics"] key = {"Manufacturing [31-33]": "mfg", "Wholesale and retail trade [41, 44-45]": "wrt"}.get(ch) if key is None: continue if geo == "Côte-Nord and Nord-du-Québec": out["combined"].setdefault(key, to_num(r["VALUE"])) continue out[key][CODE_BY_NAME[norm(geo)]] = to_num(r["VALUE"]) return out, year def fetch_wds_population(mf: Manifest, refresh: bool) -> tuple[dict, str]: p = DL / WDS_POP["file"] if refresh or not p.exists(): body = [{"productId": WDS_POP["product_id"], "coordinate": f"{m}.1.1.0.0.0.0.0.0.0", "latestN": 5} for m in WDS_POP["members"]] r = SESSION.post("https://www150.statcan.gc.ca/t1/wds/rest/getDataFromCubePidCoordAndLatestNPeriods", json=body, timeout=120) r.raise_for_status() meta = SESSION.post("https://www150.statcan.gc.ca/t1/wds/rest/getCubeMetadata", json=[{"productId": WDS_POP["product_id"]}], timeout=120).json()[0]["object"] members = {m["memberId"]: m["memberNameEn"] for m in meta["dimension"][0]["member"]} payload = {"retrieved_utc": utcnow(), "request": body, "geography_members": {str(m): members[m] for m in WDS_POP["members"]}, "cube_release_time": meta.get("releaseTime"), "data": r.json()} p.write_text(json.dumps(payload, ensure_ascii=False, indent=1), encoding="utf-8") retrieved = payload["retrieved_utc"] else: payload = json.loads(p.read_text(encoding="utf-8")); retrieved = payload["retrieved_utc"] out, ref = {}, None for item in payload["data"]: o = item["object"]; member = int(o["coordinate"].split(".")[0]) geo = payload["geography_members"][str(member)].replace(", Quebec", "") if geo == "Quebec": continue pt = o["vectorDataPoint"][-1] out[CODE_BY_NAME[norm(geo)]] = to_num(pt["value"]); ref = pt["refPer"] assert len(out) == 17 mf.add(f"dl/{WDS_POP['file']}", source_name=WDS_POP["source_name"], url="https://www150.statcan.gc.ca/t1/wds/rest/getDataFromCubePidCoordAndLatestNPeriods (POST, productId 17100150)", landing_page=WDS_POP["page"], retrieved_utc=retrieved, licence=WDS_POP["licence"], vintage=WDS_POP["vintage"], notes=WDS_POP["notes"] + f" Latest refPer in file: {ref}.", sha256=sha256_path(p), bytes=p.stat().st_size) return out, ref # -------------------------------------------------------------------------------------- # geography # -------------------------------------------------------------------------------------- def rest_layer_geojson(layer: int, fields: str, access_log: dict) -> dict: """Pull all Quebec (PRUID='24') features of a layer. The server returns HTTP 500 for a where-clause with geometry, so we page by objectIds (batches of 10) and fall back to single ids with geometryPrecision=6 for very large polygons (Côte-Nord, Nord-du-Québec...).""" ids = http_get(f"{REST_BASE}/{layer}/query", params={"where": "PRUID='24'", "returnIdsOnly": "true", "f": "json"}).json()["objectIds"] feats, notes = [], [] for i in range(0, len(ids), 10): chunk = ids[i:i + 10] got = None for variant in ({}, {"geometryPrecision": 6}): try: d = http_get(f"{REST_BASE}/{layer}/query", params={"objectIds": ",".join(map(str, chunk)), "outFields": fields, "returnGeometry": "true", "f": "geojson", **variant}, retries=2).json() if "features" in d: got = d["features"]; break except Exception: # noqa: BLE001 continue if got is None: # single ids with precision 6 got = [] for oid in chunk: d = http_get(f"{REST_BASE}/{layer}/query", params={"objectIds": oid, "outFields": fields, "returnGeometry": "true", "f": "geojson", "geometryPrecision": 6}).json() got += d["features"]; notes.append(f"objectId {oid} fetched singly with geometryPrecision=6") feats += got assert len(feats) == len(ids), (len(feats), len(ids)) access_log[f"rest_layer_{layer}"] = {"url": f"{REST_BASE}/{layer}/query", "features": len(feats), "notes": notes, "retrieved_utc": utcnow()} return {"type": "FeatureCollection", "features": feats} def content_hash(gdf, id_field: str) -> str: h = hashlib.sha256() for _, r in gdf.sort_values(id_field).iterrows(): h.update(str(r[id_field]).encode()); h.update(r.geometry.wkb) return h.hexdigest() def build_geo(mf: Manifest, access_log: dict, refresh: bool, no_geo: bool, keep_archives: bool) -> dict: import geopandas as gpd # heavy import kept local GEO.mkdir(parents=True, exist_ok=True) res: dict = {} # ---- MERN administrative regions (CC-BY 4.0) ---- spec = SOURCES["mern_bdat"] zp = DL / spec["file"] ra_gpkg = GEO / "qc_regions_admin_mern_bdat100k.gpkg" if refresh or not ra_gpkg.exists() or zp.exists(): if refresh or not zp.exists(): fetch_source("mern_bdat", spec, mf, refresh=True) else: fetch_source("mern_bdat", spec, mf, refresh=False) tmp = EXT / "_tmp_bdat" with zipfile.ZipFile(zp) as z: for n in z.namelist(): base = n.split("/")[-1] if base.startswith(("regio_s.", "mrc_s.", "munic_s.")) and not base.endswith((".lock", ".xml", ".qmd")): (tmp / base).parent.mkdir(parents=True, exist_ok=True) (tmp / base).write_bytes(z.read(n)) regio = gpd.read_file(tmp / "regio_s.shp") mrc = gpd.read_file(tmp / "mrc_s.shp") munic = gpd.read_file(tmp / "munic_s.shp") reg = regio.dissolve(by="RES_CO_REG", as_index=False)[["RES_CO_REG", "RES_NM_REG", "RES_CO_REF", "RES_CO_VER", "geometry"]] reg = reg.rename(columns={"RES_CO_REG": "region_code", "RES_NM_REG": "region_name_mern", "RES_CO_REF": "mern_ref", "RES_CO_VER": "mern_version"}) assert len(reg) == 17 and set(reg.region_code) == {c for c, _ in REGIONS} reg["region_name"] = reg.region_code.map(dict(REGIONS)) if ra_gpkg.exists(): ra_gpkg.unlink() reg.to_file(ra_gpkg, layer="regions_admin", driver="GPKG") regio.to_file(ra_gpkg, layer="regio_s_raw", driver="GPKG") mrc_tab = mrc[["MRS_CO_MRC", "MRS_NM_MRC", "MRS_CO_REG", "MRS_NM_REG", "MRS_DE_IND", "geometry"]] mrc_tab.to_file(ra_gpkg, layer="mrc_s_raw", driver="GPKG") mun_cols = [c for c in munic.columns if c.startswith("MUS_CO_") or c.startswith("MUS_NM_") or c == "MUS_DE_IND"] munic[mun_cols].to_csv(EXT / "mern_municipalites_attributes.csv", index=False, encoding="utf-8") for f in tmp.iterdir(): f.unlink() tmp.rmdir() side = DL / (spec["file"] + ".sha256") side.write_text(json.dumps({"sha256": sha256_path(zp), "bytes": zp.stat().st_size, "retrieved_utc": [e for e in mf.entries if e["file"] == f"dl/{spec['file']}"][-1]["retrieved_utc"]}), encoding="utf-8") if spec.get("delete_after_extract") and not keep_archives: zp.unlink(); log("deleted archive mern_bdat_adm_100k_shp.zip after extraction (sha256 kept in manifest + .sha256 sidecar)") else: side = DL / (spec["file"] + ".sha256") prev = mf.old.get(f"dl/{spec['file']}") or (json.loads(side.read_text(encoding="utf-8")) if side.exists() else None) if prev: mf.add(f"dl/{spec['file']}", source_name=spec["source_name"], url=spec["url"], landing_page=spec["page"], retrieved_utc=prev.get("retrieved_utc"), licence=spec["licence"], vintage=spec["vintage"], notes=spec["notes"] + " (archive deleted after extraction; checksum from .sha256 sidecar)", sha256=prev["sha256"], bytes=prev.get("bytes")) reg = gpd.read_file(ra_gpkg, layer="regions_admin") mrc = gpd.read_file(ra_gpkg, layer="mrc_s_raw") mf.add(f"geo/{ra_gpkg.name}", source_name=spec["source_name"] + " -> layers regions_admin (dissolved, 17), regio_s_raw, mrc_s_raw", url=spec["url"], landing_page=spec["page"], retrieved_utc=utcnow(), licence=spec["licence"], vintage=spec["vintage"], notes="Derived by src/ingest/external.py (dissolve by RES_CO_REG). CRS EPSG:4269.", sha256=sha256_path(ra_gpkg), content_sha256=content_hash(reg, "region_code"), bytes=ra_gpkg.stat().st_size, features=len(reg)) reg_lam = reg.to_crs(6622) # NAD83(CSRS) / Quebec Lambert, metres res["area_mern_km2"] = dict(zip(reg_lam.region_code, (reg_lam.geometry.area / 1e6).round(2))) res["mern_version"] = reg.mern_version.iloc[0] # municipalities attribute table: municipality code -> region (for CSD-level concordance) res["mern_mrc_regions"] = mrc.groupby("MRS_CO_MRC")["MRS_CO_REG"].agg(lambda s: sorted(set(s))).to_dict() # ---- StatCan 2021 layers (Open Licence) ---- if no_geo: log("--no-geo: skipping StatCan REST pulls; using existing GeoPackages if present") layers = {} rest_notes = ("Quebec subset (PRUID=24). Direct zip on www12.statcan.gc.ca returns HTTP 403 to non-browser clients " "(see access_log.json); the REST service is the official alternative listed on open.canada.ca. " "CRS EPSG:4326 as returned by f=geojson. LANDAREA attribute is StatCan's land area in km2. " "Very large polygons (e.g. Côte-Nord, Nord-du-Québec, FSA G0G) were fetched singly with " "geometryPrecision=6 (~0.1 m) because the server returns HTTP 500 at full precision.") for key, L in REST_LAYERS.items(): gp = GEO / (L["gpkg"] or f"qc_{key}_2021_statcan_cbf_points.gpkg") kind = "polygons" if L["gpkg"] else "representative points only (polygons not kept to save space)" if gp.exists() and not refresh: layers[key] = gpd.read_file(gp) prev = mf.old.get(f"geo/{gp.name}", {}) retrieved = prev.get("retrieved_utc") or datetime.fromtimestamp(gp.stat().st_mtime, timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ") log(f"cached geo layer {key}: {gp.name}") elif no_geo: continue else: gj = rest_layer_geojson(L["layer"], L["fields"], access_log) g = gpd.GeoDataFrame.from_features(gj["features"], crs="EPSG:4326") if not L["gpkg"]: g["geometry"] = g.to_crs(6622).representative_point().to_crs(4326) if gp.exists(): gp.unlink() g.to_file(gp, layer=key, driver="GPKG") layers[key] = g; retrieved = utcnow() mf.add(f"geo/{gp.name}", source_name=L["source_name"] + f" ({kind})", url=f"{REST_BASE}/{L['layer']}/query?where=PRUID='24' (paged by objectIds, f=geojson)", landing_page="https://open.canada.ca/data/dataset/ef70dc3b-1069-4037-9bce-61f47e628a1d", retrieved_utc=retrieved, licence=LIC_STATCAN, vintage="2021 Census cartographic boundary files (l***000b21a_e), release 2022", notes=rest_notes, sha256=sha256_path(gp), content_sha256=content_hash(layers[key], L["id_field"]), bytes=gp.stat().st_size, features=len(layers[key])) if "er" not in layers or "fsa" not in layers or "cd" not in layers: log("geo layers incomplete (run without --no-geo once); skipping concordance products") return res er, cd, fsa = layers["er"], layers["cd"], layers["fsa"] assert len(er) == 17, len(er) er["region_code"] = er.ERNAME.map(lambda n: CODE_BY_NAME[norm(n)]) res["er_names"] = dict(zip(er.region_code, er.ERNAME)); res["eruid"] = dict(zip(er.region_code, er.ERUID)) res["er_landarea_km2"] = dict(zip(er.region_code, er.LANDAREA)) # ---- RA (MERN) vs ER (StatCan): measured overlap, not assumed ---- # make_valid: both sources contain self-touching rings (GEOS 'side location conflict'); areas are unaffected. er["geometry"] = er.geometry.make_valid(); reg["geometry"] = reg.geometry.make_valid(); fsa["geometry"] = fsa.geometry.make_valid() er_l, reg_l = er.to_crs(6622), reg.to_crs(6622) rows = [] for _, r in reg_l.iterrows(): e = er_l[er_l.region_code == r.region_code].iloc[0] inter = r.geometry.intersection(e.geometry).area; union = r.geometry.union(e.geometry).area rows.append(dict(region_code=r.region_code, region_name=r.region_name, eruid=e.ERUID, er_name=e.ERNAME, area_mern_km2=round(r.geometry.area / 1e6, 1), area_er_km2=round(e.geometry.area / 1e6, 1), intersection_km2=round(inter / 1e6, 1), iou=round(inter / union, 4), share_of_mern_in_er=round(inter / r.geometry.area, 4), share_of_er_in_mern=round(inter / e.geometry.area, 4))) ov = pd.DataFrame(rows).sort_values("region_code") ov.to_csv(EXT / "ra_er_overlap_2021.csv", index=False, encoding="utf-8") res["iou"] = dict(zip(ov.region_code, ov.iou)); res["er_in_mern"] = dict(zip(ov.region_code, ov.share_of_er_in_mern)) # ---- CD -> ER (point in StatCan ER) and CD -> RA (point in MERN RA) ---- cd_pts = cd.copy() if (cd_pts.geom_type != "Point").any(): # fresh polygons -> representative points cd_pts["geometry"] = cd_pts.to_crs(6622).representative_point().to_crs(4326) j1 = gpd.sjoin(cd_pts[["CDUID", "CDNAME", "CDTYPE", "LANDAREA", "geometry"]], er[["ERUID", "ERNAME", "region_code", "geometry"]], how="left", predicate="within").rename(columns={"region_code": "er_region_code"}).drop(columns="index_right") j2 = gpd.sjoin(j1, reg[["region_code", "region_name", "geometry"]].to_crs(4326), how="left", predicate="within").rename(columns={"region_code": "mern_region_code", "region_name": "mern_region_name"}).drop(columns="index_right") j2["mrc_code_mern"] = j2.CDUID.str[2:] j2["same_region"] = j2.er_region_code == j2.mern_region_code j2["lon"] = j2.geometry.x.round(5); j2["lat"] = j2.geometry.y.round(5) conc = j2.drop(columns="geometry").sort_values("CDUID") conc.to_csv(EXT / "cd_to_region_concordance_2021.csv", index=False, encoding="utf-8") res["cd_mismatch"] = conc[~conc.same_region][["CDUID", "CDNAME", "er_region_code", "mern_region_code"]].to_dict("records") # ---- FSA -> region: centroid + dominant-area assignment ---- f_l = fsa.to_crs(6622); reg_l = reg.to_crs(6622); er_l = er.to_crs(6622) cen = f_l.copy(); cen["geometry"] = f_l.representative_point() a = gpd.sjoin(cen[["CFSAUID", "LANDAREA", "geometry"]], reg_l[["region_code", "region_name", "geometry"]], how="left", predicate="within").drop(columns="index_right") a = a.rename(columns={"region_code": "region_code_centroid", "region_name": "region_name_centroid"}) b = gpd.sjoin(cen[["CFSAUID", "geometry"]], er_l[["region_code", "ERUID", "geometry"]], how="left", predicate="within").drop(columns="index_right") a["er_region_code_centroid"] = a.CFSAUID.map(dict(zip(b.CFSAUID, b.region_code))) a["eruid_centroid"] = a.CFSAUID.map(dict(zip(b.CFSAUID, b.ERUID))) # dominant area share ov2 = gpd.overlay(f_l[["CFSAUID", "geometry"]], reg_l[["region_code", "geometry"]], how="intersection", keep_geom_type=False) ov2["a"] = ov2.geometry.area tot = ov2.groupby("CFSAUID")["a"].sum() ov2["share"] = ov2.a / ov2.CFSAUID.map(tot) dom = ov2.sort_values("share", ascending=False).drop_duplicates("CFSAUID").set_index("CFSAUID") a["region_code_dominant_area"] = a.CFSAUID.map(dom.region_code) a["dominant_area_share"] = a.CFSAUID.map(dom.share).round(4) a["n_regions_touched"] = a.CFSAUID.map(ov2[ov2.share > 0.001].groupby("CFSAUID").size()) pt = cen.to_crs(4326) a["lon"] = a.CFSAUID.map(dict(zip(pt.CFSAUID, pt.geometry.x.round(5)))); a["lat"] = a.CFSAUID.map(dict(zip(pt.CFSAUID, pt.geometry.y.round(5)))) a = a.drop(columns="geometry").sort_values("CFSAUID") a.to_csv(EXT / "fsa_to_region_2021.csv", index=False, encoding="utf-8") res["fsa_n"] = len(a); res["fsa_multi"] = int((a.n_regions_touched > 1).sum()) res["fsa_centroid_vs_dominant_differ"] = int((a.region_code_centroid != a.region_code_dominant_area).sum()) res["fsa_count_by_region_centroid"] = a.groupby("region_code_centroid").size().to_dict() for name, df in (("ra_er_overlap_2021.csv", ov), ("cd_to_region_concordance_2021.csv", conc), ("fsa_to_region_2021.csv", a)): p = EXT / name mf.add(name, source_name="Derived by src/ingest/external.py from the MERN and StatCan layers above", url="(derived)", retrieved_utc=utcnow(), licence="Derived product; inherits CC-BY 4.0 (MRNF) and StatCan Open Licence", vintage="MERN V2026-10 x StatCan 2021 boundaries", notes="See README.md", sha256=sha256_path(p), bytes=p.stat().st_size, rows=len(df)) return res # -------------------------------------------------------------------------------------- # access probes (recorded, so the README claims are reproducible) # -------------------------------------------------------------------------------------- def probe(url: str, access_log: dict, key: str) -> None: try: r = SESSION.head(url, timeout=60, allow_redirects=True) access_log[key] = {"url": url, "http_status": r.status_code, "probed_utc": utcnow()} except requests.RequestException as e: access_log[key] = {"url": url, "http_status": None, "error": str(e), "probed_utc": utcnow()} # -------------------------------------------------------------------------------------- # main # -------------------------------------------------------------------------------------- def main() -> int: ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--refresh", action="store_true", help="force re-download of every source") ap.add_argument("--no-geo", action="store_true", help="skip StatCan REST geometry pulls (reuse existing GeoPackages)") ap.add_argument("--strict", action="store_true", help="abort if a re-downloaded file's sha256 differs from manifest") ap.add_argument("--keep-archives", action="store_true", help="keep the 87 MB MERN zip after extraction") args = ap.parse_args() EXT.mkdir(parents=True, exist_ok=True); DL.mkdir(exist_ok=True) mf = Manifest(load_manifest(), strict=args.strict) access_log = json.loads(ACCESS_LOG.read_text(encoding="utf-8")) if ACCESS_LOG.exists() else {} # 1) tabular sources paths = {k: fetch_source(k, s, mf, args.refresh) for k, s in SOURCES.items() if k != "mern_bdat"} pop_isq, pop_isq_2021, pop_label = parse_isq_population(paths["isq_pop"]) biz = parse_isq_business(paths["isq_biz_ind"]) biz_size_total = parse_isq_business_size(paths["isq_biz_size"]) assert all(abs(biz["all"][c] - biz_size_total[c]) < 0.5 for c in biz["all"]), "ISQ business totals differ between files" inc, inc_lab = parse_isq_income(paths["isq_income"]) gdp, gdp_sheet = parse_isq_gdp(paths["isq_gdp"]) inet, inet_comb = parse_isq_eqami_ra(paths["isq_eqami2012_ra"]) cen = parse_census_er(paths["sc_98100009"]) lfs, lfs_year = parse_lfs(paths["sc_14100464"], paths["sc_14100466"]) pop_sc, pop_sc_ref = fetch_wds_population(mf, args.refresh) # 2) access probes probe(WWW12_FSA_ZIP, access_log, "www12_fsa_zip") probe(WWW12_ER_DEF, access_log, "www12_census_dictionary_er_definition") probe(SGC2016_INTRO, access_log, "statcan_sgc2016_introduction") probe(SOURCES["mern_bdat"]["url"], access_log, "mern_bdat_zip") # 3) geography geo = build_geo(mf, access_log, args.refresh, args.no_geo, args.keep_archives) # 4) assemble tidy table pop_year = re.match(r"\d{4}", pop_label).group(0) rows = [] for code, name in REGIONS: r = {"region_code": code, "region_name": name} def put(col, val, vintage, source): if isinstance(val, float) and val.is_integer(): val = int(val) r[col] = "" if val is None else val r[f"{col}_vintage"] = vintage if val is not None else "" r[f"{col}_source"] = source if val is not None else "" put("population", pop_isq.get(code), f"1 July {pop_year} ({pop_label}: provisional)", "ISQ table 3595 (StatCan annual demographic estimates adapted by ISQ)") put("population_statcan_er", pop_sc.get(code), f"{pop_sc_ref[:4]}-07-01 (2021 ER boundaries)", "StatCan 17-10-0150-01 (WDS)") put("population_census_2021_er", cen["pop"].get(code), "Census 2021 (2021 ER boundaries)", "StatCan 98-10-0009-01") put("business_locations_with_employees", biz["all"].get(code), "Dec 2025", "ISQ (StatCan 33-10-1097-01 aggregated to RA)") put("business_locations_with_employees_mfg_31_33", biz["mfg"].get(code), "Dec 2025", "ISQ (StatCan 33-10-1097-01 aggregated to RA)") put("business_locations_with_employees_retail_44_45", biz["retail"].get(code), "Dec 2025", "ISQ (StatCan 33-10-1097-01 aggregated to RA)") put("disposable_income_per_capita_cad", inc.get(code), inc_lab.get(code, "").replace("ᵖ", "p").replace("ʳ", "r"), "ISQ table 3962 (comptes économiques régionaux)") put("gdp_basic_prices_kcad", gdp["all"].get(code), gdp_sheet, "ISQ table 2299 (PIB régional)") put("gdp_mfg_31_33_kcad", gdp["mfg"].get(code), gdp_sheet, "ISQ table 2299 (PIB régional)") put("gdp_retail_44_45_kcad", gdp["retail"].get(code), gdp_sheet, "ISQ table 2299 (PIB régional)") put("employment_lfs_thousands", lfs["emp"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0464-01") put("employment_lfs_mfg_thousands", lfs["mfg"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0466-01") put("employment_lfs_wholesale_retail_thousands", lfs["wrt"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0466-01") put("unemployment_rate_pct", lfs["urate"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0464-01") put("employment_rate_pct", lfs["erate"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0464-01") comb = lfs["combined"] if code in ("09", "10") else {} put("employment_lfs_thousands_combined_09_10", comb.get("emp"), f"{lfs_year} annual average (Côte-Nord + Nord-du-Québec combined ER)", "StatCan 14-10-0464-01") put("unemployment_rate_pct_combined_09_10", comb.get("urate"), f"{lfs_year} annual average (combined ER)", "StatCan 14-10-0464-01") put("land_area_km2", cen["area"].get(code), "Census 2021 (2021 ER boundaries)", "StatCan 98-10-0009-01") put("population_density_per_km2", cen["dens"].get(code), "Census 2021 (2021 ER boundaries)", "StatCan 98-10-0009-01") put("polygon_area_mern_km2", geo.get("area_mern_km2", {}).get(code), f"MERN BDAT 100k {geo.get('mern_version', '')}", "MRNF Découpages administratifs (area computed in EPSG:6622; includes inland water)") put("households_internet_access_pct", inet.get(code), "2012 (EQAMI)", "ISQ table 3532 (Enquête québécoise sur l'accès des ménages à Internet)") put("households_internet_access_pct_combined_09_10", inet_comb.get(code), "2012 (EQAMI, Côte-Nord + Nord-du-Québec combined)", "ISQ table 3532") r["statcan_eruid_2021"] = cen["eruid"].get(code, "") r["ra_er_boundary_iou_2021"] = geo.get("iou", {}).get(code, "") r["share_of_statcan_er_land_inside_mern_ra"] = geo.get("er_in_mern", {}).get(code, "") r["n_fsa_centroids_2021"] = geo.get("fsa_count_by_region_centroid", {}).get(code, "") rows.append(r) out = pd.DataFrame(rows) out.to_csv(OUT_CSV, index=False, encoding="utf-8") mf.add("regions_external.csv", source_name="Assembled by src/ingest/external.py from the files above", url="(derived)", retrieved_utc=utcnow(), licence="Derived; see each *_source column", vintage="see *_vintage columns", notes="Empty cell = not obtained (never zero).", sha256=sha256_path(OUT_CSV), bytes=OUT_CSV.stat().st_size, rows=len(out)) mf.write(access_log) # 5) console summary log(f"wrote {OUT_CSV.relative_to(PROJECT)} ({len(out)} rows, {len(out.columns)} columns)") for col in [c for c in out.columns if not c.endswith(("_vintage", "_source")) and c not in ("region_code", "region_name")]: n = int((out[col].astype(str) != "").sum()) log(f" {col:52s} filled {n:2d}/17") if geo.get("cd_mismatch") is not None: log(f"CD->region mismatches between StatCan 2021 ER and MERN RA: {geo['cd_mismatch']}") if "fsa_n" in geo: log(f"FSA: {geo['fsa_n']} Quebec FSAs; {geo['fsa_multi']} touch >1 region; centroid vs dominant-area differ for {geo['fsa_centroid_vs_dominant_differ']}") total = sum(p.stat().st_size for p in EXT.rglob("*") if p.is_file()) log(f"data/external total size: {total / 1e6:.1f} MB") if mf.warnings: log(f"{len(mf.warnings)} checksum warning(s) - see manifest.json") return 0 if __name__ == "__main__": sys.exit(main())