0.4 Datasets, Maps and a Reproducible Workflow

Part: 0: Setup and Foundations (Chapter 0.4 of 4, final chapter)

Dataset: Bundled data/ catalogue; Northeast diabetes lattice for choropleths

Focus: Data catalogue and hygiene, first maps with spatialstats.viz, loading recipe, pitfalls, exercises

Library: spatialstats.datasets, spatialstats.viz

Companion: ch00_setup_and_foundations.ipynb

---

Table of Contents

1. Overview

2. Python Setup

3. Datasets and the Data Catalogue

4. The Four "GW" Loaders

5. First Maps with spatialstats.viz

6. A Silent Fallback

7. Putting It Together: Loading Recipe

8. Common Pitfalls

9. Exercises

10. Solutions

11. Summary and References

---

1. Overview

This final foundations notebook closes the loop: know what is in data/, draw a careful first map, and adopt a loading recipe that catches CRS and typing problems before analysis begins.

Step	What we do
Catalogue	Inspect every bundled file; synthetic vs real
Hygiene	Key-matching and categorical typing traps
Maps	Quantiles, equal interval, Fisher–Jenks; GVF; counts vs rates
Recipe	crs_audit / load_study_area
Practice	Pitfalls table and exercises with solutions

Learning objectives (this notebook). Read the data catalogue; classify and map a variable; recognise misleading maps; apply a reproducible CRS audit before analysis.

---

2. Python Setup

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

3. Datasets and the Data Catalogue

What is in data/?

Instead of a hand-written list (which goes stale), the next cell inspects every file and builds the catalogue. The CRS column is what the file declares — as section 0.5 showed, declaring is not the same as being right.

def describe(path: Path, layer=None):
    ext = path.suffix.lower()
    if ext in {".shp", ".gpkg"}:
        g = gpd.read_file(path, layer=layer) if layer else gpd.read_file(path)
        crs = g.crs
        crs_txt = "none" if crs is None else (f"EPSG:{crs.to_epsg()}" if crs.to_epsg() else crs.name)
        kinds = "/".join(sorted(g.geom_type.unique()))
        cols = [c for c in g.columns if c != "geometry"]
        return dict(file=path.name + (f"[{layer}]" if layer else ""), kind="vector", rows=len(g), geometry=kinds, CRS=crs_txt, columns=len(cols),
                    fields=", ".join(cols[:6]) + (" …" if len(cols) > 6 else ""))
    if ext == ".csv":
        d = pd.read_csv(path)
        return dict(file=path.name, kind="table", rows=len(d), geometry="—", CRS="—", columns=d.shape[1],
                    fields=", ".join(d.columns[:6]) + (" …" if d.shape[1] > 6 else ""))

files = sorted(p for p in DATA.iterdir() if p.suffix.lower() in {".shp", ".csv", ".gpkg"})
records = []
for p in files:
    if p.suffix.lower() == ".gpkg":                       # a GeoPackage can hold several layers
        records += [describe(p, layer=name) for name in gpd.list_layers(p)["name"]]
    else:
        records.append(describe(p))
catalogue = pd.DataFrame(records).set_index("file")
catalogue

Things the catalogue makes visible:

• Several files exist twice — as a shapefile and a CSV (diabetes_*) — with different column names (Dia_pct vs Diabetes_per, x vs X).
• CRSs are not uniform across files, and one file (Crime_Incidents_2018.shp) has none. Harmonise before combining layers.
• cancer_areas.shp declares a custom Lambert conformal conic with no EPSG code; that is legal, but you cannot refer to it by number.

Synthetic datasets with known truth

Real data never tell you the right answer. For learning (and for testing a method) a synthetic dataset with known structure is invaluable. load_lattice(n) returns an $n\times n$ grid of square cells with the spatially structured attribute $$ z(x,y)=\sin\!\Big(\frac{2\pi x}{n}\Big)+\cos\!\Big(\frac{2\pi y}{n}\Big)+\varepsilon,\qquad \varepsilon\sim N(0,\,0.3^2), $$ plus a binary variable, $1[z>\operatorname{median}(z)]$. The smooth part is a known spatial trend, so a good autocorrelation statistic must detect it. load_points_csr(n) draws a complete spatial randomness pattern — the null model of Chapter 1.

