Your First Python GIS Analysis: From Download to Map in One Script

On this page

Problem statement

Every tutorial teaches one operation. Read a shapefile. Do a spatial join. Make a choropleth. Each works in isolation, and then you try to answer a real question and discover that the gaps between the steps are where the work is.

The question here is concrete: which wards have the highest density of reported incidents, and where should attention go?

The data is real: Leeds' 33 electoral wards with their Census 2021 population, and every street-level crime West Yorkshire Police published for June 2026 around Leeds (the setup is under The data used on this page). Answering the question means discovering that the two files are in different coordinate systems, checking whether any points are unusable, working out that some fall outside every ward, deciding whether to map counts or rates, and producing something a colleague can read. None of that is exotic — it is the ordinary shape of a first real analysis.

This guide is one script, start to finish, with the decisions made explicitly rather than skipped.

Quick answer

Seven stages. Most tutorials cover the fourth and fifth.

import geopandas as gpd
import matplotlib.pyplot as plt

CRS = 27700                                          # a metric CRS for the area

wards = gpd.read_file("wards.gpkg").to_crs(CRS)
incidents = gpd.read_file("incidents.geojson").to_crs(CRS)

joined = gpd.sjoin(incidents, wards[["ward_code", "geometry"]],
                   how="left", predicate="within")

counts = joined.groupby("ward_code").size().rename("incidents")
wards = wards.join(counts, on="ward_code").fillna({"incidents": 0})
wards["per_1000"] = wards["incidents"] / wards["population"] * 1000

ax = wards.plot(column="per_1000", scheme="quantiles", k=5, cmap="YlOrRd",
                legend=True, figsize=(10, 11), edgecolor="white", linewidth=0.3)
ax.set_axis_off()
plt.savefig("incidents.png", dpi=200, bbox_inches="tight", facecolor="white")

Stage	What it decides
inspect	whether the data is what you think
clean	which rows survive, and why
align	one CRS, chosen for the area
join	which incident belongs to which ward
aggregate	count, or rate
map	classification, colour, missing data
export	what a reader receives

What you get — from 1. Look at what you have before writing any analysis

wards.gpkg
  33 features, Polygon, EPSG:27700
  fields: ['ward_code', 'ward_name', 'population', 'area_ha']
incidents.geojson
  15,942 features, Point, EPSG:4326
  fields: ['id', 'category', 'reported', 'lat', 'lon']

Look, clean, one CRS, join, map a rate

1. Look at what you have before writing any analysis

import geopandas as gpd
import pyogrio

for path in ["wards.gpkg", "incidents.geojson"]:
    info = pyogrio.read_info(path)
    print(f"{path}")
    print(f"  {info['features']:,} features, {info['geometry_type']}, {info['crs']}")
    print(f"  fields: {list(info['fields'])[:8]}")

wards.gpkg
  33 features, Polygon, EPSG:27700
  fields: ['ward_code', 'ward_name', 'population', 'area_ha']
incidents.geojson
  15,942 features, Point, EPSG:4326
  fields: ['id', 'category', 'reported', 'lat', 'lon']

Two things are already clear: the CRS differ, and one file is GeoJSON so it is necessarily WGS 84. Reading headers first costs milliseconds and shapes everything after.

Now read and look properly:

wards = gpd.read_file("wards.gpkg")
incidents = gpd.read_file("incidents.geojson")

print(incidents.head(3).to_string())
print(f"\ngeometry types: {incidents.geom_type.value_counts().to_dict()}")
print(f"null geometry:  {incidents.geometry.isna().sum():,}")
print(f"empty geometry: {incidents.geometry.is_empty.sum():,}")
print(f"nulls per column:\n{incidents.isna().sum().to_string()}")

          id       category reported        lat       lon                   geometry
