import ast, importlib, importlib.util, inspect, platform, sys, time, 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:
    """Walk up from ``start`` until the folder containing ``pyproject.toml`` and ``data/`` is found."""
    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:                       # package not pip-installed: fall back to the source tree
    sys.path.insert(0, str(ROOT))
    import spatialstats

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"}

print(f"working directory : {Path.cwd().name}/")
print(f"repository data   : {DATA.relative_to(ROOT)}/  ({len(list(DATA.iterdir()))} files)")
print(f"spatialstats      : version {spatialstats.__version__}")

# Lattice layer used in several later demos
ne = gpd.read_file(DATA / "diabetes_northeast.shp")
print(f"diabetes_northeast: {len(ne)} counties · CRS={ne.crs}")

from pyproj import Geod
from spatialstats.core import validate_crs, is_geographic, is_projected, to_crs, ensure_projected

inventory = []
for epsg in [4326, 3857, 5070, 3310, 26917, 3718]:
    crs = validate_crs(f"EPSG:{epsg}")
    inventory.append({
        "EPSG": epsg, "name": crs.name,
        "kind": "geographic" if is_geographic(crs) else "projected" if is_projected(crs) else "other",
        "unit": crs.axis_info[0].unit_name,
        "designed for": crs.area_of_use.name[:58] + ("…" if len(crs.area_of_use.name) > 58 else ""),
    })
pd.DataFrame(inventory).set_index("EPSG")

geod = Geod(ellps="WGS84")
ne_ll = ne.to_crs(4326)
geodesic_km2 = np.array([abs(geod.geometry_area_perimeter(g)[0]) for g in ne_ll.geometry]) / 1e6

candidates = {"Web Mercator (3857)": 3857, "CONUS Albers, equal-area (5070)": 5070,
              "UTM 18N (32618)": 32618, "Global equal-area cylindrical (6933)": 6933}
planar_km2 = {name: ne.to_crs(code).area.to_numpy() / 1e6 for name, code in candidates.items()}

ratio_table = pd.DataFrame({
    "total area (km²)": {name: a.sum() for name, a in planar_km2.items()},
    "ratio to geodesic": {name: a.sum() / geodesic_km2.sum() for name, a in planar_km2.items()},
    "worst county (ratio)": {name: (a / geodesic_km2)[np.argmax(np.abs(a / geodesic_km2 - 1))]
                             for name, a in planar_km2.items()},
})
ratio_table.loc["geodesic (truth)"] = [geodesic_km2.sum(), 1.0, 1.0]
display(ratio_table.round(4))

lat = ne_ll.geometry.representative_point().y.to_numpy()
sec2 = 1 / np.cos(np.radians(lat)) ** 2
merc_ratio = planar_km2["Web Mercator (3857)"] / geodesic_km2

fig, ax = plt.subplots(1, 2, figsize=(11.5, 4.4))
ax[0].scatter(lat, merc_ratio, s=16, color=C["blue"], label="counties (measured)")
grid_lat = np.linspace(lat.min(), lat.max(), 100)
ax[0].plot(grid_lat, 1 / np.cos(np.radians(grid_lat)) ** 2, color=C["red"], lw=2, label=r"theory: $\sec^2\varphi$")
ax[0].axhline(1, color="k", lw=0.8)
ax[0].set(xlabel="latitude (°N)", ylabel="planar area / geodesic area", title="Web Mercator inflates area by sec²φ")
ax[0].legend()
for name, a in planar_km2.items():
    ax[1].scatter(lat, a / geodesic_km2, s=10, alpha=0.7, label=name)
ax[1].axhline(1, color="k", lw=0.8)
ax[1].set(xlabel="latitude (°N)", ylabel="planar area / geodesic area", title="Area distortion by CRS", ylim=(0.98, 2.1))
ax[1].legend(fontsize=8, loc="upper left")
plt.tight_layout(); plt.show()

r = np.corrcoef(merc_ratio, sec2)[0, 1]
display(Markdown(
    f"**Take-away.** Correlation between the measured Mercator area ratio and $\\sec^2\\varphi$ is **{r:.4f}**; "
    f"summed over the region Web Mercator reports **{ratio_table.loc['Web Mercator (3857)', 'ratio to geodesic']:.2f}×** the true area. "
    f"An equal-area CRS is within {abs(ratio_table.loc['CONUS Albers, equal-area (5070)', 'ratio to geodesic'] - 1) * 100:.4f}% of the truth in total, "
    f"and within {abs(ratio_table.loc['CONUS Albers, equal-area (5070)', 'worst county (ratio)'] - 1) * 100:.3f}% in the worst county."))

