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, clark_evans, g_function, f_function, j_function,
    ripley_k, ripley_l, pair_correlation,
)

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()

rng = np.random.default_rng(3)
sub_idx = rng.choice(pp_dedup.n, 2500, replace=False)
pp_sub = pp_dedup.subset(sub_idx)

print(f"full pattern           : n = {pp_all.n:,}")
print(f"deduplicated           : n = {pp_dedup.n:,}")
print(f"working subsample       : n = {pp_sub.n:,}  (seed=3, fixed for reproducibility)")
print(f"window short side        : {pp_sub.window.short_side / 1000:.2f} km")

r_range = np.linspace(0, 1500, 60)   # metres; well under a quarter of the window's short side
print(f"distance range used below: 0 to {r_range[-1]:.0f} m, {len(r_range)} points")

full pattern           : n = 13,948
deduplicated           : n = 8,903
working subsample       : n = 2,500  (seed=3, fixed for reproducibility)
window short side        : 9.92 km
distance range used below: 0 to 1500 m, 60 points

ce = clark_evans(pp_sub, correction="donnelly")
ce_naive = clark_evans(pp_sub, correction="none")

print(ce.summary())
p_str = "< 1e-10" if ce["p_value"] < 1e-10 else f"= {ce['p_value']:.2e}"
display(Markdown(
    f"**R = {ce['R']:.3f}** (Donnelly-corrected): the mean nearest-neighbour distance is only "
    f"{100 * ce['R']:.0f}% of what CSR predicts. z = {ce['z_score']:.1f}, p {p_str} — clustering, "
    f"overwhelmingly. The uncorrected version gives R = {ce_naive['R']:.3f}: close, because the window is "
    "large relative to typical nearest-neighbour distances here, but the boundary correction is still the "
    "safer default."
))

Clark-Evans test
================
  R                    0.786722
  mean_nn              82.2063
  expected_nn          104.492
  z_score              -19.7006
  p_value              2.13018e-86
  p_clustered          1.06509e-86
  p_regular            1
  correction           donnelly
  n                    2500

g = g_function(pp_sub, r_range)
f = f_function(pp_sub, r_range)
j = j_function(pp_sub, r_range)

fig, axes = plt.subplots(1, 3, figsize=(15, 4.6))
for ax, fn, label in zip(axes, [g, f, j], ["G(r)", "F(r)", "J(r)"]):
    fn.plot(ax=ax, color=C["blue"], lw=2)
    ax.set_title(label)
axes[0].set_ylabel("probability"); axes[2].set_ylabel("J(r)")
plt.suptitle("Nearest-neighbour distribution functions (Kaplan–Meier correction)", y=1.02)
plt.tight_layout(); plt.show()

# Both G(r) and F(r) are bounded above by 1, and the expected CSR nearest-neighbour distance here is only
# about 103 m (lambda = n/area) -- so by r ~= 250 m the *theoretical* CSR curve has already climbed to
# roughly 0.99, leaving almost no room for the comparison to be informative. We therefore look at the
# short-range window where the curves have not yet saturated.
short_range = (r_range > 0) & (r_range <= 228)
above_G = np.mean(g.estimate[short_range] > g.theoretical[short_range])
above_F = np.mean(f.estimate[short_range] > f.theoretical[short_range])   # clustering pulls F BELOW theo, so this should be low

# J(r) becomes undefined once F(r) reaches (numerically) 1: with a finite grid of test locations,
# every grid point eventually lies within r of *some* event, so 1 - F(r) -> 0. Restrict the summary
# to the r where J is still finite, and use np.isfinite explicitly: a NaN < 1 comparison silently
# evaluates to False in NumPy, which would understate "share below 1" if left unmasked.
j_finite = np.isfinite(j.estimate) & (r_range >= r_range[5])
below_J = np.mean(j.estimate[j_finite] < 1)
r_last_finite = r_range[np.isfinite(j.estimate)][-1]

print(f"share of 0 < r ≤ 228 m with G(r) above the CSR curve : {above_G:.2f}  (clustering pushes G up)")
print(f"share of 0 < r ≤ 228 m with F(r) above the CSR curve : {above_F:.2f}  (clustering pulls F down, so this should be low)")
print(f"J(r) stays finite up to r = {r_last_finite:.0f} m in this sample "
      f"({j_finite.sum()} of {(r_range >= r_range[5]).sum()} points beyond r = {r_range[5]:.0f} m)")
print(f"share of those finite points with J(r) below 1                : {below_J:.2f}  (clustering pushes J below 1)")

share of 0 < r ≤ 228 m with G(r) above the CSR curve : 1.00  (clustering pushes G up)
share of 0 < r ≤ 228 m with F(r) above the CSR curve : 0.00  (clustering pulls F down, so this should be low)
J(r) stays finite up to r = 1119 m in this sample (40 of 55 points beyond r = 127 m)
share of those finite points with J(r) below 1                : 1.00  (clustering pushes J below 1)

k = ripley_k(pp_sub, r_range)
l = ripley_l(pp_sub, r_range)

fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].plot(k.r, k.estimate, color=C["blue"], lw=2, label="observed $K(r)$")
axes[0].plot(k.r, k.theoretical, "k--", lw=1, label=r"CSR: $\pi r^2$")
axes[0].set_xlabel("r (m)"); axes[0].set_ylabel("K(r)"); axes[0].legend(); axes[0].set_title("Ripley's K")

axes[1].plot(l.r, l.estimate - l.r, color=C["blue"], lw=2, label="observed $L(r) - r$")
axes[1].axhline(0, color="k", ls="--", lw=1, label="CSR")
axes[1].set_xlabel("r (m)"); axes[1].set_ylabel("L(r) − r"); axes[1].legend()
axes[1].set_title("Variance-stabilised L")
plt.tight_layout(); plt.show()

ratio_at = {int(rv): float(k.estimate[np.searchsorted(r_range, rv)] / k.theoretical[np.searchsorted(r_range, rv)])
           for rv in [150, 500, 1000, 1500]}
print("K(r) / (pi r^2) at selected distances:")
for rv, ratio in ratio_at.items():
    print(f"  r = {rv:5d} m : {ratio:.2f}× CSR")

K(r) / (pi r^2) at selected distances:
  r =   150 m : 1.97× CSR
  r =   500 m : 1.58× CSR
  r =  1000 m : 1.40× CSR
  r =  1500 m : 1.31× CSR

g_pc = pair_correlation(pp_sub, r_range[1:])   # g(r) undefined at r = 0

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(g_pc.r, g_pc.estimate, color=C["blue"], lw=2, label="observed $g(r)$")
ax.axhline(1.0, color="k", ls="--", lw=1, label="CSR: $g(r) = 1$")
ax.set_xlabel("r (m)"); ax.set_ylabel("g(r)")
ax.set_title(f"Pair correlation function (bandwidth = {g_pc['bandwidth']:.0f} m)")
ax.legend()
plt.tight_layout(); plt.show()

print(f"g(r) at r ≈ 50 m  : {g_pc.estimate[np.searchsorted(g_pc.r, 50)]:.2f}")
print(f"g(r) at r ≈ 300 m : {g_pc.estimate[np.searchsorted(g_pc.r, 300)]:.2f}")
print(f"g(r) at r ≈ 1000 m: {g_pc.estimate[np.searchsorted(g_pc.r, 1000)]:.2f}")
print(f"g(r) at r ≈ 1500 m: {g_pc.estimate[np.searchsorted(g_pc.r, 1500)]:.2f}")

g(r) at r ≈ 50 m  : 2.25
g(r) at r ≈ 300 m : 1.60
g(r) at r ≈ 1000 m: 1.30
g(r) at r ≈ 1500 m: 1.20

k_trans = ripley_k(pp_sub, r_range, correction="translation")
k_border = ripley_k(pp_sub, r_range, correction="border")
k_none = ripley_k(pp_sub, r_range, correction="none")

fig, ax = plt.subplots(figsize=(8, 5.5))
for res, label, color in [(k_trans, "translation (default)", C["blue"]),
                          (k_border, "border (reduced-sample)", C["green"]),
                          (k_none, "none (uncorrected)", C["orange"])]:
    ax.plot(res.r, res.estimate, color=color, lw=2, label=label)
ax.plot(r_range, np.pi * r_range ** 2, "k--", lw=1, label="CSR")
ax.set_xlabel("r (m)"); ax.set_ylabel("K(r)")
ax.set_title("Ripley's K under three edge corrections")
ax.legend()
plt.tight_layout(); plt.show()

with np.errstate(invalid="ignore"):
    rel_none = k_none.estimate / k_trans.estimate
    rel_border = np.where(np.isfinite(k_border.estimate), k_border.estimate / k_trans.estimate, np.nan)
print("uncorrected K as a fraction of the translation-corrected K, at increasing r:")
for i in [10, 30, 59]:
    print(f"  r = {r_range[i]:5.0f} m : none = {rel_none[i]:.3f}×   border = {rel_border[i]:.3f}×"
          + ("  (border undefined here)" if np.isnan(rel_border[i]) else ""))

# why border degrades at large r: how many anchors remain, and how much of the window is left after erosion?
b_dist = pp_sub.window.distance_to_boundary(pp_sub.coords)
r_big = r_range[59]
n_anchors = int((b_dist >= r_big).sum())
eroded_km2 = pp_sub.window.polygon.buffer(-r_big).area / 1e6
print(f"\nat r = {r_big:.0f} m: only {n_anchors} of {pp_sub.n:,} events "
      f"({100 * n_anchors / pp_sub.n:.1f}%) are far enough from the boundary to serve as a border-correction "
      f"anchor, and the window eroded by r has shrunk to {eroded_km2:.1f} km² "
      f"(from {pp_sub.window.area / 1e6:.1f} km² originally).")