0  135761323  violent-crime  2026-06  53.709572 -1.640712  POINT (-1.64071 53.70957)
1  135768947    other-crime  2026-06  53.780117 -1.713149  POINT (-1.71315 53.78012)
2  135771893  violent-crime  2026-06  53.839863 -1.778529  POINT (-1.77853 53.83986)

geometry types: {'Point': 15942}
null geometry:  0
empty geometry: 0
nulls per column:
id          0
category    0
reported    0
lat         0
lon         0
geometry    0

No empty or null geometry here: police.uk only publishes crimes it could place, and the ones it could not are simply absent from this feed — which is itself worth knowing, because the gap is invisible in the file. In a register you export yourself, a few percent of records with no coordinates is normal, and how you handle them changes the answer. The checks cost nothing when they find nothing. See null, empty, missing and invalid.

2. Clean, and record what you removed

import pandas as pd


def clean_incidents(gdf):
    """Remove what cannot be mapped, and report what went."""
    report = {"input": len(gdf)}

    gdf = gdf[gdf.geometry.notna() & ~gdf.geometry.is_empty].copy()
    report["dropped_no_geometry"] = report["input"] - len(gdf)

    # a point at exactly (0, 0) is a failed geocode, not the Gulf of Guinea
    at_null_island = (gdf.geometry.x.abs() < 1e-6) & (gdf.geometry.y.abs() < 1e-6)
    gdf = gdf[~at_null_island]
    report["dropped_null_island"] = int(at_null_island.sum())

    before = len(gdf)
    gdf = gdf.drop_duplicates(subset=["id"])
    report["dropped_duplicate_ids"] = before - len(gdf)

    gdf["category"] = gdf["category"].fillna("unknown").str.strip().str.lower()
    gdf["reported"] = pd.to_datetime(gdf["reported"], errors="coerce")
    report["unparseable_dates"] = int(gdf["reported"].isna().sum())

    report["output"] = len(gdf)
    for key, value in report.items():
        print(f"  {key:<24} {value:>8,}")
    return gdf, report

incidents, cleaning = clean_incidents(incidents)

  input                      15,942
  dropped_no_geometry             0
  dropped_null_island             0
  dropped_duplicate_ids           0
  unparseable_dates               0
  output                     15,942

Separate categories, separate counts — here all zeros, which is a result too: this feed is clean. Collapsing them into "dropped n rows" would hide why rows went when they do go: records geocoded to Null Island are a specific and fixable upstream problem, not random missingness, and a duplicated id usually means a file was appended to twice.

Keep the report. It goes on the final map, so a reader knows what the analysis covers.

3. Put everything in one metric CRS

CRS = 27700          # British National Grid — metres, designed for this area

print(f"wards      {wards.crs}  →  EPSG:{CRS}")
print(f"incidents  {incidents.crs}  →  EPSG:{CRS}")

wards = wards.to_crs(CRS)
incidents = incidents.to_crs(CRS)

assert wards.crs == incidents.crs

wards      EPSG:27700  →  EPSG:27700
incidents  EPSG:4326  →  EPSG:27700

The wards already arrive in British National Grid; the police points arrive in longitude and latitude. Two reasons this matters. A spatial join between layers in different CRS returns nothing — no error, just an empty result. And measurements in EPSG:4326 are in degrees, so any area or density calculation is meaningless.

Reproject once, here, not repeatedly later. The reasoning is in projected vs geographic CRS, and choosing the code is covered in how to choose the right projected CRS.

Check the extents overlap before joining:

from shapely.geometry import box

if not box(*wards.total_bounds).intersects(box(*incidents.total_bounds)):
    raise SystemExit(f"the layers do not overlap:\n  wards {wards.total_bounds}\n"
                     f"  incidents {incidents.total_bounds}")

The numbers steps 2 and 4 print. Knowing which step removed the rows is what makes the coverage figure honest.

4. Join, and check the cardinality

joined = gpd.sjoin(
    incidents,
    wards[["ward_code", "ward_name", "geometry"]],
    how="left",
    predicate="within",
)

