04 — GIS & Spatial Analytics

Exact-code spatial join, Queen adjacency, dan permutation test

Notebook ini adalah POC portfolio yang dapat dijalankan ulang dari artefak lokal Pasuruan Lens. Kode produksi tetap berada di python/scripts/ dan python/src/pasuruan365/; notebook berfungsi sebagai narasi analitik yang ringkas, transparan, dan mudah dipresentasikan.

from pathlib import Path
import json
import sys

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns


def find_project_root(start: Path | None = None) -> Path:
    start = (start or Path.cwd()).resolve()
    for candidate in (start, *start.parents):
        if (candidate / "package.json").exists() and (candidate / "data").exists():
            return candidate
    raise RuntimeError("Root Pasuruan Lens tidak ditemukan. Jalankan notebook dari repository ini.")


ROOT = find_project_root()
sys.path.insert(0, str(ROOT / "python" / "src"))

pd.set_option("display.max_columns", 100)
pd.set_option("display.max_colwidth", 100)
sns.set_theme(style="whitegrid", context="notebook")
COLORS = {"A": "#2563eb", "B": "#f59e0b", "C": "#10b981"}

print(f"Project root: {ROOT}")
print(f"Python: {sys.version.split()[0]}")

Project root: PASURUAN365
Python: 3.12.6

Tujuan

Geometri dipakai sebagai konteks analitik setelah model selesai. Lokasi, kecamatan, dan adjacency tidak dimasukkan ke clustering. Dengan begitu, pola spasial berfungsi sebagai bukti deskriptif eksternal—bukan target yang secara otomatis dipaksakan oleh model.

import geopandas as gpd

geo = gpd.read_file(ROOT / "public/data/gis/village_boundaries.geojson")
geo["village_id"] = geo["villageId"].astype(str)
assignments = pd.read_csv(ROOT / "data/modeling/cluster_assignments_v2.csv", dtype={"village_id": str})
edges = pd.read_csv(ROOT / "data/spatial/weights/village_adjacency.csv", dtype=str)
summary = json.loads((ROOT / "data/spatial/statistics/spatial_analysis_summary.json").read_text())

pd.Series({
    "polygon": len(geo),
    "unique GIS ID": geo.village_id.nunique(),
    "matched assignment": geo.village_id.isin(assignments.village_id).sum(),
    "CRS browser": str(geo.crs),
    "Queen directed edges": edges.weight_type.eq("QUEEN").sum(),
}).to_frame("nilai")

	nilai
polygon	365
unique GIS ID	365
matched assignment	365
CRS browser	EPSG:4326
Queen directed edges	2084

POC 1 — Visual join memakai village_id

mapped = geo.merge(assignments[["village_id", "cluster_label"]], on="village_id", validate="one_to_one")
mapped["typology"] = mapped.cluster_label.str.replace("Typology ", "", regex=False)

fig, ax = plt.subplots(figsize=(9, 9))
mapped.plot(color=mapped["typology"].map(COLORS), edgecolor="white", linewidth=.18, ax=ax)
ax.set_title("Pasuruan Lens — tipologi pada 365 poligon BIG")
ax.set_axis_off()
plt.tight_layout()

POC 2 — Hitung ulang same-type adjacency

queen = edges.loc[edges.weight_type.eq("QUEEN"), ["village_id", "neighbor_id"]].copy()
labels = assignments.set_index("village_id")["cluster_id"]
queen["same_type"] = queen.village_id.map(labels).eq(queen.neighbor_id.map(labels))
observed_directed = queen["same_type"].mean()

# Karena file adjacency menyimpan dua arah, pasangan unik dipakai untuk permutation.
queen["pair"] = queen.apply(lambda row: tuple(sorted((row.village_id, row.neighbor_id))), axis=1)
unique_pairs = queen.drop_duplicates("pair")[["village_id", "neighbor_id"]]
observed = np.mean([
    labels.at[left] == labels.at[right]
    for left, right in unique_pairs.itertuples(index=False, name=None)
])

pd.Series({
    "unique Queen edges": len(unique_pairs),
    "observed from data": observed,
    "published observed": summary["queen"]["observed"],
    "published expected": summary["queen"]["expected"],
    "published p-value": summary["queen"]["p_value"],
}).round(6).to_frame("nilai")

	nilai
unique Queen edges	1042.000000
observed from data	0.566219
published observed	0.566219
published expected	0.376074
published p-value	0.001000

POC 3 — Fixed-count permutation

rng = np.random.default_rng(365)
ids = labels.index.to_numpy()
label_values = labels.to_numpy()
left_idx = pd.Series(np.arange(len(ids)), index=ids).loc[unique_pairs.village_id].to_numpy()
right_idx = pd.Series(np.arange(len(ids)), index=ids).loc[unique_pairs.neighbor_id].to_numpy()

simulated = np.empty(999)
for i in range(999):
    permuted = rng.permutation(label_values)  # jumlah A/B/C tetap
    simulated[i] = np.mean(permuted[left_idx] == permuted[right_idx])

p_value = (1 + np.sum(simulated >= observed)) / (len(simulated) + 1)

fig, ax = plt.subplots(figsize=(8, 4))
sns.histplot(simulated, bins=28, ax=ax, color="#94a3b8")
ax.axvline(observed, color="#dc2626", linewidth=2, label=f"observed = {observed:.3f}")
ax.set(title="Null distribution: fixed-count label permutation", xlabel="Proporsi pasangan bertipologi sama")
ax.legend()
plt.tight_layout()

print(f"Expected (POC): {simulated.mean():.3f}")
print(f"Observed:       {observed:.3f}")
print(f"p-value:        {p_value:.3f}")

Expected (POC): 0.376
Observed:       0.566
p-value:        0.001

Batas interpretasi

Tipologi yang sama lebih sering bertetangga daripada ekspektasi acak dengan jumlah label tetap. Ini tidak membuktikan sebab-akibat, kualitas model, akses jalan, interaksi ekonomi, atau kedekatan waktu tempuh. Hasil juga bergantung pada unit wilayah dan definisi bobot spasial (MAUP dan specification dependence).