from spatialstats.datasets import load_lattice, load_points_csr

lattice, zcol = load_lattice(n=12, seed=0)
pts_csr, _ = load_points_csr(n=150, seed=1)
print("load_lattice   ->", type(lattice).__name__, lattice.shape, "| columns:", [c for c in lattice.columns if c != "geometry"])
print("load_points_csr ->", type(pts_csr).__name__, pts_csr.shape)

fig, ax = plt.subplots(1, 3, figsize=(13.5, 4.2))
lattice.plot(ax=ax[0], column="z", cmap="RdBu_r", edgecolor="white", linewidth=0.4, legend=True, legend_kwds={"shrink": 0.7})
ax[0].set_title("load_lattice(): continuous z")
lattice.plot(ax=ax[1], column="binary", cmap="Greys", edgecolor="white", linewidth=0.4, categorical=True, legend=True)
ax[1].set_title("load_lattice(): binary  1[z > median]")
ax[2].scatter(pts_csr.geometry.x, pts_csr.geometry.y, s=14, color=C["red"])
ax[2].set(xlim=(0, 1), ylim=(0, 1), aspect="equal", title="load_points_csr(): CSR pattern")
for a in ax[:2]:
    a.set_axis_off()
plt.tight_layout(); plt.show()

4. The Four "GW" Loaders

load_georgia, load_londonhouse, load_beijing_pm25 and load_ncovid have names of real datasets, but they are simulated. Their docstrings say so, and the numbers show it.

from spatialstats import datasets as ds
from spatialstats.gwmodel import datasets as gw_ds      # the wrappers in `datasets` carry no docstring; the originals do

for name in ["load_georgia", "load_londonhouse", "load_beijing_pm25", "load_ncovid"]:
    print(f"{name:18s}: {inspect.getdoc(getattr(gw_ds, name)).splitlines()[0]}")

gdf_g, features, target = ds.load_georgia()
print(f"\nload_georgia returned {len(gdf_g)} points; longitude range {gdf_g.x.min():.2f} to {gdf_g.x.max():.2f}, "
      f"latitude {gdf_g.y.min():.2f} to {gdf_g.y.max():.2f}; features = {features}, target = {target!r}")
print("The locations are uniform random points in a box, not the 159 Georgia county centroids.")

Use them to try out an API; never quote results from them as if they were about real places. Real data live in data/.

Data hygiene: three traps in the bundled files

1. Key columns. County codes (FIPS) identify a county across tables, but their type varies from file to file.

us_csv = pd.read_csv(DATA / "diabetes_usa.csv")
us_shp = gpd.read_file(DATA / "US_COUNTY.shp")
cancer_areas = gpd.read_file(DATA / "cancer_areas.shp")
cancer_pts = pd.read_csv(DATA / "cancer_points.csv")

print(pd.DataFrame({"FIPS dtype": [us_csv.FIPS.dtype, us_shp.FIPS.dtype, cancer_areas.FIPS.dtype, cancer_pts.FIPS.dtype]},
                   index=["diabetes_usa.csv", "US_COUNTY.shp", "cancer_areas.shp", "cancer_points.csv"]).astype(str).to_string())

# numeric keys: pandas happily joins int64 with float64 ...
m1 = us_csv.merge(us_shp[["FIPS", "STATE"]], on="FIPS")
m2 = cancer_areas.merge(cancer_pts.groupby("FIPS", as_index=False)["POP10"].sum(), on="FIPS")
print(f"\nnumeric merge  csv x shp: {len(m1)} rows | areas x points: {len(m2)} rows  (both succeed)")

# ... but string keys built from floats do not match, and leading zeros were lost upstream
as_text = lambda s: s.astype(str)
v = cancer_pts.FIPS.iloc[0]
print("string keys    areas ∩ points:", len(set(as_text(cancer_areas.FIPS)) & set(as_text(cancer_pts.FIPS))),
      f"matching keys, although the numeric overlap is {len(set(cancer_areas.FIPS) & set(cancer_pts.FIPS))}: the same county is "
      f"'{int(v)}' in one file and '{v}' in the other")