print(f"incidents in:  {len(incidents):,}")
print(f"rows out:      {len(joined):,}")
print(f"unmatched:     {joined['ward_code'].isna().sum():,}")
print(f"duplicated:    {joined.index.duplicated().sum():,}")

incidents in:  15,942
rows out:      15,942
unmatched:     6,718
duplicated:    0

Two things to understand rather than ignore.

Unmatched incidents fall inside no ward — here, because the crimes were fetched for the rectangle around Leeds, and the rectangle takes in parts of Bradford (4,419 of them), Kirklees (1,134), Wakefield (950) and Calderdale (197). how="left" kept them so they are visible; an inner join would have dropped them silently.

Duplicated rows would be points on a shared ward boundary, within both neighbours — none here, because police.uk snaps crimes to points on streets, not to boundaries. The guard costs one line:

joined = joined[~joined.index.duplicated(keep="first")]

The row-count check is two lines and catches the commonest silent error in spatial work — see spatial join cardinality explained.

Look at where the unmatched points are before dismissing them:

unmatched = joined[joined["ward_code"].isna()]
if len(unmatched):
    print(f"unmatched extent: {unmatched.total_bounds.round().tolist()}")
    print(f"study extent:     {wards.total_bounds.round().tolist()}")

unmatched extent: [413247.0, 422540.0, 446906.0, 449853.0]
study extent:     [413216.0, 422593.0, 446879.0, 450175.0]

If they cluster in one place, the ward layer has a gap. If they are scattered around the edge, they are genuinely outside the study area.

5. Aggregate — and map a rate, not a count

counts = joined.groupby("ward_code").size().rename("incidents")

wards = wards.merge(counts, left_on="ward_code", right_index=True, how="left")
wards["incidents"] = wards["incidents"].fillna(0).astype(int)

wards["area_km2"] = wards.geometry.area / 1e6
wards["per_1000_people"] = wards["incidents"] / wards["population"].replace(0, pd.NA) * 1000
wards["per_km2"] = wards["incidents"] / wards["area_km2"]

print(wards[["ward_name", "incidents", "population", "per_1000_people"]]
      .sort_values("per_1000_people", ascending=False).head(5).to_string(index=False))

                ward_name  incidents  population  per_1000_people
Little London & Woodhouse       1124       31076        36.169391
      Hunslet & Riverside        589       26474        22.248244
                   Armley        504       26725        18.858746
       Gipton & Harehills        572       30879        18.523916
        Beeston & Holbeck        418       26593        15.718422

how="left" on the merge is essential. A ward with zero incidents does not appear in counts, and an inner merge would drop it from the map entirely — leaving a hole a reader interprets as missing data rather than as zero.

And the rate is the point. A count map would largely show where people are, because a choropleth colours areas and the reader integrates colour across them. Leeds' wards are drawn to have similar populations, so here count and rate agree more than they usually do — Example 2 measures exactly how much. The full argument is in choropleth classification explained.

Sanity-check the totals:

assert wards["incidents"].sum() + len(unmatched) == len(joined), "incidents lost"
print(f"{wards['incidents'].sum():,} mapped + {len(unmatched):,} unmatched "
      f"= {len(joined):,}")

9,224 mapped + 6,718 unmatched = 15,942

6. Map it, and say what the map is

import matplotlib.pyplot as plt
import matplotlib.patheffects as pe

fig, ax = plt.subplots(figsize=(10, 11))

wards.plot(
    column="per_1000_people", scheme="quantiles", k=5, cmap="YlOrRd", ax=ax,
    edgecolor="white", linewidth=0.35, legend=True,
    legend_kwds={"loc": "lower right", "title": "Incidents per 1,000 residents",
                 "fmt": "{:,.1f}", "fontsize": 9, "framealpha": 0.9},
    missing_kwds={"color": "#e2e8f0", "hatch": "///", "label": "no population data"},
)

