import sys, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from IPython.display import Markdown, display

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)


def find_repo_root(start: Path = Path.cwd()) -> Path:
    for p in [start, *start.parents]:
        if (p / "pyproject.toml").exists() and (p / "data").is_dir():
            return p
    raise FileNotFoundError("Run this notebook from inside the PySpatialStats repository.")


ROOT = find_repo_root()
DATA = ROOT / "data"
try:
    import spatialstats
except ImportError:
    sys.path.insert(0, str(ROOT))
    import spatialstats

from spatialstats.pointpattern import (
    Window, PointPattern, envelope, nearest_neighbor_distances, g_function, ripley_k, ripley_l,
    clark_evans, cross_k,
)
from spatialstats import Queen, Moran

plt.rcParams.update({
    "figure.dpi": 100, "savefig.dpi": 100, "font.size": 10, "axes.titlesize": 11,
    "axes.titleweight": "bold", "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25, "legend.frameon": False,
})
C = {"blue": "#2b6cb0", "red": "#c0392b", "green": "#2f855a", "orange": "#dd6b20",
     "grey": "#718096", "purple": "#6b46c1", "teal": "#2c7a7b"}

print(f"spatialstats : version {spatialstats.__version__}")
print(f"data folder  : {DATA.relative_to(ROOT)}/")

spatialstats : version 0.2.0
data folder  : data/

from pyproj import Transformer

tbl = pd.read_csv(DATA / "crime_data_2018.csv")
city = gpd.read_file(DATA / "buffalo_city.shp")

CRS_CORRECT = 26917
city_fixed = city.to_crs(4326).to_crs(CRS_CORRECT)
transformer = Transformer.from_crs(4326, CRS_CORRECT, always_xy=True)
x_m, y_m = transformer.transform(tbl["long"].to_numpy(), tbl["lat"].to_numpy())
coords_all = np.column_stack([x_m, y_m])

win_city = Window.from_gdf(city_fixed)
pp_all = PointPattern(coords_all, window=win_city, marks=tbl["incident_type"].to_numpy(), clip=True)
pp_dedup = pp_all.deduplicate()

pts_gdf = gpd.GeoDataFrame(
    {"incident_type": tbl["incident_type"]}, geometry=gpd.points_from_xy(x_m, y_m), crs=CRS_CORRECT
)

print(f"crime pattern    : n = {pp_all.n:,}  (deduplicated: {pp_dedup.n:,})")
print(f"neighbourhoods   : {len(city_fixed)}")
print(f"neighbourhood fields: {[c for c in city_fixed.columns if c != 'geometry']}")

crime pattern    : n = 13,948  (deduplicated: 8,903)
neighbourhoods   : 35
neighbourhood fields: ['calcacres', 'nbhdname', 'nbhdnum', 'objectid', 'objectid_1', 'placename', 'sde_vector', 'shape_leng', 'shape_star', 'shape_stle', 'sqmiles']

joined = gpd.sjoin(pts_gdf, city_fixed[["nbhdname", "geometry"]], predicate="within", how="left")
counts = joined.groupby("nbhdname").size().reindex(city_fixed["nbhdname"]).fillna(0).astype(int)

city_fixed = city_fixed.assign(
    crime_count=counts.to_numpy(),
    area_km2=city_fixed.geometry.area / 1e6,
)
city_fixed["rate_km2"] = city_fixed["crime_count"] / city_fixed["area_km2"]

n_matched = joined["nbhdname"].notna().sum()
n_unmatched = len(joined) - n_matched
print(f"events joined to a neighbourhood: {n_matched:,} of {len(joined):,} "
      f"({100 * n_matched / len(joined):.1f}%; the remaining {n_unmatched} fall outside every neighbourhood "
      "polygon, the same city-edge points Chapter 1.1's clip=True dropped)")

city_fixed[["nbhdname", "area_km2", "crime_count", "rate_km2"]].sort_values(
    "crime_count", ascending=False
).head(5).round(2)

