TeX 55.1%
Python 44.2%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Phase 1 — geographic concordance.45Spatially joins every transaction (lat/lng) to the MRNF SDA V2026-07 boundaries6(municipality → MRC → administrative region; metropolitan community; borough),7links the MAMH municipal register (population, designation, constitution date)8and the roll's municipality code, and builds fine-location keys9(assessment neighbourhood unit when present, else 1-km / 500-m metric grid cells).10Writes data/interim/geo_transactions.parquet and results/geo/*.csv.11"""12from __future__ import annotations1314import sys15from pathlib import Path1617sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))18import geopandas as gpd19import numpy as np20import pandas as pd21import polars as pl2223from sharp_hpi import config as C2425OUT = C.RESULTS / "geo"26OUT.mkdir(parents=True, exist_ok=True)272829def main() -> None:30 tx = pl.read_parquet(C.RAW_TX, columns=["id", "date", "lat", "lng", "city", "id_provinc"])31 um = pl.read_parquet(C.RAW_UNITS_MATCHED, columns=["id_provinc", "code_mun", "municipalite", "unite_voisinage", "arrond"])32 tx = tx.join(um, on="id_provinc", how="left")33 pts = gpd.GeoDataFrame(tx.select(["id", "lat", "lng"]).to_pandas(), geometry=gpd.points_from_xy(tx["lng"], tx["lat"]), crs="EPSG:4269")3435 mun = gpd.read_file(C.SDA_GPKG, layer="munic_s")[["MUS_CO_GEO", "MUS_NM_MUN", "MUS_CO_DES", "MUS_DE_IND", "MUS_CO_MRC", "MUS_NM_MRC", "MUS_CO_REG", "MUS_NM_REG", "MUS_DA_CON", "MUS_CO_MOD", "MUS_DA_MOD", "geometry"]]36 mun = mun.to_crs("EPSG:4269")37 j = gpd.sjoin(pts, mun, how="left", predicate="within")38 j = j[~j.index.duplicated(keep="first")]39 j["join_method"] = np.where(j["MUS_CO_GEO"].notna(), "within", "unassigned")40 # nearest within 1 km for the unassigned (coastline / water slivers)41 miss = j["MUS_CO_GEO"].isna()42 if miss.any():43 p_m = pts.loc[miss].to_crs(C.CRS_METRIC)44 near = gpd.sjoin_nearest(p_m, mun.to_crs(C.CRS_METRIC), how="left", max_distance=1000, distance_col="dist_m")45 near = near[~near.index.duplicated(keep="first")]46 for c in ["MUS_CO_GEO", "MUS_NM_MUN", "MUS_CO_DES", "MUS_DE_IND", "MUS_CO_MRC", "MUS_NM_MRC", "MUS_CO_REG", "MUS_NM_REG", "MUS_DA_CON", "MUS_CO_MOD", "MUS_DA_MOD"]:47 j.loc[miss, c] = near[c].values48 j.loc[miss, "join_method"] = np.where(near["MUS_CO_GEO"].notna(), "nearest<=1km", "unassigned")4950 comet = gpd.read_file(C.SDA_GPKG, layer="comet_s")[["CMS_CO_COM", "CMS_NM_COM", "geometry"]].to_crs("EPSG:4269")51 jc = gpd.sjoin(pts, comet, how="left", predicate="within")52 jc = jc[~jc.index.duplicated(keep="first")]53 arr = gpd.read_file(C.SDA_GPKG, layer="arron_s")[["ARS_CO_ARR", "ARS_NM_ARR", "ARS_NM_MUN", "geometry"]].to_crs("EPSG:4269")54 ja = gpd.sjoin(pts, arr, how="left", predicate="within")55 ja = ja[~ja.index.duplicated(keep="first")]5657 # metric grid cells for fine location58 pm = pts.to_crs(C.CRS_METRIC)59 x, y = pm.geometry.x.values, pm.geometry.y.values60 grid1k = (np.floor(x / 1000).astype("int64").astype(str) + "_" + np.floor(y / 1000).astype("int64").astype(str))61 grid500 = (np.floor(x / 500).astype("int64").astype(str) + "_" + np.floor(y / 500).astype("int64").astype(str))62 grid2k = (np.floor(x / 2000).astype("int64").astype(str) + "_" + np.floor(y / 2000).astype("int64").astype(str))63 grid5k = (np.floor(x / 5000).astype("int64").astype(str) + "_" + np.floor(y / 5000).astype("int64").astype(str))6465 geo = pl.from_pandas(pd.DataFrame({66 "id": j["id"].values,67 "mun_code": j["MUS_CO_GEO"].values, "mun_name": j["MUS_NM_MUN"].values, "mun_design": j["MUS_CO_DES"].values,68 "mrc_code": j["MUS_CO_MRC"].values, "mrc_name": j["MUS_NM_MRC"].values,69 "region_code": j["MUS_CO_REG"].values, "region_name": j["MUS_NM_REG"].values,70 "mun_date_constituted": pd.to_datetime(j["MUS_DA_CON"]).dt.strftime("%Y-%m-%d").values,71 "mun_date_modified": pd.to_datetime(j["MUS_DA_MOD"]).dt.strftime("%Y-%m-%d").values,72 "join_method": j["join_method"].values,73 "comet_code": jc["CMS_CO_COM"].values, "comet_name": jc["CMS_NM_COM"].values,74 "borough_code": ja["ARS_CO_ARR"].values, "borough_name": ja["ARS_NM_ARR"].values,75 "x_m": x, "y_m": y, "grid500": grid500, "grid1k": grid1k, "grid2k": grid2k, "grid5k": grid5k,76 }).astype({c: object for c in ["mun_code","mun_name","mun_design","mrc_code","mrc_name","region_code","region_name","mun_date_constituted","mun_date_modified","join_method","comet_code","comet_name","borough_code","borough_name","grid500","grid1k","grid2k","grid5k"]}))77 geo = tx.select(["id", "city", "code_mun", "municipalite", "unite_voisinage", "arrond"]).join(geo, on="id", how="left")78 # MAMH register: population, designation, constitution date, official MRC/region labels79 mamh = pl.read_csv(C.MAMH_MUN_CSV, infer_schema_length=0, encoding="utf8-lossy").select(80 pl.col("mcode").alias("mun_code"), pl.col("munnom").alias("mamh_name"), pl.col("mdes").alias("mamh_design"),81 pl.col("regadm").alias("mamh_region"), pl.col("mrc").alias("mamh_mrc"), pl.col("mpopul").cast(pl.Int64, strict=False).alias("mamh_population"),82 pl.col("mdatcons").alias("mamh_date_constituted"), pl.col("msuperf").cast(pl.Float64, strict=False).alias("mamh_area_km2"),83 )84 geo = geo.join(mamh, on="mun_code", how="left")85 # fine-location keys. The roll's `unite_voisinage` turned out to be a unit-level code (379k distinct values for86 # 449k rows), so fine location is a nested metric grid (500 m ⊂ 1 km ⊂ 2 km ⊂ 5 km ⊂ municipality); the87 # estimation code picks the finest level with enough sales (hedonic.assign_loc_fe).88 geo = geo.with_columns(89 pl.concat_str([pl.lit("G1K:"), pl.col("grid1k")]).alias("loc_fine"),90 pl.concat_str([pl.lit("G1K:"), pl.col("grid1k")]).alias("loc_grid1k"),91 pl.concat_str([pl.lit("G500:"), pl.col("grid500")]).alias("loc_grid500"),92 (pl.col("code_mun") == pl.col("mun_code")).alias("roll_mun_agrees"),93 )94 geo.write_parquet(C.GEO_TX)9596 # --- diagnostics ---97 geo.group_by("join_method").agg(pl.len().alias("n")).sort("n", descending=True).write_csv(OUT / "join_methods.csv")98 agree = geo.filter(pl.col("code_mun").is_not_null()).select(99 pl.len().alias("n_with_roll_code"), pl.col("roll_mun_agrees").mean().alias("share_roll_code_agrees"),100 )101 agree.write_csv(OUT / "roll_vs_spatial_agreement.csv")102 geo.filter(pl.col("code_mun").is_not_null() & ~pl.col("roll_mun_agrees")).group_by(["code_mun", "municipalite", "mun_code", "mun_name"]).agg(pl.len().alias("n")).sort("n", descending=True).head(50).write_csv(OUT / "roll_vs_spatial_disagreements_top50.csv")103 geo.group_by(["region_code", "region_name"]).agg(pl.len().alias("n"), pl.col("mun_code").n_unique().alias("municipalities"), pl.col("mrc_code").n_unique().alias("mrcs")).sort("region_code").write_csv(OUT / "region_summary.csv")104 geo.group_by(["comet_code", "comet_name"]).agg(pl.len().alias("n")).sort("n", descending=True).write_csv(OUT / "metropolitan_communities.csv")105 # boundary changes since 2021 (municipalities constituted/modified during the sample)106 mun_df = pl.from_pandas(mun.drop(columns="geometry")).with_columns(pl.col("MUS_DA_CON").cast(pl.Utf8), pl.col("MUS_DA_MOD").cast(pl.Utf8))107 mun_df.filter((pl.col("MUS_DA_CON") >= "2021-01-01") | (pl.col("MUS_DA_MOD") >= "2021-01-01")).write_csv(OUT / "boundary_changes_since_2021.csv")108 mamh.filter(pl.col("mamh_date_constituted") >= "2021-01-01").write_csv(OUT / "mamh_constituted_since_2021.csv")109 pl.DataFrame({110 "metric": ["transactions", "share_rows_with_roll_neighbourhood_code", "distinct_grid1k", "distinct_grid2k", "distinct_grid500", "distinct_municipalities", "distinct_mrc"],111 "value": [float(geo.height), float((geo["unite_voisinage"].is_not_null() & (geo["unite_voisinage"] != "")).mean()), float(geo["grid1k"].n_unique()), float(geo["grid2k"].n_unique()), float(geo["grid500"].n_unique()), float(geo["mun_code"].n_unique()), float(geo["mrc_code"].n_unique())],112 }).write_csv(OUT / "fine_location_summary.csv")113 print(geo.group_by("join_method").agg(pl.len()).to_dicts(), agree.to_dicts())114115116if __name__ == "__main__":117 main()118