for row in wards.nlargest(6, "per_1000_people").itertuples():
    p = row.geometry.representative_point()
    ax.annotate(row.ward_name, (p.x, p.y), ha="center", fontsize=8.5,
                path_effects=[pe.withStroke(linewidth=2.5, foreground="white")])

ax.set_title("Reported incidents per 1,000 residents", fontsize=15, loc="left")
ax.annotate("By electoral ward, Leeds, June 2026", xy=(0, 1.005), xycoords="axes fraction",
            fontsize=10, color="#475569", va="bottom")
ax.set_axis_off()

ax.annotate(
    f"{wards['incidents'].sum():,} incidents mapped of {cleaning['input']:,} reported "
    f"({100 * wards['incidents'].sum() / cleaning['input']:.0f}%). "
    f"{cleaning['dropped_no_geometry']:,} had no coordinates; "
    f"{len(unmatched):,} fell outside every ward.\n"
    f"Quantile classification, 5 classes. EPSG:27700. Crimes: data.police.uk (OGL v3). "
    f"Boundaries and population: ONS, Nomis © Crown copyright.",
    xy=(0, -0.02), xycoords="axes fraction", fontsize=7.5, color="#64748b", va="top")

fig.savefig("incidents_per_1000.png", dpi=200, bbox_inches="tight", facecolor="white")

The caption is doing real work. It states the coverage, why the rest is missing, the classification scheme and the CRS. A reader can then judge the map rather than trust it — and it takes four lines.

Two numbers change at the join. Both are invisible on the finished map.

Leeds crime per ward, as one runnable script

The script below is these four steps and a main() that calls them in order.

The data used on this page

wards.gpkg is Leeds' 33 electoral wards (December 2024) from the ONS, in British National Grid, with each ward's Census 2021 population from Nomis. incidents.geojson is every street-level crime West Yorkshire Police published for June 2026 inside the rectangle around Leeds, from data.police.uk, in WGS 84 — GeoJSON's only CRS. Both are also written to data/, where the script in Example 1 looks.

import time
from pathlib import Path

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
wards = gpd.read_file(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
                      "?where=LAD24NM%3D%27Leeds%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson")
wards = (wards.set_crs(27700, allow_override=True)
              .rename(columns={"WD24CD": "ward_code", "WD24NM": "ward_name"}))
population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(wards['ward_code'])}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
wards["population"] = wards["ward_code"].map(population)
wards["area_ha"] = (wards.geometry.area / 1e4).round(1)
wards = wards[["ward_code", "ward_name", "population", "area_ha", "geometry"]]

west, south, east, north = wards.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):   # ≤10,000 a request
        for attempt in range(5):                        # the API answers bursts with 429
            r = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                              data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if r.status_code != 429:
                break
            time.sleep(2 ** attempt)
        r.raise_for_status()
        crimes += [{"id": c["id"], "category": c["category"], "reported": c["month"],
                    "lat": float(c["location"]["latitude"]), "lon": float(c["location"]["longitude"])}
                   for c in r.json()]
crimes = pd.DataFrame(crimes).sort_values("id", ignore_index=True)   # the API returns them in a varying order
incidents = gpd.GeoDataFrame(crimes, geometry=gpd.points_from_xy(crimes["lon"], crimes["lat"]), crs=4326)

Path("data").mkdir(exist_ok=True)
for folder in (Path("."), Path("data")):
    wards.to_file(folder / "wards.gpkg")
    incidents.to_file(folder / "incidents.geojson")
print(f"{len(wards)} wards, {len(incidents):,} crimes")

33 wards, 15,942 crimes

Example 1: the whole analysis as one runnable script

#!/usr/bin/env python3
"""Incidents per 1,000 residents by ward. Run: python analysis.py"""
from pathlib import Path
import json
import geopandas as gpd
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe

CRS = 27700
DATA = Path("data")
OUT = Path("output")
OUT.mkdir(exist_ok=True)