from pyproj import Transformer

albany, boston = (-73.7562, 42.6526), (-71.0589, 42.3601)              # (lon, lat)
d_true = geod.inv(*albany, *boston)[2] / 1000                          # geodesic on the ellipsoid, km

def euclid_in(epsg):
    tr = Transformer.from_crs(4326, epsg, always_xy=True)
    (x1, y1), (x2, y2) = tr.transform(*albany), tr.transform(*boston)
    return np.hypot(x2 - x1, y2 - y1) / 1000

rows = {
    "geodesic on WGS84 ellipsoid (truth)": d_true,
    "Euclidean on raw degrees × 111.32 km": np.hypot((boston[0] - albany[0]), (boston[1] - albany[1])) * 111.32,
    "Euclidean in Web Mercator (3857)": euclid_in(3857),
    "Euclidean in CONUS Albers (5070)": euclid_in(5070),
    "Euclidean in UTM 18N (32618)": euclid_in(32618),
}
dist_table = pd.DataFrame({"distance (km)": rows})
dist_table["error vs truth (%)"] = 100 * (dist_table["distance (km)"] / d_true - 1)
dist_table.round(2)

# ensure_projected: three behaviours worth knowing
geo = ne.to_crs(4326)
print("lon/lat -> default target        :", ensure_projected(geo).crs.to_epsg())
print("lon/lat -> explicit target_crs   :", ensure_projected(geo, target_crs=5070).crs.to_epsg())
print("already projected -> returned as-is :", ensure_projected(ne).crs.to_epsg())

nocrs = gpd.GeoDataFrame(ne.drop(columns="geometry"), geometry=list(ne.geometry))   # same data, CRS forgotten
print("CRS of the copy:", nocrs.crs)
try:
    ensure_projected(nocrs)
except ValueError as err:
    print("no CRS -> ValueError                :", err)

city = gpd.read_file(DATA / "buffalo_city.shp")
crimes = gpd.read_file(DATA / "Crime_location_2018.shp")
crime_tbl = pd.read_csv(DATA / "crime_data_2018.csv")            # the same events with lon/lat columns

declared = city.crs
west, south, east, north = declared.area_of_use.bounds
ll = city.to_crs(4326)
union = ll.geometry.union_all() if hasattr(ll.geometry, "union_all") else ll.geometry.unary_union
lon0, lat0 = union.centroid.x, union.centroid.y

print(f"declared CRS      : {declared.name}")
deg = lambda v, pos, neg: f"{abs(v):.1f}°{pos if v >= 0 else neg}"
print(f"designed for      : longitudes {deg(west, 'E', 'W')} to {deg(east, 'E', 'W')}, latitudes {deg(south, 'N', 'S')} to {deg(north, 'N', 'S')}")
print(f"data are located  : longitude {deg(lon0, 'E', 'W')}, latitude {deg(lat0, 'N', 'S')}  (Buffalo, New York)")
print(f"inside the CRS's area of use?  {west <= lon0 <= east and south <= lat0 <= north}")

from scipy.spatial import cKDTree

# 1. Distances: pairs of crimes 100 m - 3 km apart in the file; geodesic distance from their lon/lat
xy = np.column_stack([crimes.geometry.x, crimes.geometry.y])
lonlat = np.column_stack(Transformer.from_crs(declared, 4326, always_xy=True).transform(xy[:, 0], xy[:, 1]))
rng = np.random.default_rng(0)
sample = rng.choice(len(xy), 4000, replace=False)
pairs = cKDTree(xy[sample]).query_pairs(3000, output_type="ndarray")
pairs = sample[pairs[rng.choice(len(pairs), 3000, replace=False)]]
d_file = np.linalg.norm(xy[pairs[:, 0]] - xy[pairs[:, 1]], axis=1)
d_geo = geod.inv(lonlat[pairs[:, 0], 0], lonlat[pairs[:, 0], 1], lonlat[pairs[:, 1], 0], lonlat[pairs[:, 1], 1])[2]
keep = d_geo > 100
k_measured = (d_file[keep] / d_geo[keep]).mean()

# 2. Theory
B = np.cos(np.radians(lat0)) * np.sin(np.radians(lon0 + 117))
k_theory = 0.9996 / np.sqrt(1 - B ** 2)

# 3. Areas: planar area in the file vs geodesic area vs the layer's own 'sqmiles' attribute
area_file = city.area.sum() / 1e6
area_geo = sum(abs(geod.geometry_area_perimeter(g)[0]) for g in ll.geometry) / 1e6
area_attr = city["sqmiles"].sum() * 2.589988