print("lowest FIPS in the file:", us_csv.FIPS.min(), "-> the real code for Alabama's Autauga County is '01001'")

fips5 = lambda s: s.astype("int64").astype(str).str.zfill(5)         # the safe canonical key
print("canonical 5-digit string key:", fips5(us_csv.FIPS).iloc[:3].tolist(), "| unique in csv:", fips5(us_csv.FIPS).is_unique)

2. Categorical columns stored as text. Regression needs numbers; UrbRural holds the strings "Urban" / "Rural" and must be coded before it goes into a design matrix (Chapters 3 and 6).

3. Duplicate locations. Geocoded events pile up on street addresses; point-pattern statistics assume distinct locations (Chapter 1).

print("UrbRural values:", ne["UrbRural"].value_counts().to_dict(), "-> dummy:", "(ne.UrbRural == 'Urban').astype(float)")
ne["Urban"] = (ne["UrbRural"] == "Urban").astype(float)

xy_c = np.column_stack([crimes.geometry.x, crimes.geometry.y])
n_unique = len(np.unique(xy_c, axis=0))
display(Markdown(f"**Crime locations:** {len(xy_c):,} events but only {n_unique:,} distinct coordinates — "
                 f"{100 * (1 - n_unique / len(xy_c)):.0f}% of events share a location with an earlier event."))

5. First Maps with spatialstats.viz

Theory: turning numbers into classes

A choropleth colours each area by its value. With continuous data it is usually clearer to group values into $k$ classes with breakpoints $b_0<b_1<\dots<b_k$ and give each class one colour. How to choose the breaks is a modelling decision:

• Quantiles: $b_j=Q(j/k)$, the $j/k$ sample quantile. Every class holds about $n/k$ areas; the map is always "balanced" but breaks may split similar values.
• Equal interval: $b_j=x_{\min}+j\,\dfrac{x_{\max}-x_{\min}}{k}$. Easy to read; wasteful for skewed data (most areas fall in one or two classes).
• Natural breaks (Fisher–Jenks): choose breaks to minimise the within-class sum of squared deviations, $$ \text{SDCM}=\sum_{c=1}^{k}\;\sum_{i\in C_c}\big(x_i-\bar x_c\big)^2 . $$ The exact optimum is found by dynamic programming (Fisher, 1958; Jenks, 1967).

Classifications can be compared with the goodness of variance fit $$ \text{GVF}=1-\frac{\text{SDCM}}{\text{SDAM}},\qquad \text{SDAM}=\sum_{i=1}^{n}(x_i-\bar x)^2 , $$ the share of total variation explained by the classes ($0\le\text{GVF}\le1$). By construction Fisher–Jenks maximises GVF for a given $k$; GVF rises with $k$ with diminishing returns, so pick $k$ near the elbow — and since readers can rarely tell more than 5–7 shades apart reliably, $k\le 7$ is a sensible ceiling.

The package's choropleth(gdf, column, scheme, k, ...) accepts 'quantiles', 'equal_interval' and 'natural_breaks'.

from spatialstats.viz import choropleth

schemes = ["quantiles", "equal_interval", "natural_breaks"]
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for ax, scheme in zip(axes, schemes):
    choropleth(ne, "Dia_pct", scheme=scheme, k=5, ax=ax, cmap="viridis")
plt.tight_layout(); plt.show()

Same variable, three impressions. The colour bar of this function runs over the class index (0–4), not over data values, so to see what the classes mean we compute the breaks with mapclassify (the library the package delegates to), together with the histogram they cut.

import mapclassify

y = ne["Dia_pct"].to_numpy()
classifiers = {"quantiles": mapclassify.Quantiles(y, k=5), "equal_interval": mapclassify.EqualInterval(y, k=5),
               "natural_breaks": mapclassify.FisherJenks(y, k=5)}

fig, axes = plt.subplots(1, 3, figsize=(14, 3.6), sharey=True)
for ax, (name, cl) in zip(axes, classifiers.items()):
    ax.hist(y, bins=30, color="#cbd5e0")
    for b in cl.bins[:-1]:
        ax.axvline(b, color=C["red"], lw=1.4)
    ax.set_title(name); ax.set_xlabel("diabetes prevalence (%)")