def load(path, crs):
    gdf = gpd.read_file(path)
    if gdf.crs is None:
        raise ValueError(f"{path} has no CRS — identify it before using it")
    return gdf.to_crs(crs)


def clean_points(gdf):
    report = {"input": len(gdf)}
    gdf = gdf[gdf.geometry.notna() & ~gdf.geometry.is_empty].copy()
    report["no_geometry"] = report["input"] - len(gdf)
    origin = (gdf.geometry.x.abs() < 1e-6) & (gdf.geometry.y.abs() < 1e-6)
    report["null_island"] = int(origin.sum())
    gdf = gdf[~origin]
    before = len(gdf)
    gdf = gdf.drop_duplicates(subset=["id"])
    report["duplicate_ids"] = before - len(gdf)
    report["output"] = len(gdf)
    return gdf, report


def assign_to_areas(points, areas, key="ward_code"):
    joined = gpd.sjoin(points, areas[[key, "geometry"]], how="left",
                       predicate="within")
    duplicated = int(joined.index.duplicated().sum())
    joined = joined[~joined.index.duplicated(keep="first")]
    unmatched = int(joined[key].isna().sum())
    if len(joined) != len(points):
        raise AssertionError(f"{len(points)} in, {len(joined)} out")
    return joined, {"boundary_duplicates": duplicated, "unmatched": unmatched}


def summarise(areas, joined, key="ward_code", per=1_000):
    counts = joined.groupby(key).size().rename("incidents")
    out = areas.merge(counts, left_on=key, right_index=True, how="left")
    out["incidents"] = out["incidents"].fillna(0).astype(int)
    out["area_km2"] = out.geometry.area / 1e6
    out["per_1000_people"] = (out["incidents"]
                              / out["population"].replace(0, pd.NA) * per)
    out["per_km2"] = out["incidents"] / out["area_km2"]
    return out


def draw(areas, notes, path):
    fig, ax = plt.subplots(figsize=(10, 11))
    areas.plot(column="per_1000_people", scheme="quantiles", k=5, cmap="YlOrRd",
               ax=ax, edgecolor="white", linewidth=0.35, legend=True,
               legend_kwds={"loc": "lower right", "fmt": "{:,.1f}",
                            "title": "Incidents per 1,000 residents", "fontsize": 9},
               missing_kwds={"color": "#e2e8f0", "hatch": "///",
                             "label": "no population data"})
    for row in areas.nlargest(6, "per_1000_people").itertuples():
        p = row.geometry.representative_point()
        ax.annotate(row.ward_name, (p.x, p.y), ha="center", fontsize=8.5,
                    path_effects=[pe.withStroke(linewidth=2.5, foreground="white")])
    ax.set_title("Reported incidents per 1,000 residents", fontsize=15, loc="left")
    ax.set_axis_off()
    ax.annotate(notes, xy=(0, -0.02), xycoords="axes fraction",
                fontsize=7.5, color="#64748b", va="top")
    fig.savefig(path, dpi=200, bbox_inches="tight", facecolor="white")
    plt.close(fig)
    return path


def main():
    wards = load(DATA / "wards.gpkg", CRS)
    incidents = load(DATA / "incidents.geojson", CRS)

    incidents, cleaning = clean_points(incidents)
    joined, join_report = assign_to_areas(incidents, wards)
    wards = summarise(wards, joined)

    mapped = int(wards["incidents"].sum())
    assert mapped + join_report["unmatched"] == len(joined), "incidents lost"

    notes = (
        f"{mapped:,} of {cleaning['input']:,} reported incidents mapped "
        f"({100 * mapped / cleaning['input']:.0f}%). "
        f"{cleaning['no_geometry']:,} had no coordinates, "
        f"{cleaning['null_island']:,} were geocoded to (0, 0), "
        f"{join_report['unmatched']:,} fell outside every ward.\n"
        f"Quantiles, k=5. EPSG:{CRS}. Crimes: data.police.uk; "
        f"boundaries and population © Crown copyright.")

    draw(wards, notes, OUT / "incidents_per_1000.png")
    wards.to_file(OUT / "wards_with_rates.gpkg", driver="GPKG")
    wards.drop(columns="geometry").to_csv(OUT / "wards_with_rates.csv", index=False)
    (OUT / "run_report.json").write_text(
        json.dumps({"cleaning": cleaning, "join": join_report,
                    "mapped": mapped, "crs": CRS}, indent=2))

    print(f"\n{mapped:,} incidents across {len(wards)} wards")
    print(wards.nlargest(5, "per_1000_people")[
        ["ward_name", "incidents", "population", "per_1000_people"]
    ].to_string(index=False))
    print(f"\nwrote {OUT}/")