# 3b. Rotation: direction of the same pair-vectors in the file (UTM 11N) and in the correct zone (UTM 17N)
xy17 = np.column_stack(Transformer.from_crs(4326, 26917, always_xy=True).transform(lonlat[:, 0], lonlat[:, 1]))
pi_, pj_ = pairs[keep, 0], pairs[keep, 1]
ang = lambda P: np.arctan2(P[pj_, 1] - P[pi_, 1], P[pj_, 0] - P[pi_, 0])
dtheta = np.degrees((ang(xy) - ang(xy17) + np.pi) % (2 * np.pi) - np.pi)
rot_measured = np.median(dtheta)
gamma = lambda lam0: np.degrees(np.arctan(np.tan(np.radians(lon0 - lam0)) * np.sin(np.radians(lat0))))
rot_theory = gamma(-117) - gamma(-81)

# 4. Do the file's lon/lat (via the declared CRS) agree with the lon/lat columns of the CSV?
agree = np.abs(lonlat - crime_tbl[["long", "lat"]].to_numpy()).max()

audit = pd.DataFrame({
    "quantity": ["linear scale factor k", "area (km²)", "rotation vs. correct zone"],
    "as stored in the file": [f"{k_measured:.4f} (measured)", f"{area_file:.1f}", f"{abs(rot_measured):.1f}° (measured)"],
    "theory / truth": [f"{k_theory:.4f} (formula)", f"{area_geo:.1f} geodesic; {area_attr:.1f} from 'sqmiles' attribute",
                       f"{abs(rot_theory):.1f}° (γ formula, difference of the two zones)"],
    "distortion": [f"distances × {k_measured:.3f}", f"areas × {area_file / area_geo:.3f}", "map turned by the grid convergence"],
}).set_index("quantity")
display(audit)
print(f"lon/lat recovered through the declared CRS agree with the CSV's lon/lat to within {agree:.4f} degrees (about {agree * 111_000:.0f} m):")
print("-> the coordinates locate Buffalo correctly; what is wrong is the zone the numbers are expressed in.")

city_fixed = city.to_crs(4326).to_crs(26917)
crimes_fixed = crimes.to_crs(4326).to_crs(26917)

xy_f = np.column_stack([crimes_fixed.geometry.x, crimes_fixed.geometry.y])
k_fixed = (np.linalg.norm(xy_f[pairs[:, 0]] - xy_f[pairs[:, 1]], axis=1)[keep] / d_geo[keep]).mean()

fix = pd.DataFrame({
    "declared UTM 11N": [f"{k_measured:.4f}", f"{area_file:.1f}"],
    "re-projected to UTM 17N (EPSG:26917)": [f"{k_fixed:.4f}", f"{city_fixed.area.sum() / 1e6:.1f}"],
    "truth": ["1.0000", f"{area_attr:.1f}"],
}, index=["distance scale factor", "area (km²)"])
display(fix)

fig, ax = plt.subplots(1, 2, figsize=(11, 4))
km = lambda v, _: f"{v / 1000:,.0f}"
for a, (g, ttl) in zip(ax, [(city, f"declared CRS (UTM 11N): {area_file:.0f} km²"), (city_fixed, f"UTM 17N: {city_fixed.area.sum() / 1e6:.0f} km²")]):
    g.boundary.plot(ax=a, color="k", lw=0.6); a.set_title(ttl); a.set_aspect("equal")
    a.set_xlabel("easting (km)"); a.set_ylabel("northing (km)")
    a.xaxis.set_major_formatter(km); a.yaxis.set_major_formatter(km); a.tick_params(labelsize=8)
    a.xaxis.set_major_locator(plt.MaxNLocator(4))
plt.tight_layout(); plt.show()

counties = ne                                            # CRS: EPSG:5070 (Albers, metres)
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    joined = gpd.sjoin(crimes.head(200), counties)       # crimes: 'UTM 11N'  vs  counties: 'Albers'
print(f"spatial join of 200 crimes with counties in a different CRS -> {len(joined)} matches")
print("GeoPandas warned:", str(caught[0].message).splitlines()[0])

incidents = gpd.read_file(DATA / "Crime_Incidents_2018.shp")
print("\nCrime_Incidents_2018.shp has CRS:", incidents.crs)
print("its coordinates look like:", incidents.geometry.x.iloc[:2].round(3).tolist(), incidents.geometry.y.iloc[:2].round(3).tolist(),
      "-> plausibly degrees; the file also has 'longitude' / 'latitude' columns")
try:
    to_crs(incidents, 26917)
except ValueError as err:
    print("core.to_crs on a file with no CRS -> ValueError:", err)

