markov

LISA Markov and Spatial Markov transition analysis.

markov

LISA Markov and Spatial Markov analysis of congestion dynamics.

Computes LISA categories per segment per time period using PySAL (esda), then fits classic Markov and Spatial Markov transition models via giddy to quantify hotspot persistence and spatial contagion.

compute_lisa

compute_lisa(gdf, column='jam_factor_mean', k=8, permutations=999, significance=0.05)

Compute Local Moran's I and assign LISA categories.

Requires libpysal and esda (install with pip install traffic-congestion-pipeline[pysal]).

Returns gdf with added columns: lisa_i, lisa_p, lisa_q, lisa_cat.

Source code in src/trafficpipeline/markov.py

def compute_lisa(
    gdf: gpd.GeoDataFrame,
    column: str = "jam_factor_mean",
    k: int = 8,
    permutations: int = 999,
    significance: float = 0.05,
) -> gpd.GeoDataFrame:
    """Compute Local Moran's I and assign LISA categories.

    Requires ``libpysal`` and ``esda`` (install with
    ``pip install traffic-congestion-pipeline[pysal]``).

    Returns *gdf* with added columns: ``lisa_i``, ``lisa_p``,
    ``lisa_q``, ``lisa_cat``.
    """
    from esda import Moran_Local
    from libpysal.weights import KNN

    gdf = gdf.copy()
    w = KNN.from_dataframe(gdf, k=k)
    w.transform = "r"

    lisa = Moran_Local(gdf[column].values, w, permutations=permutations)

    gdf["lisa_i"] = lisa.Is
    gdf["lisa_p"] = lisa.p_sim
    gdf["lisa_q"] = lisa.q

    # Quadrant labels: 1=HH, 2=LH, 3=LL, 4=HL
    quad_labels = {1: "HH", 2: "LH", 3: "LL", 4: "HL"}
    gdf["lisa_cat"] = "NS"
    sig = gdf["lisa_p"] < significance
    gdf.loc[sig, "lisa_cat"] = gdf.loc[sig, "lisa_q"].map(quad_labels)

    return gdf

compute_global_moran

compute_global_moran(gdf, column='jam_factor_mean', k=8, permutations=999)

Compute Global Moran's I.

Returns dict with morans_i, z_score, p_value.

Source code in src/trafficpipeline/markov.py

def compute_global_moran(
    gdf: gpd.GeoDataFrame,
    column: str = "jam_factor_mean",
    k: int = 8,
    permutations: int = 999,
) -> dict:
    """Compute Global Moran's I.

    Returns dict with ``morans_i``, ``z_score``, ``p_value``.
    """
    from esda import Moran
    from libpysal.weights import KNN

    w = KNN.from_dataframe(gdf, k=k)
    w.transform = "r"

    mi = Moran(gdf[column].values, w, permutations=permutations)
    return {
        "morans_i": float(mi.I),
        "z_score": float(mi.z_sim),
        "p_value": float(mi.p_sim),
    }

build_lisa_matrix

build_lisa_matrix(city_code, base_dir='.', column='jam_factor_mean', k=8)

Build (n_segments × n_periods) integer matrix of LISA codes.

Aligns segments across periods by stable identifier so that Markov transitions are computed on the same physical segments throughout. Segments missing in any period are dropped.

Returns (y, gdfs) where y[i, t] ∈ {0..4} encodes LISA_STATES and each gdf is the aligned (filtered + sorted) view for that period.

Source code in src/trafficpipeline/markov.py