main()          # in analysis.py this line sits under `if __name__ == "__main__":`


9,224 incidents across 33 wards
                ward_name  incidents  population  per_1000_people
Little London & Woodhouse       1124       31076        36.169391
      Hunslet & Riverside        589       26474        22.248244
                   Armley        504       26725        18.858746
       Gipton & Harehills        572       30879        18.523916
        Beeston & Holbeck        418       26593        15.718422

wrote output/

Four properties make this a script rather than a notebook dump. Each stage is a function with a return value, so each can be tested without files. The assertion after the aggregation catches a whole class of silent error — incidents disappearing between the join and the summary. The reports are written to JSON, so the numbers on the map can be traced later. And it produces three artefacts: a map for reading, a GeoPackage for further analysis, and a CSV for a spreadsheet.

matplotlib.use("Agg") before importing pyplot means it runs on a server or in cron with no display.

Example 2: checking the answer three ways

An analysis that has not been checked is a guess with formatting:

import numpy as np

def sanity_checks(wards, joined, incidents, cleaning, unmatched_n):
    checks = []

    total = wards["incidents"].sum()
    checks.append(("all incidents accounted for",
                   total + unmatched_n == len(joined),
                   f"{total:,} + {unmatched_n:,} vs {len(joined):,}"))

    checks.append(("no negative rates",
                   (wards["per_1000_people"].dropna() >= 0).all(),
                   f"min {wards['per_1000_people'].min():.2f}"))

    # a rate above ~200 per 1,000 is possible but worth a look
    extreme = wards[wards["per_1000_people"] > 200]
    checks.append(("no implausible rates", len(extreme) == 0,
                   f"{len(extreme)} wards above 200/1,000"))

    # the count map and the rate map should disagree — if not, check the denominator
    corr = wards[["incidents", "per_1000_people"]].corr().iloc[0, 1]
    checks.append(("rate differs from raw count", corr < 0.95,
                   f"correlation {corr:.2f}"))

    # coverage
    coverage = 100 * total / cleaning["input"]
    checks.append(("coverage above 90%", coverage > 90, f"{coverage:.1f}%"))

    # every ward present
    checks.append(("no wards dropped", wards["ward_code"].notna().all(),
                   f"{len(wards)} wards"))

    for name, passed, detail in checks:
        print(f"  {'✓' if passed else '✗'} {name:<32} {detail}")
    failures = [c[0] for c in checks if not c[1]]
    if failures:
        raise AssertionError(f"failed: {failures}")
    return checks

try:
    checks = sanity_checks(wards, joined, incidents, cleaning, len(unmatched))
except AssertionError as failure:
    print(f"\nAssertionError: {failure}")

  ✓ all incidents accounted for      9,224 + 6,718 vs 15,942
  ✓ no negative rates                min 3.89
  ✓ no implausible rates             0 wards above 200/1,000
  ✗ rate differs from raw count      correlation 0.99
  ✗ coverage above 90%               57.9%
  ✓ no wards dropped                 33 wards

AssertionError: failed: ['rate differs from raw count', 'coverage above 90%']

Two checks failed, and both failures are right.

