1"""Phase 2.3 — Store master with region recovery and geographic-quality flags.23Input : data/interim/stores_raw.parquet, data/interim/listings_clean.parquet, optional data/external/<regions layer>4Output: data/interim/stores_master.parquet, experiments/store_master.json56Region logic (documented in reports/phase2_data_audit.md):7 region_db : Fabri-Ka's region (directory source / site mentions / FSA map), blank for 22 % of live stores8 region_pip : administrative region containing the FSA centroid (lat, lng) — point-in-polygon against the official9 MERN layer (if available; precision = neighbourhood, never the civic address)10 region_final : region_db if present, else region_pip, else 'Inconnue'11 region_source : db | pip | unknown12 region_conflict: region_db and region_pip both present and different (kept, flagged; sensitivity analyses drop them)13 phone_shared : phone number shared by >= 3 live stores (directory contamination; phone never used for geography)14"""15from __future__ import annotations16import glob, json17import numpy as np, pandas as pd18from src.common import ROOT, P, REGIONS, get_logger, ledger, register_output, write_json, BASE_PROC, BASE_EXP1920log = get_logger("clean.stores")212223def load_regions_layer():24 cands = sorted(glob.glob(str(ROOT / "data" / "external" / "**" / "*mern*.gpkg"), recursive=True)) # official MERN layer only (no silent fallback)25 try:26 import geopandas as gpd27 except ImportError:28 return None, None29 for c in cands:30 try:31 g = gpd.read_file(c)32 except Exception:33 continue34 # find a name column whose values cover the 17 regions35 for col in g.columns:36 if g[col].dtype == object:37 vals = set(str(v).replace("--", "–") for v in g[col].dropna().unique())38 hit = len([r for r in REGIONS if r in vals or r.replace("–", "-") in vals])39 if hit >= 15:40 g["geometry"] = g.geometry.make_valid(); g = g.to_crs(4326)41 g["_region"] = g[col].astype(str).str.replace("--", "–").map(lambda v: next((r for r in REGIONS if r == v or r.replace("–", "-") == v), v))42 return g[["_region", "geometry"]].dissolve(by="_region").reset_index(), c43 return None, None444546def city_region_lookup():47 """Map the zippopotam sector label (e.g. 'Granby Central', 'Estrie-Ouest (Fulford)') to an administrative region via the48 official MERN municipality table; only unambiguous municipality names are used."""49 import unicodedata, re50 f = ROOT / "data" / "external" / "mern_municipalites_attributes.csv"51 if not f.exists():52 return lambda c: None53 mu = pd.read_csv(f)54 norm = lambda x: re.sub(r"[^a-z]", "", unicodedata.normalize("NFKD", str(x)).encode("ascii", "ignore").decode().lower())55 mu["k"] = mu.MUS_NM_MUN.map(norm)56 amb = mu.groupby("k").MUS_NM_REG.nunique(); table = mu.drop_duplicates("k").set_index("k").MUS_NM_REG; table = table[amb[table.index] == 1]57 def look(c):58 if not isinstance(c, str) or not c: return None59 base = re.sub(r"\(.*?\)", "", c); base = re.sub(r"\b(Central|West|East|North|South|Northwest|Northeast|Southwest|Southeast|Centre|Ouest|Est|Nord|Sud)\b", "", base).strip(" -")60 k = norm(base) # the parenthesised town of rural FSAs is only a representative place, never used61 return table[k] if k and k in table.index else None62 return look636465def main():66 st = pd.read_parquet(P["interim"] / "stores_raw.parquet")67 li = pd.read_parquet(P["interim"] / "listings_clean.parquet")68 live = st[(st.enabled == 1) & (st.product_count > 0)].copy()69 live["region_db"] = live.region.fillna("").replace({"": None})70 live["region_db"] = live.region_db.where(live.region_db.isin(REGIONS), None)7172 layer, layer_path = load_regions_layer()73 live["region_pip"] = None74 if layer is not None:75 import geopandas as gpd76 pts = live[live.lat.notna() & live.lng.notna()]77 gp = gpd.GeoDataFrame(pts[["id"]], geometry=gpd.points_from_xy(pts.lng, pts.lat), crs=4326)78 j = gpd.sjoin(gp, layer, how="left", predicate="within")79 j = j.drop_duplicates("id")80 live = live.merge(j[["id", "_region"]], on="id", how="left")81 live["region_pip"] = live["_region"]; live = live.drop(columns=["_region"])82 log.info("regions layer %s — PIP assigned for %s stores", layer_path, int(live.region_pip.notna().sum()))83 else:84 log.warning("no administrative-region layer found under data/external — region_pip left empty")8586 look = city_region_lookup(); urban = live.postal_prefix.fillna("").str.len().eq(3) & ~live.postal_prefix.fillna("").str[1].eq("0")87 live["region_city"] = live.city.where(urban).map(look) # municipality lookup only for urban FSAs (second character != 0)88 live["region_city"] = live.region_city.where(live.region_city.isin(REGIONS), None)89 # Amendment A-001 (2026-11-01): official administrative geography first — municipality lookup, then FSA-centroid90 # point-in-polygon, then Fabri-Ka's directory/site label (directory regions are often tourism/brand regions, e.g.91 # 'Cantons-de-l'Est' members in Montérégie labelled Estrie). region_db kept for the db-first sensitivity variant.92 # A-001 (final form): municipality lookup (urban FSAs, exact) > Fabri-Ka label > FSA-centroid point-in-polygon93 # (centroids of rural FSAs are unreliable near region borders); pip-first kept as a sensitivity variant.94 # MRC of the FSA centroid (MERN mrc layer) to handle the 28 July 2021 transfer of Brome-Missisquoi and La Haute-Yamaska95 # from Montérégie to Estrie: platform labels predating the transfer are overridden (referee 2, must-fix)96 live["mrc_pip"] = None97 if layer_path:98 try:99 import geopandas as gpd100 mrc = gpd.read_file(layer_path, layer="mrc_s_raw"); mrc["geometry"] = mrc.geometry.make_valid(); mrc = mrc.to_crs(4326)101 ncol = next((c for c in mrc.columns if "NM_MRC" in c.upper() or c.upper().endswith("NOM") or "MRC_NM" in c.upper()), None)102 if ncol:103 pts = live[live.lat.notna()]; gp = gpd.GeoDataFrame(pts[["id"]], geometry=gpd.points_from_xy(pts.lng, pts.lat), crs=4326)104 jm = gpd.sjoin(gp, mrc[[ncol, "geometry"]], how="left", predicate="within").drop_duplicates("id")105 live = live.merge(jm[["id", ncol]].rename(columns={ncol: "mrc_pip"}), on="id", how="left", suffixes=("_x", ""))106 live = live.drop(columns=[c for c in live.columns if c.endswith("_x")])107 except Exception as ex: # noqa: BLE001108 log.warning("MRC layer not usable: %s", ex)109 transfer = live.mrc_pip.isin(["Brome-Missisquoi", "La Haute-Yamaska"]) & (live.region_db == "Montérégie") & (live.region_pip == "Estrie")110 live["region_final"] = live.region_city.fillna(live.region_db).fillna(live.region_pip).fillna("Inconnue")111 live.loc[transfer, "region_final"] = "Estrie"112 live["region_source"] = np.where(transfer, "pip_transfer2021", np.where(live.region_city.notna(), "city", np.where(live.region_db.notna(), "db", np.where(live.region_pip.notna(), "pip", "unknown"))))113 live["region_pipfirst"] = live.region_pip.fillna(live.region_city).fillna(live.region_db).fillna("Inconnue") # point-in-polygon-first sensitivity variant114 live["region_conflict"] = live.region_db.notna() & live.region_pip.notna() & (live.region_db != live.region_pip) # label vs polygon (any source)115 live["region_final_vs_polygon"] = live.region_pip.notna() & (live.region_final != "Inconnue") & (live.region_final != live.region_pip)116 live["has_geo"] = live.lat.notna()117 ph = live.phone.fillna("").str.replace(r"\D", "", regex=True)118 cnt = ph[ph != ""].value_counts()119 live["phone_shared"] = ph.map(cnt).fillna(0).ge(3)120 live["origin_class"] = live.origin_class.fillna("E")121 live["certified"] = live.n_certifications.fillna(0) > 0122 live["cert_adq"] = live.cert_labels.fillna("").str.contains("Aliments du Québec", regex=False)123 live["cert_lpq"] = live.cert_labels.fillna("").str.contains("Produits du Québec", regex=False)124 live["cert_bio"] = live.cert_labels.fillna("").str.contains("Biologique|CARTV", regex=True)125 live["store_kind"] = live.store_kind.fillna("").replace({"": "standard"})126 live["language"] = live.language.fillna("").replace({"": "unknown", "None": "unknown"})127 live["n_sources"] = live.discovery_sources_list.map(lambda s: len(json.loads(s)))128 live["first_source"] = live.discovery_sources_list.map(lambda s: (json.loads(s) or ["none"])[0])129130 a = li[li.f_analysis_assort]131 agg = a.groupby("store_id").agg(n_listings=("uid", "size"), n_categories=("category", "nunique"),132 share_autre=("category", lambda s: (s == "autre").mean()),133 median_price=("price", "median"), share_available=("available", lambda s: (s == 1).mean()),134 share_on_sale=("on_sale", "mean"), first_listing_utc=("first_seen_utc", "min"),135 last_listing_utc=("first_seen_utc", "max"))136 live = live.merge(agg, left_on="id", right_index=True, how="left")137 live["n_listings"] = live.n_listings.fillna(0).astype(int)138 top = a.groupby(["store_id", "category"]).size().reset_index(name="n").sort_values(["store_id", "n"], ascending=[True, False]).drop_duplicates("store_id")139 live = live.merge(top.rename(columns={"category": "main_category"})[["store_id", "main_category"]], left_on="id", right_on="store_id", how="left").drop(columns=["store_id"])140141 keep = ["id", "name", "platform", "city", "region_db", "region_pip", "region_city", "region_pipfirst", "mrc_pip", "region_final", "region_source", "region_final_vs_polygon", "region_conflict", "has_geo", "lat", "lng",142 "postal_prefix", "phone_shared", "origin_class", "origin_confidence", "origin_evidence", "certified", "cert_adq", "cert_lpq", "cert_bio",143 "cert_labels", "store_kind", "language", "n_sources", "first_source", "reg_status", "product_count", "n_listings", "n_categories",144 "share_autre", "median_price", "share_available", "share_on_sale", "main_category", "first_listing_utc", "last_listing_utc", "last_sync_utc", "fail_streak"]145 out = live[keep].copy()146 out.to_parquet(P["interim"] / "stores_master.parquet", index=False)147 register_output(P["interim"] / "stores_master.parquet", "clean.stores")148 summ = {"live_stores": len(out), "regions_layer": layer_path, "region_source": out.region_source.value_counts().to_dict(),149 "region_conflicts_label_vs_polygon": int(out.region_conflict.sum()), "final_vs_polygon_disagreements": int(out.region_final_vs_polygon.sum()), "transfer2021_overrides": int(transfer.sum()), "unknown_region_listings": int(out.loc[out.region_final == "Inconnue", "n_listings"].sum()),150 "phone_shared_stores": int(out.phone_shared.sum()), "certified_stores": int(out.certified.sum()),151 "region_final_counts": out.region_final.value_counts().to_dict(),152 "agreement_db_vs_pip": (round(float((out.region_db == out.region_pip)[out.region_db.notna() & out.region_pip.notna()].mean()), 4) if out.region_pip.notna().any() else None),153 "agreement_city_vs_pip": (round(float((out.region_city == out.region_pip)[out.region_city.notna() & out.region_pip.notna()].mean()), 4) if out.region_city.notna().any() else None),154 "agreement_db_vs_final": round(float((out.region_db == out.region_final)[out.region_db.notna()].mean()), 4)}155 write_json(summ, BASE_EXP / "store_master.json", "clean.stores")156 ledger("clean.stores", **summ)157 log.info(json.dumps(summ, ensure_ascii=False))158159160if __name__ == "__main__":161 main()162