def build_lisa_matrix(
    city_code: str,
    base_dir: str | Path = ".",
    column: str = "jam_factor_mean",
    k: int = 8,
) -> tuple[np.ndarray, list[gpd.GeoDataFrame]]:
    """Build (n_segments × n_periods) integer matrix of LISA codes.

    Aligns segments across periods by stable identifier so that Markov
    transitions are computed on the same physical segments throughout.
    Segments missing in any period are dropped.

    Returns ``(y, gdfs)`` where y[i, t] ∈ {0..4} encodes LISA_STATES and
    each gdf is the aligned (filtered + sorted) view for that period.
    """
    base = Path(base_dir)
    folder = base / CITIES[city_code]["traffic_output_dir"]
    raw_gdfs: list[gpd.GeoDataFrame] = []

    for period in TIME_PERIODS:
        fp = folder / f"{period}_{city_code}.gpkg"
        if not fp.exists():
            continue
        gdf = gpd.read_file(str(fp))
        gdf = compute_lisa(gdf, column=column, k=k)
        raw_gdfs.append(gdf)

    if not raw_gdfs:
        raise FileNotFoundError(
            f"No period GeoPackages found for city '{city_code}' in {folder}"
        )

    # Pick a stable identifier column that's consistent across periods.
    id_col = next(
        (c for c in ("osm_composite_id", "segment_id", "fid")
         if all(c in g.columns for g in raw_gdfs)),
        None,
    )
    if id_col is None:
        # Last-resort fallback: align by row index, but warn loudly.
        import warnings as _w
        _w.warn(
            "No stable segment ID column found across periods "
            "(checked: osm_composite_id, segment_id, fid). "
            "Falling back to row-index alignment; this may be incorrect "
            "if period GeoPackages contain different segments."
        )
        min_n = min(len(g) for g in raw_gdfs)
        gdfs = [g.iloc[:min_n].reset_index(drop=True) for g in raw_gdfs]
    else:
        # Intersect IDs across all periods, then filter + sort consistently.
        common: set = set(raw_gdfs[0][id_col])
        for g in raw_gdfs[1:]:
            common &= set(g[id_col])
        if not common:
            raise ValueError(
                f"No segments common to all periods for city '{city_code}'. "
                f"Checked {len(raw_gdfs)} periods on column '{id_col}'."
            )
        gdfs = [
            g[g[id_col].isin(common)]
              .sort_values(id_col, kind="mergesort")
              .reset_index(drop=True)
            for g in raw_gdfs
        ]

    n_seg = len(gdfs[0])
    n_per = len(gdfs)
    y = np.zeros((n_seg, n_per), dtype=int)
    for t, gdf in enumerate(gdfs):
        y[:, t] = gdf["lisa_cat"].map(LISA_CODE).values

    return y, gdfs

classic_markov

classic_markov(y)

Fit a classic (non-spatial) Markov chain.

Returns dict with transition_matrix, steady_state, and persistence (diagonal probabilities).

Source code in src/trafficpipeline/markov.py

def classic_markov(y: np.ndarray) -> dict:
    """Fit a classic (non-spatial) Markov chain.

    Returns dict with ``transition_matrix``, ``steady_state``,
    and ``persistence`` (diagonal probabilities).
    """
    from giddy.markov import Markov

    m = Markov(y)
    p = m.p
    ss = m.steady_state

    return {
        "transition_matrix": p,
        "steady_state": ss,
        "persistence": {LISA_STATES[i]: float(p[i, i]) for i in range(len(LISA_STATES))},
    }

spatial_markov

spatial_markov(y, gdf, k=8, permutations=999)

Fit a Spatial Markov model conditioned on neighbor states.

Reports the maximally-significant homogeneity test across lag classes: we return the maximum chi-squared statistic and its associated p-value, dof, and significance flag at the 0.05 level.

Returns dict with chi2, p_value, dof, significant, and per_lag (list of (chi2, p, dof) tuples per lag class).

Source code in src/trafficpipeline/markov.py

def spatial_markov(
    y: np.ndarray,
    gdf: gpd.GeoDataFrame,
    k: int = 8,
    permutations: int = 999,
) -> dict:
    """Fit a Spatial Markov model conditioned on neighbor states.

    Reports the maximally-significant homogeneity test across lag classes:
    we return the *maximum* chi-squared statistic and its associated
    p-value, dof, and significance flag at the 0.05 level.

    Returns dict with ``chi2``, ``p_value``, ``dof``, ``significant``,
    and ``per_lag`` (list of (chi2, p, dof) tuples per lag class).
    """
    from giddy.markov import Spatial_Markov
    from libpysal.weights import KNN

    w = KNN.from_dataframe(gdf, k=k)
    w.transform = "r"

    sm = Spatial_Markov(y, w, permutations=permutations)

    # giddy >=2.3 exposes per-lag homogeneity tests via .shtest as a
    # list of (chi2, p_value, dof) tuples (one per spatial-lag class).
    # Older versions sometimes returned .chi2 as a scalar or array.
    per_lag: list[tuple[float, float, int]] = []
    if hasattr(sm, "shtest") and sm.shtest is not None:
        for entry in sm.shtest:
            chi2_i = float(np.asarray(entry[0]).item())
            p_i    = float(np.asarray(entry[1]).item())
            dof_i  = int(np.asarray(entry[2]).item())
            per_lag.append((chi2_i, p_i, dof_i))

    if per_lag:
        # Pick the lag class with the largest chi² (most evidence against H0).
        idx = max(range(len(per_lag)), key=lambda i: per_lag[i][0])
        chi2_val, p_val, dof = per_lag[idx]
    else:
        # Fallback for older giddy: try .chi2 directly.
        raw = getattr(sm, "chi2", None)
        if isinstance(raw, (list, tuple)):
            chi2_val = float(max(np.asarray(v).item() for v in raw))
        elif raw is not None:
            chi2_val = float(np.asarray(raw).item())
        else:
            chi2_val = float("nan")
        p_val = float("nan")
        dof = 0

    return {
        "chi2": chi2_val,
        "p_value": p_val,
        "dof": dof,
        "significant": (p_val < 0.05),
        "per_lag": per_lag,
    }