Coverage is 57.9% because the crimes were fetched for the rectangle around Leeds, and four in ten of them are in Bradford, Kirklees, Wakefield or Calderdale. The check cannot know that; it only knows that a map claiming to show the incidents covers barely half of what was loaded. The fix is to fetch or clip to the study area first — or to change the denominator of the check to "incidents inside the study area" and say so on the map.

The correlation check is the interesting one. At 0.99, the rate map is almost the same picture as the raw count, and the denominator is doing nothing. That is not a bug in the data: Leeds' wards are drawn to have similar populations — 18,500 to 31,200 residents — so dividing by population barely reorders them. Where areas differ in size, as LSOAs or countries do, the same check passes easily. Here it is telling you that "per 1,000 residents" adds little, and that crimes per km², or a finer geography, would be the more informative map.

Raising on failure means a broken run stops rather than producing a confident wrong map.

Example 3: making it repeatable next month

import argparse
from pathlib import Path

def parse_args():
    p = argparse.ArgumentParser(description="Incident rates by ward")
    p.add_argument("--incidents", type=Path, default=Path("data/incidents.geojson"))
    p.add_argument("--areas", type=Path, default=Path("data/wards.gpkg"))
    p.add_argument("--out", type=Path, default=Path("output"))
    p.add_argument("--crs", type=int, default=27700)
    p.add_argument("--area-key", default="ward_code")
    p.add_argument("--population-col", default="population")
    p.add_argument("--scheme", default="quantiles",
                   choices=["quantiles", "natural_breaks", "equal_interval"])
    p.add_argument("--classes", type=int, default=5)
    p.add_argument("--since", help="only incidents on or after this date, YYYY-MM-DD")
    return p.parse_args()

def main():
    args = parse_args()
    args.out.mkdir(parents=True, exist_ok=True)

    areas = load(args.areas, args.crs)
    points = load(args.incidents, args.crs)

    if args.since:
        before = len(points)
        points["reported"] = pd.to_datetime(points["reported"], errors="coerce")
        points = points[points["reported"] >= pd.Timestamp(args.since)]
        print(f"date filter {args.since}: {before:,} → {len(points):,}")

    points, cleaning = clean_points(points)
    joined, join_report = assign_to_areas(points, areas, key=args.area_key)
    areas = summarise(areas, joined, key=args.area_key)
    ...

python analysis.py --since 2026-01-01 --scheme natural_breaks --out output/2026
python analysis.py --areas data/lsoa.gpkg --area-key lsoa_code --out output/lsoa

Turning constants into arguments is what separates a one-off from something reusable. The second command runs the same analysis at a different geography with no code change — and the fact that it can is a decent test that the code contains no hidden assumptions about wards.

Every argument here corresponds to a decision made earlier in this article: which CRS, which classification, which denominator, which time window. Making them explicit means a reader of the command can see what was chosen. See how to turn a GIS script into a command-line tool.

The decisions between the steps

Every stage forces a decision. Tutorials cover the operations; the decisions are the analysis.

The gap between "I can do a spatial join" and "I can answer a question with data" is not technique. It is the set of decisions that sit between the steps, and the tutorials skip them because each one is specific to a question.

Inspection comes first because everything downstream depends on it. Reading headers before reading data told us the CRS differed, that one file was necessarily WGS 84, and that a population column existed. Each of those shaped a later decision. The alternative — reading everything and discovering problems as they break things — costs more time and produces surprises in the middle rather than at the start.

Cleaning is where the analysis acquires assumptions, and they should be visible. Dropping records with no coordinates is defensible; not saying so is not. If they are concentrated in one area — a district whose addresses geocode badly — then the map systematically under-represents it, and only the caption gives a reader any chance of knowing. police.uk's feed had none to drop, but that is because crimes it cannot place never reach the feed at all, which the caption should also admit. Reporting each category separately makes the fixable ones visible: points at Null Island are a specific upstream bug, not random missingness.