same = np.allclose(incidents.geometry.x.astype(float), pd.to_numeric(incidents["longitude"], errors="coerce"), atol=1e-6, equal_nan=True)
if same:
    incidents = incidents.set_crs(4326)                   # only after checking, never blindly
    print("after verifying the coordinates equal the lon/lat columns, set_crs(4326) ->", incidents.crs.to_epsg())

from spatialstats.core import (coords_from_geometry, euclidean_distance_matrix, euclidean_distances,
                               haversine_distance_matrix, pairwise_distances, normalize_coords)

# Length of one degree of longitude at various latitudes: haversine vs the closed form (pi/180) R cos(phi)
R = 6371.0
lats = np.array([0, 20, 40, 45, 60, 80])
via_haversine = np.array([haversine_distance_matrix(np.array([[la, 0.0], [la, 1.0]]))[0, 1] for la in lats])
closed_form = np.pi / 180 * R * np.cos(np.radians(lats))
pd.DataFrame({"haversine (km)": via_haversine, "(π/180)·R·cos φ (km)": closed_form}, index=pd.Index(lats, name="latitude (°)")).round(3)

cities = {"New York": (40.7128, -74.0060), "Boston": (42.3601, -71.0589), "Atlanta": (33.7490, -84.3880),
          "Miami": (25.7617, -80.1918), "Seattle": (47.6062, -122.3321), "Denver": (39.7392, -104.9903)}
pairs_ = [("New York", "Boston"), ("New York", "Atlanta"), ("Atlanta", "Denver"), ("Miami", "Seattle")]
to_albers = Transformer.from_crs(4326, 5070, always_xy=True)

rows = []
for a, b in pairs_:
    (la1, lo1), (la2, lo2) = cities[a], cities[b]
    truth = geod.inv(lo1, la1, lo2, la2)[2] / 1000
    hav = haversine_distance_matrix(np.array([[la1, lo1], [la2, lo2]]))[0, 1]
    x1, y1 = to_albers.transform(lo1, la1); x2, y2 = to_albers.transform(lo2, la2)
    alb = np.hypot(x2 - x1, y2 - y1) / 1000
    swapped = haversine_distance_matrix(np.array([[lo1, la1], [lo2, la2]]))[0, 1]     # (lon, lat) passed by mistake
    rows.append({"pair": f"{a} – {b}", "geodesic (km)": truth, "haversine (km)": hav,
                 "haversine error (%)": 100 * (hav / truth - 1), "Albers Euclid. error (%)": 100 * (alb / truth - 1),
                 "lon/lat swapped error (%)": 100 * (swapped / truth - 1)})
pd.DataFrame(rows).set_index("pair").round(2)

xy_proj = coords_from_geometry(ne)                                                    # (x, y) in metres: one point per county
# transform *those same points* to lon/lat (a representative point is not invariant under reprojection,
# so calling coords_from_geometry() again on the lon/lat layer would pick different points)
lon, lat = Transformer.from_crs(ne.crs, 4326, always_xy=True).transform(xy_proj[:, 0], xy_proj[:, 1])
lat_lon = np.column_stack([lat, lon])                                                 # haversine wants (lat, lon)

D_eu = euclidean_distance_matrix(xy_proj) / 1000                                      # km
D_hv = pairwise_distances(lat_lon, metric="haversine")                                # km
iu = np.triu_indices(len(ne), k=1)
rel = 100 * (D_eu[iu] / D_hv[iu] - 1)

fig, ax = plt.subplots(1, 2, figsize=(11.5, 4.3))
ax[0].scatter(D_hv[iu], D_eu[iu], s=2, alpha=0.15, color=C["blue"])
mx = D_hv[iu].max()
ax[0].plot([0, mx], [0, mx], color=C["red"], lw=1.2, label="y = x")
ax[0].set(xlabel="haversine distance (km)", ylabel="Euclidean in Albers (km)", title=f"{len(iu[0]):,} county pairs"); ax[0].legend()
ax[1].scatter(D_hv[iu], rel, s=2, alpha=0.15, color=C["purple"])
ax[1].axhline(0, color="k", lw=0.8)
ax[1].set(xlabel="haversine distance (km)", ylabel="Albers vs haversine (%)", title="Relative difference")
plt.tight_layout(); plt.show()

display(Markdown(f"The two methods agree on average within **{np.abs(rel).mean():.2f}%** (worst pair {np.abs(rel).max():.2f}%). "
                 f"Median and maximum distance between two counties: {np.median(D_hv[iu]):.0f} km, {D_hv[iu].max():.0f} km."))