axes[0].set_ylabel("counties")
axes[0].yaxis.set_major_locator(plt.MaxNLocator(integer=True))
plt.tight_layout(); plt.show()

breaks = pd.DataFrame({n: np.round(cl.bins, 2) for n, cl in classifiers.items()}, index=[f"class {i + 1} upper bound" for i in range(5)])
breaks.loc["counts per class"] = [str(cl.counts.tolist()) for cl in classifiers.values()]
breaks

def gvf(values, labels):
    sdcm = sum(((values[labels == c] - values[labels == c].mean()) ** 2).sum() for c in np.unique(labels))
    return 1 - sdcm / ((values - values.mean()) ** 2).sum()

ks = range(2, 10)
curve = {name: [gvf(y, getattr(mapclassify, cls)(y, k=k).yb) for k in ks]
         for name, cls in [("quantiles", "Quantiles"), ("equal_interval", "EqualInterval"), ("natural_breaks", "FisherJenks")]}
gvf_table = pd.DataFrame(curve, index=pd.Index(list(ks), name="k"))

fig, ax = plt.subplots(figsize=(6.5, 4))
for name, vals in curve.items():
    ax.plot(list(ks), vals, marker="o", label=name)
ax.axvspan(4.5, 7.5, color="#edf2f7", zorder=0)
ax.set(xlabel="number of classes k", ylabel="goodness of variance fit", title="GVF by scheme and number of classes")
ax.legend(); plt.tight_layout(); plt.show()

display(gvf_table.loc[[3, 5, 7]].round(3))
assert (gvf_table["natural_breaks"] >= gvf_table[["quantiles", "equal_interval"]].max(axis=1) - 1e-9).all()
display(Markdown("Fisher–Jenks has the highest GVF at every $k$, as theory requires (checked by an assertion in this cell)."))

6. A Silent Fallback

spatialstats.viz delegates classification to mapclassify when it is installed. If it is not, the package falls back to a simple quantile rule — for all three schemes that need it, including "natural breaks" — without any warning. We can prove this by hiding mapclassify from the import system for a moment:

from spatialstats.viz.choropleth import _classify

saved = sys.modules.get("mapclassify")
sys.modules["mapclassify"] = None                 # makes "import mapclassify" raise ImportError
try:
    fallback_nb = _classify(y, "natural_breaks", 5)
    fallback_q = _classify(y, "quantiles", 5)
finally:
    sys.modules["mapclassify"] = saved

with_mc_nb = _classify(y, "natural_breaks", 5)
print("without mapclassify: 'natural_breaks' identical to 'quantiles'? ->", np.array_equal(fallback_nb, fallback_q))
print("with mapclassify   : 'natural_breaks' identical to 'quantiles'? ->", np.array_equal(with_mc_nb, _classify(y, "quantiles", 5)))
print(f"GVF actually achieved: fallback {gvf(y, fallback_nb):.3f} vs real Fisher-Jenks {gvf(y, with_mc_nb):.3f}")

So a map titled "natural breaks" made in an environment without the viz extra is really a quantile map. Install pyspatialstats[viz] (it includes mapclassify) and check it as in section 0.1.

Maps can mislead: counts versus rates

An area with more people has more cases even if its risk is identical. A map of raw counts therefore mostly draws the population; a map of rates (cases per person) draws the risk. Compare:

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
choropleth(ne, "Dia_count", scheme="quantiles", k=5, ax=axes[0], cmap="magma")
axes[0].set_title("Dia_count — number of cases")
choropleth(ne, "POP_Total", scheme="quantiles", k=5, ax=axes[1], cmap="magma")
axes[1].set_title("POP_Total — population")
choropleth(ne, "Dia_pct", scheme="quantiles", k=5, ax=axes[2], cmap="viridis")
axes[2].set_title("Dia_pct — prevalence (rate)")
plt.tight_layout(); plt.show()

r_counts = np.corrcoef(ne["Dia_count"], ne["POP_Total"])[0, 1]
r_rate = np.corrcoef(ne["Dia_pct"], ne["POP_Total"])[0, 1]
display(Markdown(f"Correlation with population: cases **{r_counts:.3f}**, prevalence **{r_rate:.3f}**. "
                 "The count map is essentially a map of where people live. Chapter 5 shows a second problem with maps of rates: "
                 "rates from *small* populations are noisy, and the most extreme values tend to come from the smallest areas."))