persistence_statistics

persistence_statistics(y)

Compute segment-level persistence metrics.

Returns dict mapping LISA states to {ever_pct, always_pct, avg_periods}.

Source code in src/trafficpipeline/markov.py

def persistence_statistics(y: np.ndarray) -> dict:
    """Compute segment-level persistence metrics.

    Returns dict mapping LISA states to {ever_pct, always_pct, avg_periods}.
    """
    n_seg, n_per = y.shape
    stats: dict[str, dict] = {}

    for code_idx, state in enumerate(LISA_STATES):
        mask = y == code_idx
        ever = mask.any(axis=1).sum()
        always = mask.all(axis=1).sum()
        avg = mask.sum(axis=1).mean()
        stats[state] = {
            "ever_pct": float(ever / n_seg * 100),
            "always_pct": float(always / n_seg * 100),
            "avg_periods": float(avg),
        }

    return stats

plot_transition_matrix

plot_transition_matrix(p, city_name, figures_dir='figures')

Heatmap of transition probability matrix.

Source code in src/trafficpipeline/markov.py

def plot_transition_matrix(
    p: np.ndarray,
    city_name: str,
    figures_dir: str | Path = "figures",
) -> Path:
    """Heatmap of transition probability matrix."""
    fig_dir = Path(figures_dir) / "markov"
    fig_dir.mkdir(parents=True, exist_ok=True)

    fig, ax = plt.subplots(figsize=(7, 6))
    im = ax.imshow(p, cmap="YlOrRd", vmin=0, vmax=1)
    ax.set_xticks(range(len(LISA_STATES)))
    ax.set_yticks(range(len(LISA_STATES)))
    ax.set_xticklabels(LISA_STATES)
    ax.set_yticklabels(LISA_STATES)
    ax.set_xlabel("To")
    ax.set_ylabel("From")
    ax.set_title(f"{city_name} — Transition Matrix", fontweight="bold")

    for i in range(len(LISA_STATES)):
        for j in range(len(LISA_STATES)):
            val = p[i, j]
            color = "white" if val > 0.5 else "black"
            ax.text(j, i, f"{val:.1%}", ha="center", va="center",
                    color=color, fontsize=10)

    fig.colorbar(im, ax=ax, shrink=0.8, label="Probability")
    plt.tight_layout()

    safe = city_name.lower().replace(" ", "_")
    fp = fig_dir / f"{safe}_transition_matrix.png"
    plt.savefig(fp, dpi=150, bbox_inches="tight")
    plt.close()
    return fp

plot_persistence_comparison

plot_persistence_comparison(all_persistence, figures_dir='figures')

Bar chart comparing persistence (diagonal) across cities.

Source code in src/trafficpipeline/markov.py

def plot_persistence_comparison(
    all_persistence: dict[str, dict],
    figures_dir: str | Path = "figures",
) -> Path:
    """Bar chart comparing persistence (diagonal) across cities."""
    fig_dir = Path(figures_dir) / "markov"
    fig_dir.mkdir(parents=True, exist_ok=True)

    states_to_plot = ["HH", "LL", "NS"]
    cities = list(all_persistence.keys())
    names = [CITIES[c]["name"] for c in cities]
    x = np.arange(len(cities))
    w = 0.25

    fig, ax = plt.subplots(figsize=(10, 6))
    colors = {"HH": "#e74c3c", "LL": "#3498db", "NS": "#95a5a6"}

    for i, state in enumerate(states_to_plot):
        vals = [all_persistence[c][state] * 100 for c in cities]
        bars = ax.bar(x + i * w, vals, w, label=f"P({state}→{state})",
                      color=colors[state], alpha=0.85)
        for bar in bars:
            h = bar.get_height()
            ax.text(bar.get_x() + bar.get_width() / 2, h + 0.3,
                    f"{h:.1f}%", ha="center", va="bottom", fontsize=9)

    ax.set_xticks(x + w)
    ax.set_xticklabels(names)
    ax.set_ylabel("Persistence Probability (%)")
    ax.set_title("Hotspot Persistence Comparison", fontweight="bold")
    ax.legend()
    ax.grid(axis="y", alpha=0.3)
    plt.tight_layout()

    fp = fig_dir / "persistence_comparison.png"
    plt.savefig(fp, dpi=150, bbox_inches="tight")
    plt.close()
    return fp

plot_contagion_results

plot_contagion_results(contagion, figures_dir='figures')

Bar chart of spatial contagion chi² with significance markers.

Source code in src/trafficpipeline/markov.py