uncorrected K as a fraction of the translation-corrected K, at increasing r:
  r =   254 m : none = 0.975×   border = 1.056×
  r =   763 m : none = 0.933×   border = 0.857×
  r =  1500 m : none = 0.875×   border = 0.260×

at r = 1500 m: only 48 of 2,500 events (1.9%) are far enough from the boundary to serve as a border-correction anchor, and the window eroded by r has shrunk to 3.8 km² (from 105.8 km² originally).

Step	What we do
Theory	Define second-order structure; introduce the four families of summary function
Nearest-neighbour	Clark–Evans index, and the G, F, J distribution functions
Ripley's K / L	The classical cumulative measures of clustering and inhibition
Pair correlation	A non-cumulative function that localises the scale of interaction
Edge correction	Compare translation, border, and no correction on the same data

Function	Compares	Cumulative?	Best for
Clark–Evans $R$	Mean nearest-neighbour distance	—	One-number summary
$G(r)$	Distribution of nearest-neighbour distances	Yes	Short-range clustering
$F(r)$	Distribution of distances from empty space to the nearest event	Yes	Detecting large gaps
$J(r) = \tfrac{1-G(r)}{1-F(r)}$	Combines G and F	Yes	Single diagnostic curve
$K(r)$	Expected number of further events within $r$ of a typical event	Yes	The classical standard
$L(r) = \sqrt{K(r)/\pi}$	Variance-stabilised $K$	Yes	Easier visual comparison to CSR
$g(r)$	Density of pairs at exactly distance $r$	No	Localising the scale of clustering

Correction	Idea	Cost
`'translation'` (default)	Weight each pair by $1/\gamma(\mathbf{v})$, the window's set covariance at the pair's displacement vector — exact for rectangles, a fast raster approximation for arbitrary polygons	Moderate
`'border'`	Only count events whose entire radius-$r$ neighbourhood lies inside $W$ (a "reduced sample")	Discards data at large $r$
`'none'`	No correction	Biased low, worse as $r$ grows

Reference	Notes
Ripley, B. D. (1981). Spatial Statistics. Wiley.	The original systematic treatment of K, L and edge correction
Baddeley, A., Rubak, E. & Turner, R. (2015). Spatial Point Patterns: Methodology and Applications with R. CRC Press.	Chapters 7–8: nearest-neighbour and second-order functions in depth, with the `spatstat` implementation this package's API is modelled on
Diggle, P. J. (2013). Statistical Analysis of Spatial and Spatio-Temporal Point Patterns (3rd ed.). CRC Press.	Concise treatment of G, F, J and K
Stoyan, D. & Stoyan, H. (1994). Fractals, Random Shapes and Point Fields. Wiley.	The pair correlation function in a broader stochastic-geometry context

Paper	Contribution
Clark, P. J. & Evans, F. C. (1954). Distance to nearest neighbor as a measure of spatial relationships. Ecology, 35, 445–453.	The Clark–Evans index
Donnelly, K. (1978). Simulations to determine the variance and edge effect of total nearest neighbour distance. In Simulation Methods in Archaeology.	The finite-window correction used by `clark_evans`
Ripley, B. D. (1976). The second-order analysis of stationary point processes. Journal of Applied Probability, 13, 255–266.	The K function
Van Lieshout, M. N. M. & Baddeley, A. J. (1996). A nonparametric measure of spatial interaction in point patterns. Statistica Neerlandica, 50, 344–361.	The J function

1.3 Second-order and Distance-based Summary Functions¶

Table of Contents¶

1. Overview¶

2. Second-order Structure: What These Functions Measure¶

3. Nearest-neighbour Methods¶

3.1 The Clark–Evans index¶

3.2 G and F: event and empty-space distributions¶

3.3 J: a combined diagnostic¶

4. Ripley's K, L and the Pair Correlation Function¶

4.1 Ripley's K function¶

4.2 The variance-stabilised L function¶

4.3 The pair correlation function g(r)¶

4.4 Edge corrections¶

5. Python Setup¶

6. Data: A Working Subsample¶

7. Clark–Evans Nearest-neighbour Test¶

8. G, F and J Functions¶

9. Ripley's K and L Functions¶

10. Pair Correlation Function¶

11. Comparing Edge Corrections¶

12. Summary and Conclusions¶

What we did¶

Practical takeaways¶

Next steps¶

13. Resources¶

Textbooks¶

Key papers¶

In this library¶

Object	Role
`spatialstats.pointpattern.clark_evans`	Nearest-neighbour index with Donnelly's finite-window correction
`spatialstats.pointpattern.g_function`, `f_function`, `j_function`	Nearest-event, empty-space and combined distribution functions
`spatialstats.pointpattern.ripley_k`, `ripley_l`	Cumulative second-order functions with `'translation'`, `'border'`, `'none'` corrections
`spatialstats.pointpattern.pair_correlation`	The non-cumulative pair correlation function $g(r)$
Chapter 1.1	`Window`, `PointPattern`, `deduplicate()`, the CRS correction reused here
Chapter 1.4	Formal Monte Carlo tests (`envelope`), calibration against known processes, `k_inhom`