First look, then measure

Look at the prevalence map once more. The high values are not scattered at random; they gather in some parts of the region. That impression — nearby areas look alike — is exactly what Tobler's first law of geography states ("everything is related to everything else, but near things are more related than distant things", Tobler, 1970). An impression is not evidence. Building the neighbour structure $\mathbf{W}$ and testing it with Moran's $I$ is the subject of Chapter 2.

7. Putting It Together: Loading Recipe

Everything in this chapter can be turned into a checklist that runs automatically. The function below audits a layer — CRS present? metric? appropriate for the location? geometries valid? duplicates? — and load_study_area wraps the audit around loading and reprojection, returning a SpatialData and a report. Use it (or something like it) as the first step of every analysis; it takes seconds and would have caught the Buffalo problem immediately.

def crs_audit(gdf, name="layer"):
    """Return a DataFrame of data-quality checks for a GeoDataFrame."""
    rows = []
    add = lambda check, status, detail: rows.append({"check": check, "status": status, "detail": detail})
    crs = gdf.crs
    if crs is None:
        add("CRS declared", "FAIL", "no CRS: cannot reproject and units are unknown")
    else:
        add("CRS declared", "ok", crs.name)
        unit = crs.axis_info[0].unit_name
        add("metric projected CRS", "ok" if crs.is_projected and unit == "metre" else "warn",
            f"{'geographic' if crs.is_geographic else 'projected'}, unit = {unit}")
        if crs.area_of_use is not None:
            pts = gdf.geometry.representative_point()
            if crs.is_geographic:
                lon, lat = pts.x.to_numpy(), pts.y.to_numpy()
            else:
                lon, lat = Transformer.from_crs(crs, 4326, always_xy=True).transform(pts.x.to_numpy(), pts.y.to_numpy())
            w, s, e, n = crs.area_of_use.bounds
            in_lon = (lon >= w) & (lon <= e) if w <= e else (lon >= w) | (lon <= e)
            frac = float(np.mean(in_lon & (lat >= s) & (lat <= n)))
            deg = lambda v, pos, neg: f"{abs(v):.0f}°{pos if v >= 0 else neg}"
            add("inside the CRS's area of use", "ok" if frac > 0.99 else "FAIL",
                f"{100 * frac:.0f}% of features; CRS is for {deg(w, 'E', 'W')} to {deg(e, 'E', 'W')}, "
                f"data at {abs(np.median(lon)):.1f}°{'E' if np.median(lon) >= 0 else 'W'}, {abs(np.median(lat)):.1f}°{'N' if np.median(lat) >= 0 else 'S'}")
    bad = int((~gdf.geometry.is_valid).sum()); empty = int(gdf.geometry.is_empty.sum()); miss = int(gdf.geometry.isna().sum())
    add("valid, non-empty geometries", "ok" if bad + empty + miss == 0 else "warn", f"{bad} invalid, {empty} empty, {miss} missing")
    types = gdf.geom_type.value_counts().to_dict()
    add("geometry types", "ok" if len(types) == 1 else "note", str(types))
    dup = int(gdf.geometry.to_wkb().duplicated().sum())
    add("repeated geometries", "ok" if dup == 0 else "note", f"{dup} of {len(gdf)} features repeat an earlier geometry")
    add("unique index", "ok" if gdf.index.is_unique else "warn", f"{int(gdf.index.duplicated().sum())} duplicated index labels")
    return pd.DataFrame(rows).set_index("check")


