⚙️ How-To Getting started Maps & visualisation Updated

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

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?

Answering it means downloading two datasets, discovering that one is in a different coordinate system, finding that 4% of the points have no coordinates at all, working out that some points 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

Flow from acquire through inspect, clean, join, aggregate and map to export.
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

Step-by-step solution

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
  215 features, MultiPolygon, EPSG:27700
  fields: ['ward_code', 'ward_name', 'population', 'area_ha']
incidents.geojson
  8,436 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))
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()}")
geometry types: {'Point': 8436}
null geometry:  0
empty geometry: 341
nulls per column:
id            0
category      12
reported      0
lat         341
lon         341
geometry      0

341 records with empty geometry and no coordinates — 4% of the data. That is a finding, not a nuisance, and how you handle it changes the answer. See null, empty, missing and invalid.

2. Clean, and record what you removed 

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"] = gpd.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                       8,436
  dropped_no_geometry           341
  dropped_null_island            18
  dropped_duplicate_ids           4
  unparseable_dates               0
  output                      8,073

Three separate categories, three separate counts. Collapsing them into "dropped 363 rows" would hide that 18 records were geocoded to Null Island — a specific and fixable upstream problem, not random missingness.

Keep the report. It goes on the final map, so a reader knows the analysis covers 96% of reports rather than all of them.

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

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

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:  8,073
rows out:      8,075
unmatched:     112
duplicated:    2

Two findings, both worth understanding rather than ignoring.

112 unmatched incidents fall inside no ward — outside the study area, or in a gap in the boundaries. how="left" kept them so they are visible; an inner join would have dropped them silently.

2 duplicated rows are points on a shared ward boundary, which are within both neighbours. Resolve deliberately rather than letting the first one win:

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: {[round(v) for v in unmatched.total_bounds]}")
    print(f"study extent:     {[round(v) for v in wards.total_bounds]}")

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
   City Centre         412       11,204            36.77
      Ardwick          188        9,882            19.02
     Ancoats           174       12,441            13.99
   Piccadilly          166       14,208            11.68
    Longsight          142       17,004             8.35

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 show where people are, because a choropleth colours areas and the reader integrates colour across them. Two of the top five wards by rate are not in the top five by count. 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):,}")

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, 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. "
    f"Source: incident register, ward boundaries © 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 (96%), 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.

Scene showing incidents joined to wards with unmatched points and boundary duplicates called out.
Two numbers change at the join. Both are invisible on the finished map.

Code examples

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}. Source: incident register; "
        f"boundaries © 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}/")


if __name__ == "__main__":
    main()
8,073 incidents across 215 wards
  ward_name  incidents  population  per_1000_people
City Centre        412      11,204            36.77
    Ardwick        188       9,882            19.02
    Ancoats        174      12,441            13.99
 Piccadilly        166      14,208            11.68
  Longsight        142      17,004             8.35

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

sanity_checks(wards, joined, incidents, cleaning, join_report["unmatched"])
  ✓ all incidents accounted for      8,073 + 112 vs 8,185
  ✓ no negative rates                min 0.00
  ✓ no implausible rates             0 wards above 200/1,000
  ✓ rate differs from raw count      correlation 0.61
  ✓ coverage above 90%               96.0%
  ✓ no wards dropped                 215 wards

The correlation check is the interesting one. If the rate map correlated above 0.95 with the raw count, the denominator would be doing nothing — usually because population is roughly uniform across the areas, or because the wrong column was used. At 0.61 the rate genuinely tells a different story from the count, which is the whole reason for computing it.

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.

Explanation

Grid mapping each pipeline stage to the decision it forces and what goes wrong if skipped.
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 341 records with no coordinates is defensible; not saying so is not. If those 341 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. Reporting each category separately also makes the fixable ones visible: 18 points at Null Island is 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 — here, points on shared boundaries. 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 112 unmatched points visible instead of silently dropping them, turning "some incidents are outside the study area" from an unknown into a number.

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.

Edge cases or notes

  • 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.

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.