events joined to a neighbourhood: 13,948 of 14,001 (99.6%; the remaining 53 fall outside every neighbourhood polygon, the same city-edge points Chapter 1.1's clip=True dropped)

fig, axes = plt.subplots(1, 2, figsize=(14, 6.4))
city_fixed.plot(column="crime_count", cmap="YlOrRd", edgecolor="k", linewidth=0.4, legend=True,
                legend_kwds={"label": "crime count", "shrink": 0.75}, ax=axes[0])
city_fixed.plot(column="rate_km2", cmap="YlOrRd", edgecolor="k", linewidth=0.4, legend=True,
                legend_kwds={"label": "crimes / km²", "shrink": 0.75}, ax=axes[1])
axes[0].set_title("Raw crime count by neighbourhood"); axes[1].set_title("Crime rate by neighbourhood")
for ax in axes:
    ax.set_axis_off()
plt.tight_layout(); plt.show()

print(city_fixed[["area_km2", "crime_count", "rate_km2"]].describe().round(2))

       area_km2  crime_count  rate_km2
count     35.00        35.00     35.00
mean       3.02       398.51    151.96
std        1.85       200.64     79.67
min        0.86        84.00     34.74
25%        1.89       226.50    105.57
50%        2.33       402.00    138.36
75%        3.72       516.00    182.90
max       10.74       942.00    415.41

w = Queen(city_fixed)
print(w.summary())

mi_count = Moran(city_fixed["crime_count"].to_numpy(), w, permutations=999, seed=1)
mi_rate = Moran(city_fixed["rate_km2"].to_numpy(), w, permutations=999, seed=1)

comparison = pd.DataFrame([
    {"variable": "crime_count (raw)", "Moran's I": mi_count["I"], "z": mi_count["z_score"],
     "p (999 perms)": mi_count["p_value"]},
    {"variable": "rate_km2 (area-normalised)", "Moran's I": mi_rate["I"], "z": mi_rate["z_score"],
     "p (999 perms)": mi_rate["p_value"]},
]).set_index("variable")
comparison.round(4)

Spatial Weights (W)
===================
  Observations        35
  Transform           row
  Min neighbors       2
  Mean neighbors      4.800
  Max neighbors       8
  Islands             0
  Components          1
  Nonzero weights     168

corr_area_count = np.corrcoef(city_fixed["area_km2"], city_fixed["crime_count"])[0, 1]
corr_area_rate = np.corrcoef(city_fixed["area_km2"], city_fixed["rate_km2"])[0, 1]
print(f"correlation(area, raw count) : {corr_area_count:+.3f}")
print(f"correlation(area, rate)      : {corr_area_rate:+.3f}")
display(Markdown(
    "Raw count is meaningfully correlated with polygon area; the rate, by construction, is not. **Always "
    "normalise counts before testing for spatial pattern across areas of unequal size** — the single "
    "most important habit to carry from Part 1 into Part 2."
))

correlation(area, raw count) : +0.401
correlation(area, rate)      : -0.426

coords_wrong = tbl[["x_m", "y_m"]].to_numpy()   # native declared-CRS coordinates, used as-is
win_wrong = Window.from_gdf(city)               # the city boundary, also NOT reprojected

rng = np.random.default_rng(5)
idx = rng.choice(len(tbl), 1500, replace=False)
sub_right = PointPattern(coords_all[idx], window=win_city, clip=True)
sub_wrong = PointPattern(coords_wrong[idx], window=win_wrong, clip=True)

nn_right = nearest_neighbor_distances(sub_right)
nn_wrong = nearest_neighbor_distances(sub_wrong)
ce_right = clark_evans(sub_right)
ce_wrong = clark_evans(sub_wrong)

print(f"mean nearest-neighbour distance,  corrected CRS (EPSG:26917) : {nn_right.mean():.1f} m")
print(f"mean nearest-neighbour distance,  declared CRS  (as-is)      : {nn_wrong.mean():.1f} m")
print(f"ratio (declared / corrected)                                  : {nn_wrong.mean() / nn_right.mean():.3f}"
      "  (Part 0 §0.5's independent CRS audit measured a 1.121× distance scale factor for this file)")
print()
print(f"Clark–Evans R,  corrected CRS : {ce_right['R']:.3f}")
print(f"Clark–Evans R,  declared CRS  : {ce_wrong['R']:.3f}   (a ratio of two distances — nearly unaffected)")

mean nearest-neighbour distance,  corrected CRS (EPSG:26917) : 93.3 m
mean nearest-neighbour distance,  declared CRS  (as-is)      : 105.4 m
ratio (declared / corrected)                                  : 1.129  (Part 0 §0.5's independent CRS audit measured a 1.121× distance scale factor for this file)

Clark–Evans R,  corrected CRS : 0.687
Clark–Evans R,  declared CRS  : 0.693   (a ratio of two distances — nearly unaffected)

theft = pp_dedup.by_mark("Theft")
rng_ex1 = np.random.default_rng(9)
theft_sample_coords = theft.subset(rng_ex1.choice(theft.n, 1200, replace=False)).coords

win_hull_ex1 = Window.convex_hull(theft_sample_coords)
pp_city_ex1 = PointPattern(theft_sample_coords, window=win_city, clip=True)
pp_hull_ex1 = PointPattern(theft_sample_coords, window=win_hull_ex1, clip=True)

r_ex1 = np.linspace(0, 1200, 30)
env_city = envelope(pp_city_ex1, "L", r=r_ex1, nsim=99, seed=1)
env_hull = envelope(pp_hull_ex1, "L", r=r_ex1, nsim=99, seed=1)

print(f"city-boundary window : p_mad = {env_city.test('mad'):.3f}, p_dclf = {env_city.test('dclf'):.3f}")
print(f"convex-hull window   : p_mad = {env_hull.test('mad'):.3f}, p_dclf = {env_hull.test('dclf'):.3f}")

fig, ax = plt.subplots(figsize=(8, 5.5))
ax.plot(r_ex1, env_city.observed - r_ex1, color=C["red"], lw=2, label="city boundary window")
ax.plot(r_ex1, env_hull.observed - r_ex1, color=C["blue"], lw=2, label="convex hull window")
ax.axhline(0, color="k", ls="--", lw=1, label="CSR")
ax.set_xlabel("r (m)"); ax.set_ylabel("L(r) − r")
ax.set_title("Same events, two windows: L(r) − r")
ax.legend()
plt.tight_layout(); plt.show()

city-boundary window : p_mad = 0.010, p_dclf = 0.010
convex-hull window   : p_mad = 0.010, p_dclf = 0.010

r_ex2 = np.linspace(0, 1000, 30)
k_right = ripley_k(sub_right, r_ex2)
k_wrong = ripley_k(sub_wrong, r_ex2)

i500 = np.searchsorted(r_ex2, 500)
print(f"K(500 m), corrected CRS : {k_right.estimate[i500]:,.0f}")
print(f"K(500 m), declared CRS  : {k_wrong.estimate[i500]:,.0f}")
print(f"ratio (declared / corrected) : {k_wrong.estimate[i500] / k_right.estimate[i500]:.3f}")

K(500 m), corrected CRS : 1,405,313
K(500 m), declared CRS  : 1,441,986
ratio (declared / corrected) : 1.026

rng_ex3 = np.random.default_rng(2)
sub_ex3 = pp_all.subset(rng_ex3.choice(pp_all.n, 3000, replace=False))
sub_dedup_ex3 = sub_ex3.deduplicate()
sub_jitter_ex3 = sub_ex3.jitter(scale=5.0, seed=1)

print(f"n (duplicates left in) : {sub_ex3.n:,}  ({sub_ex3.n_duplicates} sharing a location, "
      f"{100 * sub_ex3.n_duplicates / sub_ex3.n:.1f}%)")
print(f"n (deduplicated)       : {sub_dedup_ex3.n:,}")
print(f"n (jittered ±5 m)     : {sub_jitter_ex3.n:,}  (every event kept)")

r_short = np.linspace(0, 60, 25)
i10 = np.searchsorted(r_short, 10)
for label, pat in [("duplicates left in", sub_ex3), ("deduplicated", sub_dedup_ex3), ("jittered ±5 m", sub_jitter_ex3)]:
    gg = g_function(pat, r_short)
    nn = nearest_neighbor_distances(pat)
    print(f"{label:20s}: G(10 m) = {gg.estimate[i10]:.3f}   exact-zero NN distances = {100 * np.mean(nn == 0):.1f}%   "
          f"min NN distance = {nn.min():.2f} m")

n (duplicates left in) : 3,000  (534 sharing a location, 17.8%)
n (deduplicated)       : 2,466
n (jittered ±5 m)     : 3,000  (every event kept)

duplicates left in  : G(10 m) = 0.288   exact-zero NN distances = 26.8%   min NN distance = 0.00 m
deduplicated        : G(10 m) = 0.031   exact-zero NN distances = 0.0%   min NN distance = 1.64 m
jittered ±5 m       : G(10 m) = 0.286   exact-zero NN distances = 0.0%   min NN distance = 0.11 m

rng_ex4 = np.random.default_rng(3)
sub_ex4 = pp_dedup.subset(rng_ex4.choice(pp_dedup.n, 2000, replace=False))
mask = np.isin(sub_ex4.marks, ["Assault", "Theft of Vehicle"])
pp_two = sub_ex4.subset(mask)
print(pd.Series(pp_two.marks).value_counts())

r_ex4 = np.linspace(0, 1500, 40)

def cross_k_stat(pattern, r):
    """A plain ndarray, as `envelope` requires from a custom statistic."""
    return cross_k(pattern, "Assault", "Theft of Vehicle", r).estimate

def random_labelling(rng):
    """Same locations, marks reshuffled -- the correct null for genuine type-to-type interaction."""
    shuffled_marks = rng.permutation(pp_two.marks)
    return PointPattern(pp_two.coords, window=pp_two.window, marks=shuffled_marks)

env_cross = envelope(pp_two, statistic=cross_k_stat, r=r_ex4, nsim=99, seed=21, simulator=random_labelling)

print(f"\nrandom-labelling test: p_mad = {env_cross.test('mad'):.3f}, p_dclf = {env_cross.test('dclf'):.3f}")
share_below_mean = np.mean(env_cross.observed < env_cross.mean_sim)
print(f"share of r with observed cross-K below the random-labelling mean: {share_below_mean:.2f}")

fig, ax = plt.subplots(figsize=(8, 5.5))
ax.fill_between(r_ex4, env_cross.lo, env_cross.hi, color="0.8", label="random-labelling envelope")
ax.plot(r_ex4, env_cross.mean_sim, "k--", lw=1, label="random-labelling mean")
ax.plot(r_ex4, env_cross.observed, color=C["purple"], lw=2, label="observed cross-K")
ax.set_xlabel("r (m)"); ax.set_ylabel(r"$K_{ij}(r)$")
ax.set_title("Assault ↔ Theft of Vehicle: observed vs. random relabelling of the same locations")
ax.legend()
plt.tight_layout(); plt.show()

Assault             398
Theft of Vehicle    173
Name: count, dtype: int64

random-labelling test: p_mad = 0.010, p_dclf = 0.010
share of r with observed cross-K below the random-labelling mean: 0.85

Step	What we do
Aggregation	Count crime events within each of Buffalo's 35 neighbourhoods
MAUP	Explain why the aggregated result depends on the areal units chosen, not just the data
Bridge to Part 2	Compute Moran's I on the aggregated counts and on the aggregated rates — and show they disagree
CRS pitfall	Directly measure how much the wrong-CRS shapefile would have distorted a distance-based result
Pitfall review	Consolidate every data-quality lesson from Chapters 1.1–1.4 in one place
Exercises	Four solved problems: window choice, CRS impact, deduplication, and a rigorous re-test of Chapter 1.4's cross-K finding

	nbhdname	area_km2	crime_count	rate_km2
22	Broadway Fillmore	5.83	942	161.55
12	Kensington-Bailey	3.04	778	255.73
3	Central	5.97	749	125.56
24	North Park	4.68	648	138.36
0	Genesee-Moselle	3.83	642	167.46

	Moran's I	z	p (999 perms)
variable
crime_count (raw)	-0.0431	-0.1309	0.718
rate_km2 (area-normalised)	0.4442	4.6754	0.001

Pitfall	Where it showed up	What was measured
Wrong CRS	Chapter 1.1 §5.1, this chapter §8	Declared UTM 11N vs. correct UTM 17N: +25.7% area, ~12–13% distance inflation
Window choice	Chapter 1.1 §10	Bounding box 1.25× the true city area; convex hull leaks 3.9 km² outside city limits despite a similar total area
Duplicate locations	Chapter 1.1 §9, Chapter 1.3 §6, Exercise 3 below	36% of events share a geocoded location; left uncorrected, short-range $G(r)$ is inflated by an order of magnitude
Bandwidth / parameter sensitivity	Chapter 1.2 §8	Scott's rule and cross-validation bandwidths differed by a factor of 4.3×, with visibly different maps
Edge correction	Chapter 1.3 §11	The uncorrected $K(r)$ estimate fell up to 12% below the translation-corrected one by $r=1{,}500$ m; the border correction became unusable at the same distance (only 1.9% of events qualified as anchors)
Cumulative-statistic ceiling effects	Chapter 1.3 §8	$G(r)$ and $F(r)$ both approach 1 at long range, making "above/below CSR" comparisons meaningless there; $J(r)$ became numerically undefined beyond ~1.1 km
Reading a pointwise envelope as a test	Chapter 1.4 §4	Explained why only `EnvelopeResult.test()` — not the shaded band — is a valid significance test
Intensity trend confounding "clustering"	Chapter 1.4 §10	The known citywide intensity surface explained up to 70% of the apparent excess in the homogeneous K function at long range
Shared intensity confounding cross-type interaction	Chapter 1.4 §11, Exercise 4 below	Flagged as an open question in Chapter 1.4; resolved with a proper test below
Raw counts vs rates at the areal level	This chapter §7	Moran's I on raw neighbourhood counts: not significant (p = 0.72); on the area-normalised rate: highly significant (p = 0.001)

Reference	Notes
Openshaw, S. (1984). The Modifiable Areal Unit Problem. CATMOG 38, Geo Books.	The foundational statement of MAUP
Bailey, T. C. & Gatrell, A. C. (1995). Interactive Spatial Data Analysis. Longman.	Bridges point pattern analysis to areal and continuous-field methods
Waller, L. A. & Gotway, C. A. (2004). Applied Spatial Statistics for Public Health Data. Wiley.	Small-area rate instability, aggregation, and disease mapping (previews Part 5)
Cressie, N. (1993). Statistics for Spatial Data. Wiley.	All three spatial data types treated in one unified framework

1.5 Integration with Areal Analysis, Pitfalls and Exercises¶

Table of Contents¶

1. Overview¶

2. From Points to Areas: Why Aggregate?¶

3. Small-area Instability and the Modifiable Areal Unit Problem¶

4. Python Setup¶

5. Data: Rebuild the Point Pattern and the Neighbourhood Polygons¶

6. Aggregating Points to Polygons¶

7. Counts vs Rates: A Bridge to Part 2's Autocorrelation Methods¶

8. Pitfall: Computing Distances in the Wrong CRS¶

9. Pitfall Review: Every Trap Met in Part 1¶

10. Exercises¶

10.1 Window Choice: City Boundary vs Convex Hull¶

10.2 Quantifying the CRS Distortion in K(r)¶

10.3 Deduplicate vs Jitter at Short Range¶

10.4 Do Assault and Vehicle Theft Really Co-locate?¶

11. Summary and Conclusions¶

What we did¶

Practical takeaways¶

Where Part 1 leaves off¶

12. Resources¶

Textbooks¶

Key papers¶

In this library¶

Paper	Contribution
Loosmore, N. B. & Ford, E. D. (2006). Statistical inference using the G or K point pattern spatial statistics. Ecology, 87, 1925–1931.	Random labelling as a null for marked-pattern interaction (Exercise 4)
Diggle, P. J. & Cox, T. F. (1983). Some distance-based tests of independence for sparsely-sampled multivariate spatial point patterns. International Statistical Review, 51, 11–23.	Formal tests of interaction between marked types
Openshaw, S. & Taylor, P. J. (1979). A million or so correlated coefficients: three experiments on the modifiable areal unit problem. In Statistical Applications in the Spatial Sciences.	Classic empirical demonstration of MAUP

Object	Role
`geopandas.sjoin`	Point-in-polygon join used to aggregate events to neighbourhoods
`spatialstats.weights.Queen`, `spatialstats.esda.Moran`	Spatial weights and global autocorrelation (Part 2's core tools, previewed here)
`spatialstats.pointpattern.envelope(..., simulator=...)`	Custom null models for Monte Carlo tests, used for the random-labelling test in Exercise 4
Chapters 1.1–1.4	Every method and finding this chapter built on: CRS correction, windows, intensity, second-order functions, formal tests
Part 2	Spatial weights, Moran's I, and the Modifiable Areal Unit Problem in full
Part 5	Empirical Bayes and Bayesian smoothing of exactly the small-area rate instability previewed in Section 3