The CRS decision is load-bearing in two directions. A join between mismatched CRS returns nothing, with no error. And measurements in degrees are meaningless, so a density computed in EPSG:4326 is a number with no units. Reprojecting once, early, to a CRS chosen for the area resolves both — and doing it once rather than repeatedly matters, since to_crs is among the most expensive operations available.

The join is where row counts change, and checking them is the cheapest insurance in spatial work. More rows out than in means duplication — points on shared boundaries, when there are any. Fewer means loss. Both are usually accidents, both are invisible on a map, and both take two lines to detect. The how="left" choice matters equally: it kept the 6,718 unmatched points visible instead of silently dropping them, turning "the download covered more than Leeds" from an unknown into a number.

Which is the next paragraph: the caption is part of the analysis, not decoration.

The count-versus-rate decision is the one that changes the conclusion. A choropleth encodes value as colour over area, and a reader's eye integrates colour across area whether or not that was intended. So mapping raw counts produces a picture of where the denominator is large — usually population. Two of the top five wards by rate here are not in the top five by count, which means the count map and the rate map support different recommendations from identical data.

Finally, the caption is part of the analysis, not decoration. Coverage, exclusions, classification scheme and CRS are four choices that determine what the map appears to say. A reader who can see them can evaluate the argument; one who cannot is being asked to trust four decisions they have no way to inspect. It costs four lines, and it is what separates a chart from a finding.

Empty geometry, (0, 0) points and boundary points

Read headers before data. pyogrio.read_info costs milliseconds and prevents a wrong read.
Empty geometry is not null geometry. Test for both — notna() and ~is_empty.
Points at exactly (0, 0) are failed geocodes, not data. Filter them explicitly.
how="left" on both the join and the merge, or unmatched incidents and zero-incident wards vanish.
Compare row counts across every join. More means duplication, fewer means loss.
A point on a shared boundary is within both polygons. Deduplicate deliberately.
fillna(0) on a count is correct; on a rate or a mean it is not. Zero incidents is a fact; unknown population is not.
Reproject once, at the start. to_crs builds a new geometry per feature.
matplotlib.use("Agg") before importing pyplot for anything that runs unattended.
Write the report as JSON, so the numbers in the caption can be traced back later.

Internal links

How to explore a spatial dataset you have never seen before — the inspection step in depth
The Python GIS data cleaning checklist — the cleaning step in depth
How to perform a spatial join in Python (GeoPandas) — the join step in depth
Spatial join cardinality explained — why the row count changed
How to create a choropleth map in Python with GeoPandas — the mapping step in depth
Choropleth classification explained — count versus rate, and the scheme
How to choose the right projected CRS for your study area — the CRS decision
How to turn a GIS script into a command-line tool with argparse — making it repeatable

FAQ

Where do I start with a dataset I have never seen?

Read the headers first — feature count, geometry type, CRS, field names. Then check for null and empty geometry, duplicates and missing values before writing any analysis.

Why does my spatial join return nothing?

Almost always because the two layers are in different coordinate systems. Print both crs and total_bounds; degrees next to metres is the giveaway.

Should I map counts or rates?

Rates, essentially always. A choropleth colours areas, so a count map shows where the denominator is large — usually population — regardless of what the legend says.

What do I do with rows that have no coordinates?

Remove them and report how many. If they cluster geographically, the map under-represents that area, and only the caption gives a reader a chance to know.

My join returned more rows than I put in. Why?

Points on a shared boundary are within both neighbouring polygons. Deduplicate on the index deliberately rather than letting the first match win by accident.

Why do wards with zero incidents disappear?

An inner merge drops them, because they are absent from the count. Use how="left" and fillna(0) — a hole in a choropleth reads as missing data, not as zero.

What should the map caption say?

Coverage, what was excluded and why, the classification scheme and the CRS. Those four choices determine what the map appears to say, and a reader cannot judge it without them.