rng = np.random.default_rng(0)
for name, D in [("Euclidean (Albers, km)", D_eu), ("haversine (km)", D_hv)]:
    i, j, k = rng.integers(0, len(D), size=(3, 20000))
    violation = D[i, k] - (D[i, j] + D[j, k])
    print(f"{name:24s} symmetric: {np.allclose(D, D.T)} | zero diagonal: {np.allclose(np.diag(D), 0)} | "
          f"min off-diagonal: {D[np.triu_indices(len(D), 1)].min():.3f} | worst triangle-inequality violation: {violation.max():.2e}")

from scipy.spatial import cKDTree

rng = np.random.default_rng(1)
rows = []
for n in [500, 1000, 2000, 3107, 4000]:
    pts = rng.random((n, 2)) * 1e5
    t0 = time.perf_counter(); D = euclidean_distance_matrix(pts); t_dense = time.perf_counter() - t0
    t0 = time.perf_counter(); cKDTree(pts).query(pts, k=9); t_tree = time.perf_counter() - t0
    rows.append({"n": n, "dense matrix (MB)": D.nbytes / 1e6, "dense time (s)": t_dense, "8-NN by k-d tree (s)": t_tree})
    del D
scale = pd.DataFrame(rows).set_index("n")
display(scale.round(3))
big = {"US counties": 3107, "Buffalo crimes": 14001, "CA PM2.5 station-days": 59566}
display(Markdown("**Dense-matrix footprint at real sizes:** " + "; ".join(f"{k} (n = {n:,}) → {8 * n * n / 1e9:.2f} GB" for k, n in big.items()) + "."))

def scaling_diagnostics(xy, region):
    rx, ry = np.ptp(xy[:, 0]), np.ptp(xy[:, 1])
    D0 = euclidean_distance_matrix(xy)
    iu_ = np.triu_indices(len(xy), k=1)
    A = D0.copy(); np.fill_diagonal(A, np.inf)
    rows = []
    for label, Z in [("normalize_coords (per axis)", normalize_coords(xy)),
                     ("single common range", (xy - xy.min(axis=0)) / max(rx, ry))]:
        Dn = euclidean_distance_matrix(Z)
        ratio = Dn[iu_] / D0[iu_]
        B = Dn.copy(); np.fill_diagonal(B, np.inf)
        rows.append({"region": region, "E-W / N-S extent": rx / ry, "scaling": label,
                     "corr with true distance": np.corrcoef(Dn[iu_], D0[iu_])[0, 1],
                     "distance ratio, max / min": ratio.max() / ratio.min(),
                     "same nearest neighbour (%)": 100 * np.mean(A.argmin(axis=1) == B.argmin(axis=1))})
    return rows

us_xy_ = pd.read_csv(DATA / "diabetes_usa.csv")[["X", "Y"]].to_numpy()
tab = pd.DataFrame(scaling_diagnostics(xy_proj, "Northeast (217 counties)") + scaling_diagnostics(us_xy_, "Contiguous US (3,107 counties)"))
tab.set_index(["region", "scaling"]).round(3)

Step	What we do
Theory	Scale factors, Mercator $\sec^2\varphi$, grid convergence
Measure	Area and distance inflation under Web Mercator vs equal-area
Case study	Buffalo files tagged UTM 11N (wrong zone for NY)
Distances	Euclidean, haversine, pairwise; metric checks; matrix cost; normalisation trap

Function	Purpose	Input order
`coords_from_geometry(geom)`	$(x,y)$ per feature (representative point for polygons)	—
`euclidean_distance_matrix(c)`	$n\times n$ Euclidean distances	$(x,y)$
`euclidean_distances(a, b)`	rectangular $m\times n$ distances between two sets	$(x,y)$
`haversine_distance_matrix(ll)`	great-circle distances in km	(latitude, longitude)
`pairwise_distances(c, metric)`	dispatcher: `euclidean`, `haversine`, `network`, or any SciPy metric	depends on metric
`normalize_coords(c)`	rescale each axis to $[0,1]$	$(x,y)$

0.3 Coordinate Reference Systems and Distances¶

Table of Contents¶

1. Overview¶

2. Python Setup¶

3. Coordinate Reference Systems¶

Theory: from an ellipsoid to a flat map¶

4. How Much Does the CRS Matter?¶

5. Distances Suffer Too¶

6. Case Study: Valid but Wrong CRS¶

Two more CRS failure modes¶

7. Distances¶

Theory: what is a distance?¶

The functions¶

8. Five Ways to Measure the Same Distance¶

Distance matrices for real data¶

9. Metric Axioms¶

10. Cost of an $n\times n$ Matrix¶

11. Normalising Coordinates¶

12. Summary and Where Next¶