def plot_contagion_results(
    contagion: dict[str, dict],
    figures_dir: str | Path = "figures",
) -> Path:
    """Bar chart of spatial contagion chi² with significance markers."""
    fig_dir = Path(figures_dir) / "markov"
    fig_dir.mkdir(parents=True, exist_ok=True)

    cities = list(contagion.keys())
    names = [CITIES[c]["name"] for c in cities]
    chi2_vals = [contagion[c]["chi2"] for c in cities]
    p_vals = [contagion[c]["p_value"] for c in cities]

    fig, ax = plt.subplots(figsize=(8, 5))
    bars = ax.bar(names, chi2_vals, color=[CITIES[c]["color"] for c in cities],
                  alpha=0.85)

    for bar, pv in zip(bars, p_vals):
        h = bar.get_height()
        sig = "***" if pv < 0.001 else "**" if pv < 0.01 else "*" if pv < 0.05 else "ns"
        ax.text(bar.get_x() + bar.get_width() / 2, h + 0.1,
                f"p={pv:.3f} {sig}", ha="center", va="bottom", fontsize=10)

    ax.set_ylabel("Chi-squared statistic")
    ax.set_title("Spatial Contagion Test", fontweight="bold")
    ax.axhline(y=3.84, ls="--", color="gray", alpha=0.5, label="χ² critical (p=0.05)")
    ax.legend()
    ax.grid(axis="y", alpha=0.3)
    plt.tight_layout()

    fp = fig_dir / "spatial_contagion_test.png"
    plt.savefig(fp, dpi=150, bbox_inches="tight")
    plt.close()
    return fp

run_analysis

run_analysis(base_dir='.', figures_dir='figures', output_dir='analysis_results')

Run LISA Markov and Spatial Markov analysis for all cities.

Source code in src/trafficpipeline/markov.py

def run_analysis(
    base_dir: str | Path = ".",
    figures_dir: str | Path = "figures",
    output_dir: str | Path = "analysis_results",
) -> None:
    """Run LISA Markov and Spatial Markov analysis for all cities."""
    out_dir = Path(output_dir)
    out_dir.mkdir(exist_ok=True)

    all_classic: dict[str, dict] = {}
    all_contagion: dict[str, dict] = {}
    all_persistence_diag: dict[str, dict] = {}
    summary_rows: list[dict] = []

    for code in CITIES:
        name = CITIES[code]["name"]
        print(f"\n{'=' * 50}")
        print(f"LISA Markov Analysis: {name}")
        print(f"{'=' * 50}")

        # Build LISA matrix
        print("  Computing LISA across time periods …")
        y, gdfs = build_lisa_matrix(code, base_dir)

        # Global Moran's I (evening peak)
        ep_idx = TIME_PERIODS.index("evening_peak")
        if ep_idx < len(gdfs):
            mi = compute_global_moran(gdfs[ep_idx])
            print(f"  Global Moran's I = {mi['morans_i']:.4f} "
                  f"(z={mi['z_score']:.2f}, p={mi['p_value']:.3f})")

        # Classic Markov
        print("  Fitting classic Markov chain …")
        cm = classic_markov(y)
        all_classic[code] = cm
        all_persistence_diag[code] = cm["persistence"]
        plot_transition_matrix(cm["transition_matrix"], name, figures_dir)

        diag = cm["persistence"]
        print(f"    P(HH→HH)={diag['HH']:.1%}  "
              f"P(LL→LL)={diag['LL']:.1%}  "
              f"P(NS→NS)={diag['NS']:.1%}")

        # Persistence statistics
        ps = persistence_statistics(y)
        print(f"    HH ever: {ps['HH']['ever_pct']:.1f}%  "
              f"always: {ps['HH']['always_pct']:.1f}%  "
              f"avg: {ps['HH']['avg_periods']:.1f} periods")

        # Spatial Markov
        print("  Fitting Spatial Markov model …")
        sm = spatial_markov(y, gdfs[0])
        all_contagion[code] = sm
        sig_str = "SIGNIFICANT" if sm["significant"] else "not significant"
        print(f"    χ²={sm['chi2']:.2f}  p={sm['p_value']:.3f}  ({sig_str})")

        summary_rows.append({
            "city": name,
            "n_segments": y.shape[0],
            "P_HH_HH": diag["HH"],
            "P_LL_LL": diag["LL"],
            "P_NS_NS": diag["NS"],
            "contagion_chi2": sm["chi2"],
            "contagion_p": sm["p_value"],
            "contagion_significant": sm["significant"],
            "morans_i": mi["morans_i"] if ep_idx < len(gdfs) else np.nan,
            "morans_p": mi["p_value"] if ep_idx < len(gdfs) else np.nan,
        })

    # Comparative figures
    plot_persistence_comparison(all_persistence_diag, figures_dir)
    plot_contagion_results(all_contagion, figures_dir)

    # Save summary CSV
    pd.DataFrame(summary_rows).to_csv(
        out_dir / "markov_analysis_results.csv", index=False
    )
    print("\nMarkov analysis complete.")