def load_study_area(source, y=None, X=None, target_crs=None):
    """Load, audit and (optionally) reproject a layer. Returns (SpatialData, audit DataFrame)."""
    gdf = source if isinstance(source, gpd.GeoDataFrame) else gpd.read_file(source)
    report = crs_audit(gdf)
    roles = ([y] if isinstance(y, str) else list(y or [])) + list(X or [])
    n_missing = {c: int(gdf[c].isna().sum()) for c in roles}
    report.loc["no missing values in roles"] = ("ok" if not any(n_missing.values()) else "warn",
                                                 str(n_missing) if roles else "no roles assigned")
    if target_crs is not None:
        if gdf.crs is None:
            report.loc["reprojection"] = ("FAIL", "skipped: layer has no CRS")
        else:
            gdf = to_crs(gdf, target_crs)
            report.loc["reprojection"] = ("ok", f"-> EPSG:{gdf.crs.to_epsg()}")
    return SpatialData(gdf, y=y, X=X), report


def show(report):
    colour = {"ok": "#c6f6d5", "note": "#fefcbf", "warn": "#feebc8", "FAIL": "#fed7d7"}
    return report.style.map(lambda v: f"background-color: {colour.get(v, '')}", subset=["status"])

_, rep_ne = load_study_area(DATA / "diabetes_northeast.shp", y="Dia_pct", X=["Obesity", "PhysInact", "SVI"])
display(Markdown("**diabetes_northeast.shp** — the well-behaved case"))
show(rep_ne)

_, rep_buf = load_study_area(DATA / "Crime_location_2018.shp")
display(Markdown("**Crime_location_2018.shp** — a valid CRS, wrong for the place"))
display(show(rep_buf))

_, rep_inc = load_study_area(DATA / "Crime_Incidents_2018.shp", target_crs=26917)
display(Markdown("**Crime_Incidents_2018.shp** — no CRS at all"))
show(rep_inc)

The recipe flags what a human eye cannot see on a map: the Buffalo file outside its CRS's zone, and the incident file without any CRS (and therefore un-reprojectable). Note too the "repeated geometries" line for crimes — the duplicate-location problem of section 0.7.

8. Common Pitfalls

Pitfall	Symptom	Guard
Geographic (lon/lat) coordinates used with a distance method	Bandwidths in degrees; distortion of tens of percent	is_projected check; ensure_projected(..., target_crs=...) with a suitable CRS
CRS valid but wrong for the location (§0.5)	Areas and distances biased with no error	Compare data location with crs.area_of_use; crs_audit
CRS missing or different between layers	Empty spatial joins; nonsense overlays	Check .crs of every layer; to_crs to one common CRS
Web Mercator (3857) used for measurement	Areas inflated by $\sec^2\varphi$	Equal-area / local metric CRS
(lat, lon) vs (lon, lat) order (§0.6)	Silently wrong haversine distances	Remember: haversine takes latitude first, GeoPandas is (x, y)
Centroid assumed = coords_array() (§0.4)	Km-scale location differences	Set the coords role explicitly
Per-axis coordinate normalisation (§0.6)	Distorted neighbourhoods	Scale with one common range
Dense $n\times n$ distance matrix on large $n$	Memory errors	k-d tree / sparse weights
Keys as floats or strings (§0.7)	Joins with zero matches; lost leading zeros	Canonical 5-digit string key
natural_breaks without mapclassify (§0.8)	Quantile map with the wrong title	Install the viz extra; verify
Mapping counts instead of rates (§0.8)	Map reproduces population	Use rates; small-area smoothing (Chapter 5)
Simulated "real-looking" datasets (§0.7)	Results reported as if about real places	Read docstrings; use data/ for real work

9. Exercises

Try each before opening the solutions cell that follows.

1. A third CRS. Build a SpatialData from crime_data_2018.csv using its long/lat columns (from_frame), project it to UTM 17N (EPSG:26917), and compute the average ratio of inter-point distances to those in the shapefile's declared CRS. What would ensure_projected() with its default target give instead, and why is it also unsuitable?
2. Scaling up. How much memory would a dense distance matrix for all 3,107 US counties need (diabetes_usa.csv, columns X, Y)? Compute the 8 nearest neighbours of each county instead and report how long it takes.
3. Skew. POP_Total is strongly right-skewed. Compare quantile and equal-interval classification of POP_Total for $k=5$: how many counties fall in each class, and what is each scheme's GVF? Does a $\log$ transform change the answer? Is the highest GVF always the best map?
4. Choosing a projection. The California PM2.5 stations were plotted in CONUS Albers (5070). Re-project them to California Albers (EPSG:3310) and measure how the distance between the two most distant stations differs. Which CRS would you choose for a kriging analysis of California, and why?
5. Which data type? For each — (i) wildfire ignition points in 2020, (ii) tract-level asthma rates, (iii) hourly ozone at 60 monitors, (iv) county unemployment — name the data type (§0.3), the quantity of interest, and the chapter that treats it.

