Merge pull request #38 from henryspatialanalysis/feature/static_data_viz

Add static viz and table creation scripts.

Author: njhenry

scripts/conflation/plot_source_contributions.py

@@ -0,0 +1,231 @@
+#!/usr/bin/env python
+"""
+Plot confidence-weighted source contributions to the conflated dataset.
+
+Reads conflated.parquet and draws a horizontal stacked bar chart. Each row
+is a shared_label; the 15 most common shared labels (by POI count) in the
+most recent database are shown, ordered descending so the largest sits at
+the top. Each bar is split left-to-right into the three provenance classes
+of the conflated dataset:
+
+    Overture only (left)  — source == "overture"
+    Both          (mid)   — source == "matched"
+    OSM only      (right) — source == "osm"
+
+Bar lengths are **confidence-weighted** observation counts: each POI is
+weighted by its combined confidence score (``conf_mean``, the final
+post-change-detection blended confidence carried in the published dataset)
+rather than counted as 1. Each colored sub-bar is annotated, in white, with
+its confidence-weighted contribution rounded roughly to the nearest
+thousand.
+
+Blank/null shared labels are excluded — they are not a meaningful label
+type — so the 15 rows are the 15 most common *named* labels.
+
+Config keys used (config.yaml):
+    conflation.conflated — input GeoParquet path (conflated.parquet)
+    conflation            — directory; output PNG lands in its viz/ subdir
+
+Prerequisites:
+    Run scripts/conflation/conflate.py (and apply_change_detection.py) first.
+
+Output file (in conflation/<version>/viz/):
+    source_contributions.png
+"""
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+from config_versioned import Config
+
+from pathlib import Path  # noqa: E402
+
+import matplotlib
+matplotlib.use("Agg")  # noqa: E402
+import matplotlib.font_manager as fm  # noqa: E402
+import matplotlib.pyplot as plt  # noqa: E402
+from matplotlib.ticker import FuncFormatter, MultipleLocator  # noqa: E402
+
+# ----------------------------------------------------------------------------------------
+# Configuration constants
+# ----------------------------------------------------------------------------------------
+
+config = Config("~/repos/openpois/config.yaml")
+INPUT_PATH = config.get_file_path("conflation", "conflated")
+VIZ_DIR = config.get_dir_path("conflation") / "viz"
+OUTPUT_PATH = VIZ_DIR / "source_contributions.png"
+
+TOP_N = 14
+
+# Shared labels to drop before ranking: every "Other ..." catch-all plus a
+# few named labels we don't want surfaced here.
+EXCLUDE_EXACT = {
+    "Home Service", "Swimming Pool", "Real Estate", "Specialty Store",
+    "Hotel", "Recreation",
+}
+
+# Register the Figtree variable font from the first candidate path that
+# exists, then make it the default family. Falls back to matplotlib's
+# default if none are found.
+FIGTREE_CANDIDATES = [
+    Path("/mnt/c/Users/nathe/AppData/Local/Microsoft/Windows/Fonts/"
+         "Figtree-VariableFont_wght.ttf"),
+    Path("/mnt/d/Users/Lenovo/AppData/Local/Microsoft/Windows/Fonts/"
+         "Figtree-VariableFont_wght.ttf"),
+]
+for _font in FIGTREE_CANDIDATES:
+    if _font.exists():
+        fm.fontManager.addfont(str(_font))
+        plt.rcParams["font.family"] = "Figtree"
+        break
+
+# Base font size, bumped 20% over the matplotlib default of 10.
+plt.rcParams["font.size"] = 12
+
+# source value -> (legend label, color). Listed left-to-right.
+SOURCES = [
+    ("overture", "Overture only", "#2e86c9"),
+    ("matched", "Both", "#3d00a5"),
+    ("osm", "OSM only", "#a0d787"),
+]
+
+# Skip the in-bar number when a sub-bar is narrower than this fraction of the
+# widest total bar — the text would overflow its segment and collide.
+LABEL_MIN_FRAC = 0.018
+
+# ----------------------------------------------------------------------------------------
+# Helpers
+# ----------------------------------------------------------------------------------------
+
+
+def fmt_weighted(value: float) -> str:
+    """Format a confidence-weighted count roughly to the nearest thousand.
+
+    >= 1k  -> integer thousands with a "k" suffix ("837k", "17k", "9k")
+    < 1k   -> raw integer ("751")
+    """
+    if value >= 1_000:
+        return f"{round(value / 1000):d}k"
+    return f"{round(value):d}"
+
+
+# ----------------------------------------------------------------------------------------
+# Main workflow
+# ----------------------------------------------------------------------------------------
+
+if __name__ == "__main__":
+    print(f"Reading {INPUT_PATH} ...")
+    df = pd.read_parquet(
+        INPUT_PATH,
+        columns = ["shared_label", "source", "conf_mean"],
+    )
+    print(f"  {len(df):,} rows")
+
+    # Drop blank/null shared labels — not a meaningful label type — plus the
+    # excluded "Other ..." catch-alls and named labels.
+    label = df["shared_label"]
+    keep = (
+        label.notna()
+        & (label != "")
+        & ~label.str.startswith("Other ")
+        & ~label.isin(EXCLUDE_EXACT)
+    )
+    df = df[keep]
+
+    # Confidence-weighted sum and raw count per (shared_label, source).
+    agg = df.groupby(["shared_label", "source"], observed = True).agg(
+        wsum = ("conf_mean", "sum"),
+        n = ("conf_mean", "size"),
+    )
+    wsum = agg["wsum"].unstack(fill_value = 0.0)
+    counts = agg["n"].unstack(fill_value = 0)
+    for source, _, _ in SOURCES:
+        if source not in wsum.columns:
+            wsum[source] = 0.0
+            counts[source] = 0
+
+    # Top N shared labels by total confidence-weighted POI count; largest at
+    # the top of the chart.
+    top_labels = (
+        wsum.sum(axis = 1).sort_values(ascending = False).head(TOP_N).index
+    )
+    wsum = wsum.loc[top_labels]
+
+    # ----------------------------------------------------------------------------
+    # Draw the horizontal stacked bar chart
+    # ----------------------------------------------------------------------------
+    y = np.arange(len(top_labels))
+    max_total = wsum.sum(axis = 1).max()
+    label_threshold = max_total * LABEL_MIN_FRAC
+
+    fig, ax = plt.subplots(figsize = (13.33, 7.5))
+
+    left = np.zeros(len(top_labels))
+    for source, _, color in SOURCES:
+        widths = wsum[source].to_numpy()
+        ax.barh(y, widths, left = left, color = color, height = 0.90)
+        # White in-bar annotation at each segment's center.
+        for yi, w, x0 in zip(y, widths, left):
+            if w >= label_threshold:
+                ax.text(
+                    x0 + w / 2,
+                    yi,
+                    fmt_weighted(w),
+                    ha = "center",
+                    va = "center",
+                    color = "white",
+                    fontsize = 10,
+                    fontweight = "bold",
+                )
+        left += widths
+
+    ax.set_yticks(y)
+    ax.set_yticklabels(top_labels)
+    ax.invert_yaxis()  # largest label at the top
+    ax.set_xlabel("Confidence weighted POI count")
+    # Title 20% larger than matplotlib's default (1.2x base -> 1.44x base).
+    ax.set_title(
+        "Confidence-weighted records in the conflated dataset, by source",
+        fontsize = plt.rcParams["font.size"] * 1.44,
+    )
+
+    # Light-grey panel with white gridlines at each 100k tick, drawn behind
+    # the bars. Ticks labelled "100k", "200k", ... rather than in millions.
+    ax.set_xlim(0, left.max() * 1.02)
+    ax.xaxis.set_major_locator(MultipleLocator(100_000))
+    ax.xaxis.set_major_formatter(
+        FuncFormatter(lambda v, _: "0" if v == 0 else f"{v / 1000:g}k")
+    )
+    ax.set_facecolor("white")
+    ax.set_axisbelow(True)
+    ax.grid(axis = "x", color = "#DDDDDD", linewidth = 1.0)
+    ax.grid(axis = "y", visible = False)
+
+    # Drop the border (all four spines) and tick marks.
+    for spine in ax.spines.values():
+        spine.set_visible(False)
+    ax.tick_params(length = 0)
+
+    # Single-row legend at the bottom.
+    handles = [
+        plt.Rectangle((0, 0), 1, 1, color = color)
+        for _, _, color in SOURCES
+    ]
+    labels = [legend for _, legend, _ in SOURCES]
+    ax.legend(
+        handles,
+        labels,
+        loc = "upper center",
+        bbox_to_anchor = (0.5, -0.065),
+        ncol = len(SOURCES),
+        frameon = False,
+    )
+
+    # Low padding so the plot fills the fixed 13.33 x 7.5" canvas.
+    fig.subplots_adjust(
+        left = 0.14, right = 0.99, top = 0.94, bottom = 0.11,
+    )
+    VIZ_DIR.mkdir(parents = True, exist_ok = True)
+    # No bbox_inches="tight" — keep the canvas at the exact 13.33 x 7.5".
+    fig.savefig(OUTPUT_PATH, dpi = 300)
+    print(f"\nSaved to {OUTPUT_PATH}")

