1#!/usr/bin/env python2"""3WP38 - External regional covariates for the 17 Quebec administrative regions.45Re-downloads (idempotently, with SHA-256 checksums) every authoritative public6source listed in SOURCES, derives the geography products, and rebuilds78 data/external/regions_external.csv (one row per region, tidy)9 data/external/manifest.json (one entry per file: url, sha256, vintage, licence, notes)10 data/external/geo/*.gpkg (MERN admin regions, StatCan 2021 ER / FSA for Quebec)11 data/external/fsa_to_region_2021.csv (FSA -> region assignment, centroid + dominant-area)12 data/external/cd_to_region_concordance_2021.csv13 data/external/ra_er_overlap_2021.csv (does StatCan ER == Quebec RA ? measured, not assumed)14 data/external/access_log.json (HTTP status of blocked / alternative endpoints)1516Rules (CLAUDE.md rule 1): nothing is typed in by hand. Every number in the CSV17is parsed from a downloaded file whose checksum is in the manifest. Cells that18could not be obtained stay EMPTY (never 0).1920Usage21-----22 .venv/bin/python src/ingest/external.py # download what is missing, rebuild everything23 .venv/bin/python src/ingest/external.py --refresh # force re-download of every source24 .venv/bin/python src/ingest/external.py --no-geo # skip the (slow) StatCan REST geometry pulls25 .venv/bin/python src/ingest/external.py --strict # fail if a re-downloaded file's sha256 changed2627Only official hosts are contacted: statistique.quebec.ca (ISQ), www150.statcan.gc.ca28(StatCan WDS), geo.statcan.gc.ca (StatCan boundary REST service), www12.statcan.gc.ca29(probe only - returns 403 to non-browser clients), diffusion.mern.gouv.qc.ca (MRNF,30via donneesquebec.ca catalogue), transformation-numerique.ulaval.ca (ATN / NETendances).31"""32from __future__ import annotations3334import argparse35import csv36import hashlib37import io38import json39import re40import sys41import time42import unicodedata43import zipfile44from datetime import datetime, timezone45from pathlib import Path4647import pandas as pd48import requests4950PROJECT = Path(__file__).resolve().parents[2]51EXT = PROJECT / "data" / "external"52DL = EXT / "dl"53GEO = EXT / "geo"54MANIFEST = EXT / "manifest.json"55ACCESS_LOG = EXT / "access_log.json"56OUT_CSV = EXT / "regions_external.csv"5758UA = ("Mozilla/5.0 (Macintosh; Intel Mac OS X 14_0) AppleWebKit/537.36 "59 "(KHTML, like Gecko) Chrome/124.0.0.0 Safari/537.36 WP38-UQO-research/1.0")60SESSION = requests.Session()61SESSION.headers.update({"User-Agent": UA})6263# --------------------------------------------------------------------------------------64# The 17 administrative regions (official codes and spellings; en-dash in 02 and 11)65# --------------------------------------------------------------------------------------66REGIONS = [67 ("01", "Bas-Saint-Laurent"), ("02", "Saguenay–Lac-Saint-Jean"), ("03", "Capitale-Nationale"),68 ("04", "Mauricie"), ("05", "Estrie"), ("06", "Montréal"), ("07", "Outaouais"),69 ("08", "Abitibi-Témiscamingue"), ("09", "Côte-Nord"), ("10", "Nord-du-Québec"),70 ("11", "Gaspésie–Îles-de-la-Madeleine"), ("12", "Chaudière-Appalaches"), ("13", "Laval"),71 ("14", "Lanaudière"), ("15", "Laurentides"), ("16", "Montérégie"), ("17", "Centre-du-Québec"),72]73CODE_BY_NAME = {}747576def norm(s: str) -> str:77 """Accent/dash/space-insensitive key used to match region names across sources."""78 s = str(s)79 s = s.replace("\n", " ")80 s = re.sub(r"[\u2013\u2014\-]+", "-", s) # en/em dash and double hyphen -> '-' (before ASCII fold)81 s = unicodedata.normalize("NFKD", s).encode("ascii", "ignore").decode()82 s = re.sub(r"\s+", " ", s).strip().lower()83 s = re.sub(r"\s*\(\d+\)\s*$", "", s) # trailing "(01)"84 s = re.sub(r"[\d,]+$", "", s) # ISQ footnote digits "Estrie2,3"85 s = s.replace(", quebec", "")86 s = s.replace(" ", "")87 return s888990for _c, _n in REGIONS:91 CODE_BY_NAME[norm(_n)] = _c92# StatCan spellings (double hyphen) and ISQ footnoted variants map through norm() already.9394# --------------------------------------------------------------------------------------95# Sources96# --------------------------------------------------------------------------------------97LIC_ISQ = ("ISQ - Licence d'utilisation des produits de l'Institut (reproduction permitted with "98 "attribution; see https://statistique.quebec.ca/fr/institut/politique-de-confidentialite "99 "and the 'Conditions d'utilisation' page)")100LIC_STATCAN = "Statistics Canada Open Licence (https://www.statcan.gc.ca/en/reference/licence)"101LIC_MERN = "Creative Commons Attribution 4.0 (CC-BY 4.0), as stated on donneesquebec.ca dataset page"102LIC_ATN = "Copyright Académie de la transformation numérique (Université Laval); report PDF kept for citation only"103104SOURCES = {105 # ---------------- ISQ ----------------106 "isq_pop": dict(107 file="isq_population_ra_1986_2025_fichier_complet_3595.xlsx",108 url="https://statistique.quebec.ca/docs-ken/multimedia/Fichier_complet_3595.xlsx",109 source_name="ISQ - Estimations de la population et indicateurs d'accroissement des régions administratives, "110 "Québec, 1er juillet 1986 à 2025 (tableau no 3595, 'Fichier complet')",111 page="https://statistique.quebec.ca/fr/produit/tableau/estimations-population-regions-administratives",112 licence=LIC_ISQ,113 vintage="1 July 2025 (provisional, 'p'); 2022-2024 revised; ISQ update 14 Jan 2026; RA boundaries as of 1 July 2025",114 notes="Source of the ISQ table: Statistics Canada, Annual Demographic Estimates (subprovincial areas, January 2026), "115 "adapted by ISQ. Estrie/Montérégie reflect the 28 July 2021 boundary change (MRC La Haute-Yamaska and "116 "Brome-Missisquoi moved from Montérégie to Estrie); Estrie/Chaudière-Appalaches reflect the 1 Jan 2024 "117 "Courcelles merger.",118 ),119 "isq_biz_ind": dict(120 file="isq_entreprises_actives_ra_industrie_dec2025.xlsx",121 url="https://statistique.quebec.ca/fr/fichier/entreprises-actives-regions-adminsitratives-industrie-decembre-2025-donnees-semstrielles.xlsx",122 source_name="ISQ - Nombre d'entreprises actives (avec employés) dans les régions administratives du Québec selon "123 "l'industrie (SCIAN 2-4 chiffres), décembre 2025 (données semestrielles)",124 page="https://statistique.quebec.ca/fr/document/nombre-entreprises-actives-quebec",125 licence=LIC_ISQ,126 vintage="December 2025 (Business Register snapshot); ISQ compilation of StatCan table 33-10-1097-01, data captured 2 March 2026",127 notes="Unit = business 'statistical locations' WITH EMPLOYEES only (employer locations). Locations without "128 "employees are published by ISQ at the province level only. ISQ aggregates StatCan census-subdivision "129 "counts (33-10-1097-01) to administrative regions; 703 employer locations are 'non classifiées par "130 "municipalité' and are excluded from the regional totals.",131 ),132 "isq_biz_size": dict(133 file="isq_entreprises_actives_ra_municipalites_dec2025.xlsx",134 url="https://statistique.quebec.ca/fr/fichier/entreprises-actives-regions-administratives-municipalites-decembre-2025-donnees-semestrielles.xlsx",135 source_name="ISQ - Nombre d'entreprises avec employés par région administrative et par tranche d'effectif, "136 "ensemble des industries, décembre 2025",137 page="https://statistique.quebec.ca/fr/document/nombre-entreprises-actives-quebec",138 licence=LIC_ISQ,139 vintage="December 2025",140 notes="Same universe as isq_biz_ind (employer locations); used as cross-check of regional totals and for the "141 "employment-size distribution.",142 ),143 "isq_income": dict(144 file="isq_revenu_disponible_ra_tableau_3962.csv",145 url="https://statistique.quebec.ca/pls/ken/ken411_data_explt_v2.p_retrn_data",146 params={"p_id_tabl": "3962",147 "p_champs": "co_tertr,tertr,de_codfc_refrn,mesr,titre_23,perd_23,titre_22,perd_22,titre_21,perd_21,"148 "titre_20,perd_20,titre_19,perd_19,titre_18,perd_18,titre_17,perd_17,titre_16,perd_16,"149 "perd_17_sign,titre_15,perd_15,titre_14,perd_14,titre_13,perd_13,titre_12,perd_12,titre_11,"150 "perd_11,titre_10,perd_10,titre_9,perd_9,titre_8,perd_8,titre_7,perd_7,titre_6,perd_6,"151 "titre_5,perd_5,titre_4,perd_4,titre_3,perd_3,titre_2,perd_2,titre_1,perd_1",152 "p_tri": "tri_tertr asc,tri_codfc asc"},153 source_name="ISQ - Revenu disponible et ses composantes, régions administratives et ensemble du Québec, "154 "2002-2024 (tableau dynamique no 3962; data endpoint used by the ISQ web table)",155 page="https://statistique.quebec.ca/fr/document/revenu-disponible/tableau/revenu-disponible-et-composantes-regions-administratives-et-ensemble-du-quebec",156 licence=LIC_ISQ,157 vintage="2024 provisional ('p'); 2022-2023 revised; ISQ release 1 June 2026",158 notes="Disposable income = household-sector aggregate of the System of National Accounts (ISQ regional "159 "economic accounts), per capita in current dollars ('$/hab'). Not comparable with MFQ tax statistics.",160 ),161 "isq_gdp": dict(162 file="isq_pib_regional_ensemble_ra.xlsx",163 url="https://statistique.quebec.ca/docs-ken/multimedia/pib-regional-ensemble-ra.xlsx",164 source_name="ISQ - Produit intérieur brut aux prix de base par industrie, régions administratives, Québec "165 "(tableau no 2299, 'Fichier complet RA')",166 page="https://statistique.quebec.ca/fr/produit/tableau/produit-interieur-brut-aux-prix-de-base-par-industrie-regions-administratives-du-quebec",167 licence=LIC_ISQ,168 vintage="2024 estimated ('e'); 2023 estimated-revised ('er'); 2011-2022 revised; ISQ update 15 Sep 2026",169 notes="GDP at basic prices, thousands of current dollars, by industry (NAICS). Provincial totals from "170 "StatCan 36-10-0487-01 (Nov 2025) for 2011-2022; 'X' = confidential (suppressed) cells left empty.",171 ),172 "isq_eqami2016": dict(173 file="isq_eqami_2016_ensemble_resultats.xlsx",174 url="https://statistique.quebec.ca/fr/fichier/xls-ensemble-des-resultats-de-lenquete-quebecoise-sur-lacces-des-menages-a-internet-2016.xlsx",175 source_name="ISQ - Enquête québécoise sur l'accès des ménages à Internet (EQAMI) 2016, ensemble des résultats",176 page="https://statistique.quebec.ca/fr/document/acces-a-internet",177 licence=LIC_ISQ,178 vintage="2016 (and 2012 comparisons)",179 notes="Checked: 2016 tables are broken down by income quartile, children, education and RMR / hors RMR only; "180 "NO administrative-region breakdown in 2016. Kept as evidence of that limitation.",181 ),182 "isq_eqami2012_ra": dict(183 file="isq_eqami_menages_branches_situation_geographique_tableau_3532.csv",184 url="https://statistique.quebec.ca/pls/ken/ken411_data_explt_v2.p_retrn_data",185 params={"p_id_tabl": "3532",186 "p_champs": "co_tertr,de_codfc_refrn,an_2012,an_2012_sign,an_2016,an_2016_sign"},187 source_name="ISQ - Proportion des ménages branchés à Internet, selon la situation géographique "188 "(régions administratives et MRC), Québec, 2012 et 2016 (tableau dynamique no 3532)",189 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",190 licence=LIC_ISQ,191 vintage="2012 (regional values); 2016 column is '..' (not available) at RA/MRC level",192 notes="Only regional (RA-level) household Internet-access source found. Côte-Nord and Nord-du-Québec are "193 "published combined ('09-10'). Quality code 'A' = excellent.",194 ),195 # ---------------- Statistics Canada (WDS full-table CSV) ----------------196 "sc_98100009": dict(197 file="statcan_98100009-eng.zip",198 url="https://www150.statcan.gc.ca/n1/tbl/csv/98100009-eng.zip",199 source_name="Statistics Canada, Table 98-10-0009-01 Population and dwelling counts: Canada, provinces and "200 "territories, and economic regions (Census 2021)",201 page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=9810000901",202 licence=LIC_STATCAN,203 vintage="Census 2021 (released 9 Feb 2022); 2021 Census economic-region boundaries",204 notes="Land area (km2) and population density by economic region. CAUTION: 2021 ER boundaries for Estrie "205 "(2430) and Montérégie (2435) predate the 28 July 2021 RA boundary change (see README).",206 ),207 "sc_14100464": dict(208 file="statcan_14100464-eng.zip",209 url="https://www150.statcan.gc.ca/n1/tbl/csv/14100464-eng.zip",210 source_name="Statistics Canada, Table 14-10-0464-01 Labour force characteristics by province, territory and "211 "economic region, annual",212 page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1410046401",213 licence=LIC_STATCAN,214 vintage="2025 annual averages (released 9 Jan 2026)",215 notes="Labour Force Survey; persons in thousands; population 15+. Côte-Nord and Nord-du-Québec are a single "216 "LFS economic region ('Côte-Nord and Nord-du-Québec').",217 ),218 "sc_14100466": dict(219 file="statcan_14100466-eng.zip",220 url="https://www150.statcan.gc.ca/n1/tbl/csv/14100466-eng.zip",221 source_name="Statistics Canada, Table 14-10-0466-01 Employment characteristics by economic region, annual",222 page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1410046601",223 licence=LIC_STATCAN,224 vintage="2025 annual averages (released 9 Jan 2026)",225 notes="LFS employment by industry (NAICS sectors); manufacturing and 'wholesale and retail trade' extracted. "226 "Retail alone (44-45) is not separated from wholesale (41) in this table.",227 ),228 "sc_22100134": dict(229 file="statcan_22100134-eng.zip",230 url="https://www150.statcan.gc.ca/n1/tbl/csv/22100134-eng.zip",231 source_name="Statistics Canada, Table 22-10-0134-01 Access to the Internet at home by geography (CIUS 2022)",232 page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=2210013401",233 licence=LIC_STATCAN,234 vintage="2022",235 notes="Evidence that the Canadian Internet Use Survey is published only for province x (CMA/CA vs outside "236 "CMA/CA); no Quebec administrative-region breakdown. Not used for regional cells.",237 ),238 # ---------------- MRNF / Données Québec ----------------239 "mern_bdat": dict(240 file="mern_bdat_adm_100k_shp.zip",241 url="https://diffusion.mern.gouv.qc.ca/diffusion/RGQ/Vectoriel/Theme/Regional/SDA_100k/SHP/BDAT(adm)_SHP.zip",242 source_name="MRNF (Québec) - Découpages administratifs 1/100 000 (BDAT adm), Shapefile; catalogued on Données Québec",243 page="https://www.donneesquebec.ca/recherche/dataset/decoupages-administratifs",244 licence=LIC_MERN,245 vintage="V2026-10 (file last modified 27 Oct 2026 on diffusion.mern.gouv.qc.ca; catalogue modified 2026-10-27)",246 notes="Layers regio_s (régions administratives, 21 polygons for 17 regions) and mrc_s (MRC) extracted to "247 "geo/qc_regions_admin_mern_bdat100k.gpkg; the 87 MB archive itself is deleted after extraction unless "248 "--keep-archives (its sha256 is recorded here). CRS EPSG:4269 (NAD83).",249 delete_after_extract=True,250 ),251 # ---------------- ATN (ex-CEFRIO) ----------------252 "atn_netendances": dict(253 file="atn_netendances_2023_commerce_electronique.pdf",254 url="https://transformation-numerique.ulaval.ca/wp-content/uploads/2024/04/netendances-2023-commerce-electronique.pdf",255 source_name="Académie de la transformation numérique (Université Laval), NETendances 2023 - Le commerce "256 "électronique au Québec",257 page="https://transformation-numerique.ulaval.ca/",258 licence=LIC_ATN,259 vintage="Collection January 2024 (edition 2023)",260 notes="Regional results are published for 5 grouped zones (RMR Montréal, RMR Québec, Est, Centre, Ouest du "261 "Québec), not for the 17 administrative regions. Not used for regional cells; kept for citation.",262 ),263}264265# StatCan WDS coordinate pulls (population estimates by ER, 2021 boundaries)266WDS_POP = dict(267 key="sc_17100150",268 file="statcan_17100150_quebec_er_total_population.json",269 product_id=17100150,270 members=list(range(21, 39)), # 21 = Quebec, 22..38 = the 17 Quebec ERs271 source_name="Statistics Canada, Table 17-10-0150-01 Population estimates, July 1, by economic region, 2021 boundaries",272 page="https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1710015001",273 licence=LIC_STATCAN,274 vintage="1 July 2025 (released 14 Jan 2026); 2021 Census ER boundaries",275 notes="Pulled through the WDS endpoint getDataFromCubePidCoordAndLatestNPeriods (total gender, all ages, "276 "latest 5 periods). Same underlying estimates as the ISQ table, but on 2021 ER boundaries (Estrie / "277 "Montérégie differ from the ISQ RA figures for that reason).",278)279280# StatCan 2021 cartographic boundary files served by the official Esri REST service281REST_BASE = "https://geo.statcan.gc.ca/geo_wa/rest/services/2021/Cartographic_boundary_files/MapServer"282REST_LAYERS = {283 "er": dict(layer=2, fields="ERUID,DGUID,ERNAME,LANDAREA,PRUID", id_field="ERUID",284 gpkg="qc_er_2021_statcan_cbf.gpkg",285 source_name="Statistics Canada, 2021 Census Economic Region cartographic boundary file (ler_000b21a_e), "286 "Quebec subset via geo.statcan.gc.ca REST layer 2"),287 "cd": dict(layer=4, fields="CDUID,DGUID,CDNAME,CDTYPE,LANDAREA,PRUID", id_field="CDUID",288 gpkg=None, # polygons not kept (size); only representative points + concordance289 source_name="Statistics Canada, 2021 Census Division cartographic boundary file (lcd_000b21a_e), "290 "Quebec subset via geo.statcan.gc.ca REST layer 4"),291 "fsa": dict(layer=14, fields="CFSAUID,DGUID,PRUID,PRNAME,LANDAREA", id_field="CFSAUID",292 gpkg="qc_fsa_2021_statcan_cbf.gpkg",293 source_name="Statistics Canada, 2021 Census Forward Sortation Area cartographic boundary file "294 "(lfsa000b21a_e), Quebec subset via geo.statcan.gc.ca REST layer 14"),295}296WWW12_FSA_ZIP = ("https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/boundary-limites/"297 "files-fichiers/lfsa000b21a_e.zip")298WWW12_ER_DEF = "https://www12.statcan.gc.ca/census-recensement/2021/ref/dict/az/Definition-eng.cfm?ID=geo022"299SGC2016_INTRO = "https://www.statcan.gc.ca/en/subjects/standard/sgc/2016/introduction"300301302# --------------------------------------------------------------------------------------303# helpers304# --------------------------------------------------------------------------------------305def utcnow() -> str:306 return datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")307308309def sha256_path(p: Path) -> str:310 h = hashlib.sha256()311 with open(p, "rb") as f:312 for chunk in iter(lambda: f.read(1 << 20), b""):313 h.update(chunk)314 return h.hexdigest()315316317def sha256_bytes(b: bytes) -> str:318 return hashlib.sha256(b).hexdigest()319320321def log(msg: str) -> None:322 print(f"[external] {msg}", flush=True)323324325def load_manifest() -> dict:326 if MANIFEST.exists():327 return json.loads(MANIFEST.read_text(encoding="utf-8"))328 return {"generated_utc": None, "files": []}329330331def http_get(url: str, params=None, timeout=300, retries=3) -> requests.Response:332 last = None333 for i in range(retries):334 try:335 r = SESSION.get(url, params=params, timeout=timeout)336 if r.status_code == 200:337 return r338 last = RuntimeError(f"HTTP {r.status_code} for {r.url}")339 except requests.RequestException as e: # noqa: PERF203340 last = e341 time.sleep(2 * (i + 1))342 raise last343344345class Manifest:346 def __init__(self, old: dict, strict: bool):347 self.old = {e["file"]: e for e in old.get("files", [])}348 self.entries: list[dict] = []349 self.strict = strict350 self.warnings: list[str] = []351352 def add(self, file: str, **kw) -> None:353 kw["file"] = file354 prev = self.old.get(file)355 if prev and "sha256" in kw and prev.get("sha256") and prev["sha256"] != kw["sha256"] and file.startswith("dl/"):356 # only raw downloads are checked; derived files legitimately change when the script changes357 msg = f"sha256 changed for {file}: {prev['sha256'][:12]}... -> {kw['sha256'][:12]}..."358 self.warnings.append(msg)359 log("WARNING " + msg)360 kw["previous_sha256"] = prev["sha256"]361 kw["previous_retrieved_utc"] = prev.get("retrieved_utc")362 if self.strict:363 raise SystemExit("--strict: checksum mismatch, aborting")364 self.entries.append(kw)365366 def write(self, access_log: dict) -> None:367 out = {"generated_utc": utcnow(), "project": "UQO WP38", "note":368 "Every entry is a file in data/external (relative path). sha256 is of the file as stored. "369 "For GeoPackages the file hash changes on each rebuild (gpkg_contents.last_change timestamp), "370 "so content_sha256 (sorted feature ids + WKB) is the stable checksum.",371 "warnings": self.warnings, "files": self.entries}372 MANIFEST.write_text(json.dumps(out, ensure_ascii=False, indent=2), encoding="utf-8")373 ACCESS_LOG.write_text(json.dumps(access_log, ensure_ascii=False, indent=2), encoding="utf-8")374375376def fetch_source(key: str, spec: dict, mf: Manifest, refresh: bool) -> Path:377 p = DL / spec["file"]378 retrieved = None379 if refresh or not p.exists():380 log(f"downloading {key}: {spec['url']}")381 r = http_get(spec["url"], params=spec.get("params"))382 p.write_bytes(r.content)383 retrieved = utcnow()384 else:385 prev = mf.old.get(f"dl/{spec['file']}", {})386 retrieved = prev.get("retrieved_utc") or datetime.fromtimestamp(p.stat().st_mtime, timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")387 log(f"cached {key}: {p.name}")388 full_url = spec["url"] + ("?" + "&".join(f"{k}={requests.utils.quote(str(v))}" for k, v in spec["params"].items()) if spec.get("params") else "")389 mf.add(f"dl/{spec['file']}", source_name=spec["source_name"], url=full_url, landing_page=spec.get("page"),390 retrieved_utc=retrieved, licence=spec["licence"], vintage=spec["vintage"], notes=spec["notes"],391 sha256=sha256_path(p), bytes=p.stat().st_size)392 return p393394395# --------------------------------------------------------------------------------------396# parsers (each returns {region_code: value}) -- all values parsed, never typed397# --------------------------------------------------------------------------------------398def to_num(x):399 """ISQ/StatCan cell -> float or None. '..', 'X', 'x', 'F', '' -> None."""400 if x is None or (isinstance(x, float) and pd.isna(x)):401 return None402 if isinstance(x, (int, float)):403 return float(x)404 s = str(x).strip().replace(" ", "").replace(" ", "")405 if s in {"", "..", "...", "X", "x", "F", "—", "–", "-", "nan", "NaN"}:406 return None407 s = s.replace(",", ".")408 s = re.sub(r"[^\d.\-]", "", s) # strip 'p','r','ᵖ','ʳ' flags409 try:410 return float(s)411 except ValueError:412 return None413414415def parse_isq_population(p: Path) -> tuple[dict, dict, str]:416 df = pd.read_excel(p, header=None)417 hdr = df.iloc[1].tolist()418 # latest column = last header that starts with 4 digits419 cols = [(i, str(h)) for i, h in enumerate(hdr) if re.match(r"^\d{4}", str(h))]420 i_last, label = cols[-1]421 i_2021 = next(i for i, h in cols if str(h).startswith("2021"))422 out, out21 = {}, {}423 for _, row in df.iloc[3:].iterrows():424 code = row[0]425 if isinstance(code, (int, float)) and not pd.isna(code) and 1 <= int(code) <= 17:426 c = f"{int(code):02d}"427 assert CODE_BY_NAME[norm(row[1])] == c, (row[1], c)428 out[c] = to_num(row[i_last]); out21[c] = to_num(row[i_2021])429 return out, out21, label # label like '2025p'430431432def parse_isq_business(p: Path) -> dict:433 df = pd.read_excel(p, sheet_name="Qc RA_industrie", header=None)434 hdr = df.iloc[4].tolist()435 colmap = {}436 for i, h in enumerate(hdr):437 m = re.search(r"\((\d{2})\)\s*$", str(h))438 if m:439 colmap[m.group(1)] = i440 assert len(colmap) == 17, colmap441 rows = {"all": "Ensemble des industries", "mfg": "Fabrication [31-33]", "retail": "Commerce de détail [44-45]"}442 out = {k: {} for k in rows}443 seen = set()444 for _, row in df.iterrows():445 lab = str(row[0]).strip()446 for k, want in rows.items():447 if lab == want and k not in seen: # first occurrence = 2-digit NAICS block448 seen.add(k)449 for c, i in colmap.items():450 out[k][c] = to_num(row[i])451 assert seen == set(rows), seen452 return out453454455def parse_isq_business_size(p: Path) -> dict:456 df = pd.read_excel(p, sheet_name="Qc par RA_taille", header=None)457 out = {}458 for _, row in df.iterrows():459 m = re.search(r"\((\d{2})\)\s*$", str(row[0]))460 if m:461 out[m.group(1)] = to_num(row[1])462 return out463464465def parse_isq_income(p: Path) -> tuple[dict, dict]:466 rows = list(csv.reader(open(p, encoding="utf-8"), delimiter=";"))467 hdr = rows[0]468 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_")]469 best: dict[str, tuple[int, str, float]] = {}470 for r in rows[1:]:471 if r[2].strip() != "Revenu disponible" or r[3].strip() != "$/hab":472 continue473 code = r[0].strip()474 if code not in {c for c, _ in REGIONS}:475 continue476 for it, iv in pairs:477 lab = r[it].strip()478 m = re.match(r"^(\d{4})", lab)479 if not m:480 continue481 v = to_num(r[iv])482 if v is None:483 continue484 y = int(m.group(1))485 if code not in best or y > best[code][0]:486 best[code] = (y, lab, v)487 vals = {c: best[c][2] for c in best}488 labels = {c: best[c][1] for c in best}489 return vals, labels490491492def parse_isq_gdp(p: Path) -> tuple[dict, str]:493 x = pd.ExcelFile(p)494 sheet = x.sheet_names[-1] # most recent year (e.g. '2024e')495 df = x.parse(sheet, header=None)496 hdr = df.iloc[2].tolist()497 colmap = {}498 for i, h in enumerate(hdr[1:], start=1):499 k = norm(h)500 if k in CODE_BY_NAME:501 colmap[CODE_BY_NAME[k]] = i502 assert len(colmap) == 17, colmap503 rows = {"all": "Ensemble des industries (T001)", "mfg": "Fabrication (31-33)", "retail": "Commerce de détail (44-45)",504 "goods": "Secteur de production de biens (T002)", "services": "Secteur des services (T003)"}505 out = {k: {} for k in rows}506 for _, row in df.iterrows():507 lab = str(row[0]).strip()508 for k, want in rows.items():509 if lab == want:510 for c, i in colmap.items():511 out[k][c] = to_num(row[i])512 for k in rows:513 assert out[k], k514 return out, sheet515516517def parse_isq_eqami_ra(p: Path) -> tuple[dict, dict]:518 rows = list(csv.reader(open(p, encoding="utf-8"), delimiter=";"))519 out, combined = {}, {}520 for r in rows[1:]:521 code, name = r[0].strip(), r[1].strip()522 k = norm(name)523 if code == "09-10":524 combined["09"] = to_num(r[2]); combined["10"] = to_num(r[2])525 elif len(code) == 2 and k in CODE_BY_NAME and CODE_BY_NAME[k] == code and code not in out:526 out[code] = to_num(r[2]) # an_2012527 return out, combined528529530def read_statcan_zip(p: Path) -> pd.DataFrame:531 with zipfile.ZipFile(p) as z:532 name = [n for n in z.namelist() if n.endswith(".csv") and "MetaData" not in n][0]533 return pd.read_csv(io.BytesIO(z.read(name)), low_memory=False)534535536def parse_census_er(p: Path) -> dict:537 d = read_statcan_zip(p)538 d.columns = [c.split(": ")[-1].split(" [")[0] if c.startswith("Population and") else c for c in d.columns]539 q = d[d["DGUID"].astype(str).str.startswith("2021S050024")]540 out = {"pop": {}, "area": {}, "dens": {}, "eruid": {}}541 for _, r in q.iterrows():542 c = CODE_BY_NAME[norm(r["GEO"])]543 out["pop"][c] = to_num(r["Population, 2021"])544 out["area"][c] = to_num(r["Land area in square kilometres, 2021"])545 out["dens"][c] = to_num(r["Population density per square kilometre, 2021"])546 out["eruid"][c] = str(r["DGUID"])[-4:]547 assert len(out["pop"]) == 17548 return out549550551def parse_lfs(p464: Path, p466: Path) -> tuple[dict, int]:552 d = read_statcan_zip(p464)553 year = int(d["REF_DATE"].max())554 q = d[(d["REF_DATE"] == year) & d["GEO"].str.endswith(", Quebec")]555 out = {"emp": {}, "urate": {}, "erate": {}, "pop15": {}, "combined": {}}556 for _, r in q.iterrows():557 geo = r["GEO"].replace(", Quebec", "")558 ch = r["Labour force characteristics"]559 key = {"Employment": "emp", "Unemployment rate": "urate", "Employment rate": "erate", "Population": "pop15"}.get(ch)560 if key is None:561 continue562 if geo == "Côte-Nord and Nord-du-Québec":563 out["combined"].setdefault(key, to_num(r["VALUE"]))564 continue565 out[key][CODE_BY_NAME[norm(geo)]] = to_num(r["VALUE"])566 d2 = read_statcan_zip(p466)567 y2 = int(d2["REF_DATE"].max()); assert y2 == year568 q2 = d2[(d2["REF_DATE"] == year) & d2["GEO"].str.endswith(", Quebec")]569 out["mfg"], out["wrt"] = {}, {}570 for _, r in q2.iterrows():571 geo = r["GEO"].replace(", Quebec", "")572 ch = r["Employment characteristics"]573 key = {"Manufacturing [31-33]": "mfg", "Wholesale and retail trade [41, 44-45]": "wrt"}.get(ch)574 if key is None:575 continue576 if geo == "Côte-Nord and Nord-du-Québec":577 out["combined"].setdefault(key, to_num(r["VALUE"]))578 continue579 out[key][CODE_BY_NAME[norm(geo)]] = to_num(r["VALUE"])580 return out, year581582583def fetch_wds_population(mf: Manifest, refresh: bool) -> tuple[dict, str]:584 p = DL / WDS_POP["file"]585 if refresh or not p.exists():586 body = [{"productId": WDS_POP["product_id"], "coordinate": f"{m}.1.1.0.0.0.0.0.0.0", "latestN": 5}587 for m in WDS_POP["members"]]588 r = SESSION.post("https://www150.statcan.gc.ca/t1/wds/rest/getDataFromCubePidCoordAndLatestNPeriods",589 json=body, timeout=120)590 r.raise_for_status()591 meta = SESSION.post("https://www150.statcan.gc.ca/t1/wds/rest/getCubeMetadata",592 json=[{"productId": WDS_POP["product_id"]}], timeout=120).json()[0]["object"]593 members = {m["memberId"]: m["memberNameEn"] for m in meta["dimension"][0]["member"]}594 payload = {"retrieved_utc": utcnow(), "request": body, "geography_members": {str(m): members[m] for m in WDS_POP["members"]},595 "cube_release_time": meta.get("releaseTime"), "data": r.json()}596 p.write_text(json.dumps(payload, ensure_ascii=False, indent=1), encoding="utf-8")597 retrieved = payload["retrieved_utc"]598 else:599 payload = json.loads(p.read_text(encoding="utf-8")); retrieved = payload["retrieved_utc"]600 out, ref = {}, None601 for item in payload["data"]:602 o = item["object"]; member = int(o["coordinate"].split(".")[0])603 geo = payload["geography_members"][str(member)].replace(", Quebec", "")604 if geo == "Quebec":605 continue606 pt = o["vectorDataPoint"][-1]607 out[CODE_BY_NAME[norm(geo)]] = to_num(pt["value"]); ref = pt["refPer"]608 assert len(out) == 17609 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)",610 landing_page=WDS_POP["page"], retrieved_utc=retrieved, licence=WDS_POP["licence"], vintage=WDS_POP["vintage"],611 notes=WDS_POP["notes"] + f" Latest refPer in file: {ref}.", sha256=sha256_path(p), bytes=p.stat().st_size)612 return out, ref613614615# --------------------------------------------------------------------------------------616# geography617# --------------------------------------------------------------------------------------618def rest_layer_geojson(layer: int, fields: str, access_log: dict) -> dict:619 """Pull all Quebec (PRUID='24') features of a layer. The server returns HTTP 500 for a620 where-clause with geometry, so we page by objectIds (batches of 10) and fall back to621 single ids with geometryPrecision=6 for very large polygons (Côte-Nord, Nord-du-Québec...)."""622 ids = http_get(f"{REST_BASE}/{layer}/query", params={"where": "PRUID='24'", "returnIdsOnly": "true", "f": "json"}).json()["objectIds"]623 feats, notes = [], []624 for i in range(0, len(ids), 10):625 chunk = ids[i:i + 10]626 got = None627 for variant in ({}, {"geometryPrecision": 6}):628 try:629 d = http_get(f"{REST_BASE}/{layer}/query", params={"objectIds": ",".join(map(str, chunk)), "outFields": fields,630 "returnGeometry": "true", "f": "geojson", **variant}, retries=2).json()631 if "features" in d:632 got = d["features"]; break633 except Exception: # noqa: BLE001634 continue635 if got is None: # single ids with precision 6636 got = []637 for oid in chunk:638 d = http_get(f"{REST_BASE}/{layer}/query", params={"objectIds": oid, "outFields": fields, "returnGeometry": "true",639 "f": "geojson", "geometryPrecision": 6}).json()640 got += d["features"]; notes.append(f"objectId {oid} fetched singly with geometryPrecision=6")641 feats += got642 assert len(feats) == len(ids), (len(feats), len(ids))643 access_log[f"rest_layer_{layer}"] = {"url": f"{REST_BASE}/{layer}/query", "features": len(feats), "notes": notes, "retrieved_utc": utcnow()}644 return {"type": "FeatureCollection", "features": feats}645646647def content_hash(gdf, id_field: str) -> str:648 h = hashlib.sha256()649 for _, r in gdf.sort_values(id_field).iterrows():650 h.update(str(r[id_field]).encode()); h.update(r.geometry.wkb)651 return h.hexdigest()652653654def build_geo(mf: Manifest, access_log: dict, refresh: bool, no_geo: bool, keep_archives: bool) -> dict:655 import geopandas as gpd # heavy import kept local656657 GEO.mkdir(parents=True, exist_ok=True)658 res: dict = {}659 # ---- MERN administrative regions (CC-BY 4.0) ----660 spec = SOURCES["mern_bdat"]661 zp = DL / spec["file"]662 ra_gpkg = GEO / "qc_regions_admin_mern_bdat100k.gpkg"663 if refresh or not ra_gpkg.exists() or zp.exists():664 if refresh or not zp.exists():665 fetch_source("mern_bdat", spec, mf, refresh=True)666 else:667 fetch_source("mern_bdat", spec, mf, refresh=False)668 tmp = EXT / "_tmp_bdat"669 with zipfile.ZipFile(zp) as z:670 for n in z.namelist():671 base = n.split("/")[-1]672 if base.startswith(("regio_s.", "mrc_s.", "munic_s.")) and not base.endswith((".lock", ".xml", ".qmd")):673 (tmp / base).parent.mkdir(parents=True, exist_ok=True)674 (tmp / base).write_bytes(z.read(n))675 regio = gpd.read_file(tmp / "regio_s.shp")676 mrc = gpd.read_file(tmp / "mrc_s.shp")677 munic = gpd.read_file(tmp / "munic_s.shp")678 reg = regio.dissolve(by="RES_CO_REG", as_index=False)[["RES_CO_REG", "RES_NM_REG", "RES_CO_REF", "RES_CO_VER", "geometry"]]679 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"})680 assert len(reg) == 17 and set(reg.region_code) == {c for c, _ in REGIONS}681 reg["region_name"] = reg.region_code.map(dict(REGIONS))682 if ra_gpkg.exists():683 ra_gpkg.unlink()684 reg.to_file(ra_gpkg, layer="regions_admin", driver="GPKG")685 regio.to_file(ra_gpkg, layer="regio_s_raw", driver="GPKG")686 mrc_tab = mrc[["MRS_CO_MRC", "MRS_NM_MRC", "MRS_CO_REG", "MRS_NM_REG", "MRS_DE_IND", "geometry"]]687 mrc_tab.to_file(ra_gpkg, layer="mrc_s_raw", driver="GPKG")688 mun_cols = [c for c in munic.columns if c.startswith("MUS_CO_") or c.startswith("MUS_NM_") or c == "MUS_DE_IND"]689 munic[mun_cols].to_csv(EXT / "mern_municipalites_attributes.csv", index=False, encoding="utf-8")690 for f in tmp.iterdir():691 f.unlink()692 tmp.rmdir()693 side = DL / (spec["file"] + ".sha256")694 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")695 if spec.get("delete_after_extract") and not keep_archives:696 zp.unlink(); log("deleted archive mern_bdat_adm_100k_shp.zip after extraction (sha256 kept in manifest + .sha256 sidecar)")697 else:698 side = DL / (spec["file"] + ".sha256")699 prev = mf.old.get(f"dl/{spec['file']}") or (json.loads(side.read_text(encoding="utf-8")) if side.exists() else None)700 if prev:701 mf.add(f"dl/{spec['file']}", source_name=spec["source_name"], url=spec["url"], landing_page=spec["page"],702 retrieved_utc=prev.get("retrieved_utc"), licence=spec["licence"], vintage=spec["vintage"],703 notes=spec["notes"] + " (archive deleted after extraction; checksum from .sha256 sidecar)",704 sha256=prev["sha256"], bytes=prev.get("bytes"))705 reg = gpd.read_file(ra_gpkg, layer="regions_admin")706 mrc = gpd.read_file(ra_gpkg, layer="mrc_s_raw")707 mf.add(f"geo/{ra_gpkg.name}", source_name=spec["source_name"] + " -> layers regions_admin (dissolved, 17), regio_s_raw, mrc_s_raw",708 url=spec["url"], landing_page=spec["page"], retrieved_utc=utcnow(), licence=spec["licence"], vintage=spec["vintage"],709 notes="Derived by src/ingest/external.py (dissolve by RES_CO_REG). CRS EPSG:4269.", sha256=sha256_path(ra_gpkg),710 content_sha256=content_hash(reg, "region_code"), bytes=ra_gpkg.stat().st_size, features=len(reg))711 reg_lam = reg.to_crs(6622) # NAD83(CSRS) / Quebec Lambert, metres712 res["area_mern_km2"] = dict(zip(reg_lam.region_code, (reg_lam.geometry.area / 1e6).round(2)))713 res["mern_version"] = reg.mern_version.iloc[0]714 # municipalities attribute table: municipality code -> region (for CSD-level concordance)715 res["mern_mrc_regions"] = mrc.groupby("MRS_CO_MRC")["MRS_CO_REG"].agg(lambda s: sorted(set(s))).to_dict()716717 # ---- StatCan 2021 layers (Open Licence) ----718 if no_geo:719 log("--no-geo: skipping StatCan REST pulls; using existing GeoPackages if present")720 layers = {}721 rest_notes = ("Quebec subset (PRUID=24). Direct zip on www12.statcan.gc.ca returns HTTP 403 to non-browser clients "722 "(see access_log.json); the REST service is the official alternative listed on open.canada.ca. "723 "CRS EPSG:4326 as returned by f=geojson. LANDAREA attribute is StatCan's land area in km2. "724 "Very large polygons (e.g. Côte-Nord, Nord-du-Québec, FSA G0G) were fetched singly with "725 "geometryPrecision=6 (~0.1 m) because the server returns HTTP 500 at full precision.")726 for key, L in REST_LAYERS.items():727 gp = GEO / (L["gpkg"] or f"qc_{key}_2021_statcan_cbf_points.gpkg")728 kind = "polygons" if L["gpkg"] else "representative points only (polygons not kept to save space)"729 if gp.exists() and not refresh:730 layers[key] = gpd.read_file(gp)731 prev = mf.old.get(f"geo/{gp.name}", {})732 retrieved = prev.get("retrieved_utc") or datetime.fromtimestamp(gp.stat().st_mtime, timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")733 log(f"cached geo layer {key}: {gp.name}")734 elif no_geo:735 continue736 else:737 gj = rest_layer_geojson(L["layer"], L["fields"], access_log)738 g = gpd.GeoDataFrame.from_features(gj["features"], crs="EPSG:4326")739 if not L["gpkg"]:740 g["geometry"] = g.to_crs(6622).representative_point().to_crs(4326)741 if gp.exists():742 gp.unlink()743 g.to_file(gp, layer=key, driver="GPKG")744 layers[key] = g; retrieved = utcnow()745 mf.add(f"geo/{gp.name}", source_name=L["source_name"] + f" ({kind})",746 url=f"{REST_BASE}/{L['layer']}/query?where=PRUID='24' (paged by objectIds, f=geojson)",747 landing_page="https://open.canada.ca/data/dataset/ef70dc3b-1069-4037-9bce-61f47e628a1d",748 retrieved_utc=retrieved, licence=LIC_STATCAN, vintage="2021 Census cartographic boundary files (l***000b21a_e), release 2022",749 notes=rest_notes, sha256=sha256_path(gp), content_sha256=content_hash(layers[key], L["id_field"]),750 bytes=gp.stat().st_size, features=len(layers[key]))751 if "er" not in layers or "fsa" not in layers or "cd" not in layers:752 log("geo layers incomplete (run without --no-geo once); skipping concordance products")753 return res754755 er, cd, fsa = layers["er"], layers["cd"], layers["fsa"]756 assert len(er) == 17, len(er)757 er["region_code"] = er.ERNAME.map(lambda n: CODE_BY_NAME[norm(n)])758 res["er_names"] = dict(zip(er.region_code, er.ERNAME)); res["eruid"] = dict(zip(er.region_code, er.ERUID))759 res["er_landarea_km2"] = dict(zip(er.region_code, er.LANDAREA))760761 # ---- RA (MERN) vs ER (StatCan): measured overlap, not assumed ----762 # make_valid: both sources contain self-touching rings (GEOS 'side location conflict'); areas are unaffected.763 er["geometry"] = er.geometry.make_valid(); reg["geometry"] = reg.geometry.make_valid(); fsa["geometry"] = fsa.geometry.make_valid()764 er_l, reg_l = er.to_crs(6622), reg.to_crs(6622)765 rows = []766 for _, r in reg_l.iterrows():767 e = er_l[er_l.region_code == r.region_code].iloc[0]768 inter = r.geometry.intersection(e.geometry).area; union = r.geometry.union(e.geometry).area769 rows.append(dict(region_code=r.region_code, region_name=r.region_name, eruid=e.ERUID, er_name=e.ERNAME,770 area_mern_km2=round(r.geometry.area / 1e6, 1), area_er_km2=round(e.geometry.area / 1e6, 1),771 intersection_km2=round(inter / 1e6, 1), iou=round(inter / union, 4),772 share_of_mern_in_er=round(inter / r.geometry.area, 4),773 share_of_er_in_mern=round(inter / e.geometry.area, 4)))774 ov = pd.DataFrame(rows).sort_values("region_code")775 ov.to_csv(EXT / "ra_er_overlap_2021.csv", index=False, encoding="utf-8")776 res["iou"] = dict(zip(ov.region_code, ov.iou)); res["er_in_mern"] = dict(zip(ov.region_code, ov.share_of_er_in_mern))777778 # ---- CD -> ER (point in StatCan ER) and CD -> RA (point in MERN RA) ----779 cd_pts = cd.copy()780 if (cd_pts.geom_type != "Point").any(): # fresh polygons -> representative points781 cd_pts["geometry"] = cd_pts.to_crs(6622).representative_point().to_crs(4326)782 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")783 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")784 j2["mrc_code_mern"] = j2.CDUID.str[2:]785 j2["same_region"] = j2.er_region_code == j2.mern_region_code786 j2["lon"] = j2.geometry.x.round(5); j2["lat"] = j2.geometry.y.round(5)787 conc = j2.drop(columns="geometry").sort_values("CDUID")788 conc.to_csv(EXT / "cd_to_region_concordance_2021.csv", index=False, encoding="utf-8")789 res["cd_mismatch"] = conc[~conc.same_region][["CDUID", "CDNAME", "er_region_code", "mern_region_code"]].to_dict("records")790791 # ---- FSA -> region: centroid + dominant-area assignment ----792 f_l = fsa.to_crs(6622); reg_l = reg.to_crs(6622); er_l = er.to_crs(6622)793 cen = f_l.copy(); cen["geometry"] = f_l.representative_point()794 a = gpd.sjoin(cen[["CFSAUID", "LANDAREA", "geometry"]], reg_l[["region_code", "region_name", "geometry"]], how="left", predicate="within").drop(columns="index_right")795 a = a.rename(columns={"region_code": "region_code_centroid", "region_name": "region_name_centroid"})796 b = gpd.sjoin(cen[["CFSAUID", "geometry"]], er_l[["region_code", "ERUID", "geometry"]], how="left", predicate="within").drop(columns="index_right")797 a["er_region_code_centroid"] = a.CFSAUID.map(dict(zip(b.CFSAUID, b.region_code)))798 a["eruid_centroid"] = a.CFSAUID.map(dict(zip(b.CFSAUID, b.ERUID)))799 # dominant area share800 ov2 = gpd.overlay(f_l[["CFSAUID", "geometry"]], reg_l[["region_code", "geometry"]], how="intersection", keep_geom_type=False)801 ov2["a"] = ov2.geometry.area802 tot = ov2.groupby("CFSAUID")["a"].sum()803 ov2["share"] = ov2.a / ov2.CFSAUID.map(tot)804 dom = ov2.sort_values("share", ascending=False).drop_duplicates("CFSAUID").set_index("CFSAUID")805 a["region_code_dominant_area"] = a.CFSAUID.map(dom.region_code)806 a["dominant_area_share"] = a.CFSAUID.map(dom.share).round(4)807 a["n_regions_touched"] = a.CFSAUID.map(ov2[ov2.share > 0.001].groupby("CFSAUID").size())808 pt = cen.to_crs(4326)809 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))))810 a = a.drop(columns="geometry").sort_values("CFSAUID")811 a.to_csv(EXT / "fsa_to_region_2021.csv", index=False, encoding="utf-8")812 res["fsa_n"] = len(a); res["fsa_multi"] = int((a.n_regions_touched > 1).sum())813 res["fsa_centroid_vs_dominant_differ"] = int((a.region_code_centroid != a.region_code_dominant_area).sum())814 res["fsa_count_by_region_centroid"] = a.groupby("region_code_centroid").size().to_dict()815 for name, df in (("ra_er_overlap_2021.csv", ov), ("cd_to_region_concordance_2021.csv", conc), ("fsa_to_region_2021.csv", a)):816 p = EXT / name817 mf.add(name, source_name="Derived by src/ingest/external.py from the MERN and StatCan layers above", url="(derived)",818 retrieved_utc=utcnow(), licence="Derived product; inherits CC-BY 4.0 (MRNF) and StatCan Open Licence",819 vintage="MERN V2026-10 x StatCan 2021 boundaries", notes="See README.md", sha256=sha256_path(p), bytes=p.stat().st_size, rows=len(df))820 return res821822823# --------------------------------------------------------------------------------------824# access probes (recorded, so the README claims are reproducible)825# --------------------------------------------------------------------------------------826def probe(url: str, access_log: dict, key: str) -> None:827 try:828 r = SESSION.head(url, timeout=60, allow_redirects=True)829 access_log[key] = {"url": url, "http_status": r.status_code, "probed_utc": utcnow()}830 except requests.RequestException as e:831 access_log[key] = {"url": url, "http_status": None, "error": str(e), "probed_utc": utcnow()}832833834# --------------------------------------------------------------------------------------835# main836# --------------------------------------------------------------------------------------837def main() -> int:838 ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)839 ap.add_argument("--refresh", action="store_true", help="force re-download of every source")840 ap.add_argument("--no-geo", action="store_true", help="skip StatCan REST geometry pulls (reuse existing GeoPackages)")841 ap.add_argument("--strict", action="store_true", help="abort if a re-downloaded file's sha256 differs from manifest")842 ap.add_argument("--keep-archives", action="store_true", help="keep the 87 MB MERN zip after extraction")843 args = ap.parse_args()844845 EXT.mkdir(parents=True, exist_ok=True); DL.mkdir(exist_ok=True)846 mf = Manifest(load_manifest(), strict=args.strict)847 access_log = json.loads(ACCESS_LOG.read_text(encoding="utf-8")) if ACCESS_LOG.exists() else {}848849 # 1) tabular sources850 paths = {k: fetch_source(k, s, mf, args.refresh) for k, s in SOURCES.items() if k != "mern_bdat"}851 pop_isq, pop_isq_2021, pop_label = parse_isq_population(paths["isq_pop"])852 biz = parse_isq_business(paths["isq_biz_ind"])853 biz_size_total = parse_isq_business_size(paths["isq_biz_size"])854 assert all(abs(biz["all"][c] - biz_size_total[c]) < 0.5 for c in biz["all"]), "ISQ business totals differ between files"855 inc, inc_lab = parse_isq_income(paths["isq_income"])856 gdp, gdp_sheet = parse_isq_gdp(paths["isq_gdp"])857 inet, inet_comb = parse_isq_eqami_ra(paths["isq_eqami2012_ra"])858 cen = parse_census_er(paths["sc_98100009"])859 lfs, lfs_year = parse_lfs(paths["sc_14100464"], paths["sc_14100466"])860 pop_sc, pop_sc_ref = fetch_wds_population(mf, args.refresh)861862 # 2) access probes863 probe(WWW12_FSA_ZIP, access_log, "www12_fsa_zip")864 probe(WWW12_ER_DEF, access_log, "www12_census_dictionary_er_definition")865 probe(SGC2016_INTRO, access_log, "statcan_sgc2016_introduction")866 probe(SOURCES["mern_bdat"]["url"], access_log, "mern_bdat_zip")867868 # 3) geography869 geo = build_geo(mf, access_log, args.refresh, args.no_geo, args.keep_archives)870871 # 4) assemble tidy table872 pop_year = re.match(r"\d{4}", pop_label).group(0)873 rows = []874 for code, name in REGIONS:875 r = {"region_code": code, "region_name": name}876877 def put(col, val, vintage, source):878 if isinstance(val, float) and val.is_integer():879 val = int(val)880 r[col] = "" if val is None else val881 r[f"{col}_vintage"] = vintage if val is not None else ""882 r[f"{col}_source"] = source if val is not None else ""883884 put("population", pop_isq.get(code), f"1 July {pop_year} ({pop_label}: provisional)", "ISQ table 3595 (StatCan annual demographic estimates adapted by ISQ)")885 put("population_statcan_er", pop_sc.get(code), f"{pop_sc_ref[:4]}-07-01 (2021 ER boundaries)", "StatCan 17-10-0150-01 (WDS)")886 put("population_census_2021_er", cen["pop"].get(code), "Census 2021 (2021 ER boundaries)", "StatCan 98-10-0009-01")887 put("business_locations_with_employees", biz["all"].get(code), "Dec 2025", "ISQ (StatCan 33-10-1097-01 aggregated to RA)")888 put("business_locations_with_employees_mfg_31_33", biz["mfg"].get(code), "Dec 2025", "ISQ (StatCan 33-10-1097-01 aggregated to RA)")889 put("business_locations_with_employees_retail_44_45", biz["retail"].get(code), "Dec 2025", "ISQ (StatCan 33-10-1097-01 aggregated to RA)")890 put("disposable_income_per_capita_cad", inc.get(code), inc_lab.get(code, "").replace("ᵖ", "p").replace("ʳ", "r"), "ISQ table 3962 (comptes économiques régionaux)")891 put("gdp_basic_prices_kcad", gdp["all"].get(code), gdp_sheet, "ISQ table 2299 (PIB régional)")892 put("gdp_mfg_31_33_kcad", gdp["mfg"].get(code), gdp_sheet, "ISQ table 2299 (PIB régional)")893 put("gdp_retail_44_45_kcad", gdp["retail"].get(code), gdp_sheet, "ISQ table 2299 (PIB régional)")894 put("employment_lfs_thousands", lfs["emp"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0464-01")895 put("employment_lfs_mfg_thousands", lfs["mfg"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0466-01")896 put("employment_lfs_wholesale_retail_thousands", lfs["wrt"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0466-01")897 put("unemployment_rate_pct", lfs["urate"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0464-01")898 put("employment_rate_pct", lfs["erate"].get(code), f"{lfs_year} annual average", "StatCan 14-10-0464-01")899 comb = lfs["combined"] if code in ("09", "10") else {}900 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")901 put("unemployment_rate_pct_combined_09_10", comb.get("urate"), f"{lfs_year} annual average (combined ER)", "StatCan 14-10-0464-01")902 put("land_area_km2", cen["area"].get(code), "Census 2021 (2021 ER boundaries)", "StatCan 98-10-0009-01")903 put("population_density_per_km2", cen["dens"].get(code), "Census 2021 (2021 ER boundaries)", "StatCan 98-10-0009-01")904 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)")905 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)")906 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")907 r["statcan_eruid_2021"] = cen["eruid"].get(code, "")908 r["ra_er_boundary_iou_2021"] = geo.get("iou", {}).get(code, "")909 r["share_of_statcan_er_land_inside_mern_ra"] = geo.get("er_in_mern", {}).get(code, "")910 r["n_fsa_centroids_2021"] = geo.get("fsa_count_by_region_centroid", {}).get(code, "")911 rows.append(r)912 out = pd.DataFrame(rows)913 out.to_csv(OUT_CSV, index=False, encoding="utf-8")914 mf.add("regions_external.csv", source_name="Assembled by src/ingest/external.py from the files above", url="(derived)",915 retrieved_utc=utcnow(), licence="Derived; see each *_source column", vintage="see *_vintage columns",916 notes="Empty cell = not obtained (never zero).", sha256=sha256_path(OUT_CSV), bytes=OUT_CSV.stat().st_size, rows=len(out))917 mf.write(access_log)918919 # 5) console summary920 log(f"wrote {OUT_CSV.relative_to(PROJECT)} ({len(out)} rows, {len(out.columns)} columns)")921 for col in [c for c in out.columns if not c.endswith(("_vintage", "_source")) and c not in ("region_code", "region_name")]:922 n = int((out[col].astype(str) != "").sum())923 log(f" {col:52s} filled {n:2d}/17")924 if geo.get("cd_mismatch") is not None:925 log(f"CD->region mismatches between StatCan 2021 ER and MERN RA: {geo['cd_mismatch']}")926 if "fsa_n" in geo:927 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']}")928 total = sum(p.stat().st_size for p in EXT.rglob("*") if p.is_file())929 log(f"data/external total size: {total / 1e6:.1f} MB")930 if mf.warnings:931 log(f"{len(mf.warnings)} checksum warning(s) - see manifest.json")932 return 0933934935if __name__ == "__main__":936 sys.exit(main())937