10. Solutions

# --- Exercise 1 -----------------------------------------------------------------------------
sd_csv = SpatialData.from_frame(crime_tbl, x_col="long", y_col="lat", crs=4326)
sd_17 = sd_csv.to_crs(26917)
# from_frame() set coords=['long','lat'] and to_crs() does not touch those columns (the trap of section 0.4),
# so take the coordinates from the GEOMETRY of the reprojected layer:
xy_of = lambda obj: np.column_stack([obj.gdf.geometry.x, obj.gdf.geometry.y])
c17 = xy_of(sd_17)
c11 = np.column_stack([crimes.geometry.x, crimes.geometry.y])              # declared UTM 11N
ii, jj = pairs[keep, 0], pairs[keep, 1]                                  # the pairs used in section 0.5 (> 100 m apart)
ratio = (np.linalg.norm(c11[ii] - c11[jj], axis=1) / np.linalg.norm(c17[ii] - c17[jj], axis=1)).mean()
c3857 = xy_of(sd_csv.ensure_projected())
ratio_3857 = (np.linalg.norm(c3857[ii] - c3857[jj], axis=1) / np.linalg.norm(c17[ii] - c17[jj], axis=1)).mean()
print(f"Ex.1  distances in the declared file are {ratio:.3f}x those in UTM 17N;")
print(f"      Web Mercator (the ensure_projected default) gives {ratio_3857:.3f}x  (theory: sec(42.9°) = {1 / np.cos(np.radians(42.9)):.3f}).")

# --- Exercise 2 -----------------------------------------------------------------------------
us_xy = us_csv[["X", "Y"]].to_numpy()
n = len(us_xy)
t0 = time.perf_counter(); dist8, idx8 = cKDTree(us_xy).query(us_xy, k=9); t8 = time.perf_counter() - t0
print(f"\nEx.2  dense matrix for n = {n:,}: {8 * n * n / 1e6:,.0f} MB; 8-NN by k-d tree: {t8 * 1000:.1f} ms, "
      f"median 8th-neighbour distance {np.median(dist8[:, -1]) / 1000:.0f} km")

# --- Exercise 3 -----------------------------------------------------------------------------
pop = ne["POP_Total"].to_numpy()
rows = []
for label, v in [("POP_Total", pop), ("log(POP_Total)", np.log(pop))]:
    for scheme, cls in [("quantiles", mapclassify.Quantiles), ("equal_interval", mapclassify.EqualInterval)]:
        c = cls(v, k=5)
        rows.append({"variable": label, "scheme": scheme, "class counts": str(c.counts.tolist()), "GVF": gvf(v, c.yb)})
ex3 = pd.DataFrame(rows).set_index(["variable", "scheme"])
print("\nEx.3")
display(ex3.round(3))
top = ne["POP_Total"].sort_values(ascending=False)
ei = mapclassify.EqualInterval(pop, k=5)
display(Markdown(
    f"On the raw, right-skewed population, equal-interval reaches a *higher* GVF than quantiles ({ex3.loc[('POP_Total', 'equal_interval'), 'GVF']:.2f} vs "
    f"{ex3.loc[('POP_Total', 'quantiles'), 'GVF']:.2f}) because a handful of huge counties dominate the total variance — yet it puts {ei.counts[0]} of {len(pop)} counties "
    f"({100 * ei.counts[0] / len(pop):.0f}%) in the lowest class, so the map shows almost nothing. Quantiles spread the counties evenly but mix very different sizes within the top class. "
    f"After a log transform both schemes give balanced classes and a high GVF (≥ {ex3.loc['log(POP_Total)', 'GVF'].min():.2f}). "
    "**GVF is a useful yardstick, not the whole story: look at the class counts too.**"))