scripts/models/city_confidence_quartiles.py

@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+"""
+Confidence quartiles by city x shared_label for cities over 100k population.
+
+Uses the per-POI ``conf_mean`` written onto the rated OSM snapshot by the best
+random-effects turnover model (``2026-06-05-nationwide-full``; ``conf_mean`` =
+1 - p_change at the POI's age, for that POI's shared_label x MSA x urbanicity
+cell). Each POI's ``model_group`` encodes ``"shared_label | msa_code |
+urban_rural"``; the shared_label is the first field.
+
+Steps
+-----
+1. Census places (cb_2023_us_place_500k) joined to 2020 decennial population,
+   filtered to places with population > 100,000.
+2. Each rated POI is assigned to a place by its centroid (point-in-polygon).
+3. For each city, compute confidence quartiles (min, Q1, median, Q3, max) plus
+   count and mean, both overall (all POIs in the city, ``shared_label =
+   __ALL__``) and per shared_label.
+
+Output: one row per city x shared_label (plus one __ALL__ row per city).
+"""
+import os
+import sys
+from pathlib import Path
+
+os.environ.setdefault("PROJ_DATA", str(Path(sys.prefix) / "share" / "proj"))
+os.environ.setdefault("PROJ_LIB", str(Path(sys.prefix) / "share" / "proj"))
+
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+
+CENSUS_DIR = Path("~/data/openpois/census_areas").expanduser()
+PLACE_SHP = CENSUS_DIR / "cb_2023_us_place_500k.shp"
+POP_CSV = CENSUS_DIR / "place_population_2020.csv"
+RATED_SNAPSHOT = Path(
+    "~/data/openpois/snapshots/osm/20260521/osm_snapshot_rated.parquet"
+).expanduser()
+POP_THRESHOLD = 100_000
+OUT_CSV = Path(
+    "~/data/openpois/osm_turnover_model/2026-06-05-nationwide-full/"
+    "city_shared_label_confidence_quartiles.csv"
+).expanduser()
+
+ALL_LABEL = "__ALL__"
+
+
+def load_big_cities() -> gpd.GeoDataFrame:
+    """Census places with 2020 population > threshold, in EPSG:4326."""
+    places = gpd.read_file(PLACE_SHP)[
+        ["GEOID", "NAME", "NAMELSAD", "STUSPS", "STATE_NAME", "geometry"]
+    ]
+    pop = pd.read_csv(POP_CSV, dtype = {"place_geoid": str})
+    pop["population"] = pd.to_numeric(pop["population"], errors = "coerce")
+    places = places.merge(
+        pop, left_on = "GEOID", right_on = "place_geoid", how = "inner"
+    )
+    big = places[places["population"] > POP_THRESHOLD].copy()
+    big = big.to_crs("EPSG:4326")
+    print(f"{len(big)} places with population > {POP_THRESHOLD:,}")
+    return big
+
+
+def quantile_table(df: pd.DataFrame, group_cols: list[str]) -> pd.DataFrame:
+    """count / min / Q1 / median / Q3 / max / mean of conf_mean per group."""
+    g = df.groupby(group_cols, observed = True)["conf_mean"]
+    out = g.agg(
+        n_pois = "size",
+        conf_min = "min",
+        conf_q1 = lambda s: s.quantile(0.25),
+        conf_median = "median",
+        conf_q3 = lambda s: s.quantile(0.75),
+        conf_max = "max",
+        conf_mean = "mean",
+    ).reset_index()
+    return out
+
+
+def main() -> None:
+    big = load_big_cities()
+
+    print(f"Reading rated snapshot from {RATED_SNAPSHOT} ...")
+    pois = gpd.read_parquet(
+        RATED_SNAPSHOT, columns = ["geometry", "conf_mean", "model_group"]
+    )
+    print(f"  {len(pois):,} POIs")
+
+    # Assign every POI a representative point (building polygons -> centroid).
+    pois["geometry"] = pois.geometry.centroid
+    pois["shared_label"] = (
+        pois["model_group"].str.split(" | ", regex = False).str[0]
+    )
+    pois = pois.drop(columns = "model_group")
+
+    print("Spatial join POIs -> big cities ...")
+    joined = gpd.sjoin(
+        pois, big[["GEOID", "NAME", "STUSPS", "population", "geometry"]],
+        how = "inner", predicate = "within",
+    )
+    print(f"  {len(joined):,} POIs fall within a >100k city")
+    joined = pd.DataFrame(joined.drop(columns = "geometry"))
+
+    city_cols = ["GEOID", "NAME", "STUSPS", "population"]
+
+    # Per city x shared_label, plus an __ALL__ row per city (all POIs in city).
+    by_label = quantile_table(joined, city_cols + ["shared_label"])
+    overall = quantile_table(joined, city_cols)
+    overall["shared_label"] = ALL_LABEL
+
+    result = pd.concat([overall, by_label], ignore_index = True)
+    result = result.rename(
+        columns = {"GEOID": "city_geoid", "NAME": "city_name", "STUSPS": "state"}
+    )
+    # Order: __ALL__ first within each city, then labels by descending median.
+    result["_is_all"] = (result["shared_label"] == ALL_LABEL).astype(int)
+    result = result.sort_values(
+        ["city_name", "state", "_is_all", "conf_median"],
+        ascending = [True, True, False, False],
+    ).drop(columns = "_is_all")
+
+    cols = [
+        "city_geoid", "city_name", "state", "population", "shared_label",
+        "n_pois", "conf_min", "conf_q1", "conf_median", "conf_q3", "conf_max",
+        "conf_mean",
+    ]
+    result = result[cols]
+    OUT_CSV.parent.mkdir(parents = True, exist_ok = True)
+    result.to_csv(OUT_CSV, index = False)
+    print(f"\nWrote {len(result):,} rows to {OUT_CSV}")
+
+    # Brief summary: highest / lowest median-confidence shared_label per city,
+    # aggregated nationally for a quick read.
+    lab = by_label.rename(columns = {"NAME": "city_name"})
+    lab = lab[lab["n_pois"] >= 20]
+    nat = (
+        lab.groupby("shared_label")
+        .apply(
+            lambda d: np.average(d["conf_median"], weights = d["n_pois"]),
+            include_groups = False,
+        )
+        .rename("wtd_median_conf")
+        .reset_index()
+        .sort_values("wtd_median_conf", ascending = False)
+    )
+    print("\nMost stable shared_labels (pop-weighted median conf, n>=20/city):")
+    print(nat.head(8).to_string(index = False))
+    print("\nLeast stable shared_labels:")
+    print(nat.tail(8).to_string(index = False))
+
+
+if __name__ == "__main__":
+    main()