# --- Exercise 4 -----------------------------------------------------------------------------
sites_5070 = sites
sites_3310 = sites.to_crs(3310)
D_a = euclidean_distance_matrix(np.column_stack([sites_5070.geometry.x, sites_5070.geometry.y])) / 1000
D_b = euclidean_distance_matrix(np.column_stack([sites_3310.geometry.x, sites_3310.geometry.y])) / 1000
i, j = np.unravel_index(D_b.argmax(), D_b.shape)
ll_sites = sites.to_crs(4326)
d_true = geod.inv(ll_sites.geometry.x.iloc[i], ll_sites.geometry.y.iloc[i], ll_sites.geometry.x.iloc[j], ll_sites.geometry.y.iloc[j])[2] / 1000
print(f"Ex.4  most distant station pair: geodesic {d_true:.1f} km | CONUS Albers {D_a[i, j]:.1f} km ({100 * (D_a[i, j] / d_true - 1):+.2f}%) "
      f"| California Albers {D_b[i, j]:.1f} km ({100 * (D_b[i, j] / d_true - 1):+.2f}%)")
print("      A regional CRS (3310) has smaller distortion because its standard parallels are chosen for California;")
print("      for kriging, where variogram distances matter, use it (or the local UTM zone).")

Exercise 5 — answers. (i) Wildfire ignition points: a point pattern, quantity of interest the intensity $\lambda(\mathbf{s})$ (Chapter 1). (ii) Tract-level asthma rates: lattice data, spatial dependence through $\mathbf{W}$; small tracts need smoothing (Chapters 2, 5). (iii) Hourly ozone at 60 monitors: geostatistical (a field $Z(\mathbf{s},t)$ observed at stations; interpolate in Chapter 7, add time in Chapter 9). (iv) County unemployment: lattice data (Chapters 2–4, 6).

11. Summary

• Environment. Extras keep the core small; a missing extra raises an ImportError naming what to install, only when the feature is used. Record versions with your results.
• Structure. core is the bottom layer; the analysis modules share a SpatialResult interface (summary, stats, to_frame, plot).
• Data types. Fields $Z(\mathbf{s})$, lattices $\{Z_i\}$ with weights $\mathbf{W}$, and point patterns with intensity $\lambda(\mathbf{s})$ call for different models; the areal unit is part of the data (MAUP).
• CRS. Every projection distorts: Web Mercator by $\sec\varphi$ in distance and $\sec^2\varphi$ in area. Use a projected, metric CRS suited to the region for anything involving distance or area — and verify that the declared CRS fits the location.
• Distances. Euclidean (projected) or haversine (lat first!); matrices cost $8n^2$ bytes, so prefer k-d trees for large $n$; rescale coordinates with a single common range.
• Maps. Classification is a choice (GVF quantifies it); confirm mapclassify is installed; map rates, not counts.
• Habit. Audit every layer before you analyse it.

Where next

You have finished Part 0. Continue with:

• 0.1 · 0.2 · 0.3 — earlier foundations (if you jumped ahead)
• 1.1 Fundamentals of Point Patterns — Part 1
• Or Chapter 2 — Spatial autocorrelation turns the visual impression from §0.8 into a test: it builds spatial weights matrices (queen, rook, distance band, $k$-nearest neighbours), then computes and interprets global and local Moran's $I$. Chapter 1 (point patterns) can be read independently.

References

• Brewer, C. A. (2015). Designing Better Maps: A Guide for GIS Users (2nd ed.). Esri Press.
• Cressie, N. (1993). Statistics for Spatial Data. Wiley.
• Fisher, W. D. (1958). On grouping for maximum homogeneity. Journal of the American Statistical Association, 53, 789–798.
• Jenks, G. F. (1967). The data model concept in statistical mapping. International Yearbook of Cartography, 7, 186–190.
• Karney, C. F. F. (2013). Algorithms for geodesics. Journal of Geodesy, 87, 43–55.
• Openshaw, S. (1984). The Modifiable Areal Unit Problem. CATMOG 38, Geo Books.
• Sinnott, R. W. (1984). Virtues of the haversine. Sky and Telescope, 68(2), 159.
• Snyder, J. P. (1987). Map Projections — A Working Manual. USGS Professional Paper 1395.
• Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Economic Geography, 46, 234–240.
• Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2nd ed.). Springer.