How to Batch-Generate a Map Image for Every Region in Python

Problem statement

Someone needs one PNG per ward. Greater Manchester has 215 of them (the setup is under The data used on this page). The loop is obvious:

for name, part in gdf.groupby("ward_name"):
    fig, ax = plt.subplots(figsize=(8, 8))
    part.plot(ax=ax)
    fig.savefig(f"out/{name}.png")

Run it and four things go wrong.

Memory climbs with every ward, because every figure is still open. Bury has a ward called St. Mary's and Oldham one called St Mary's; Manchester and Trafford each have a Brooklands; so some files overwrite others. Every map is at a different scale, so a small ward looks the same size as a large one and nobody can compare them. And the map shows the ward alone on a white background, with no surrounding context to say where it is.

Batch cartography is a batch job that happens to draw pictures, and it needs the same things any batch job needs — bounded resources, safe filenames, per-item error handling, a report — plus one cartographic requirement the others do not have: consistency between the outputs.

Quick answer

from pathlib import Path
import re
import geopandas as gpd
import matplotlib
matplotlib.use("Agg")                      # no display needed, and much faster
import matplotlib.pyplot as plt

def safe(name):
    return re.sub(r"[^\w\-]+", "_", str(name)).strip("_")[:80]

gdf = gpd.read_file("wards.gpkg").to_crs(27700)
out = Path("out"); out.mkdir(exist_ok=True)

vmin, vmax = gdf["density"].min(), gdf["density"].max()    # one shared scale

for row in gdf.itertuples():
    fig, ax = plt.subplots(figsize=(7, 7))
    gdf.plot(ax=ax, facecolor="#f1f5f9", edgecolor="#e2e8f0", linewidth=0.3)
    gpd.GeoSeries([row.geometry], crs=gdf.crs).plot(
        ax=ax, facecolor="#ef4444", alpha=0.6, edgecolor="#991b1b", linewidth=1.5)
    minx, miny, maxx, maxy = row.geometry.bounds
    pad = 0.35 * max(maxx - minx, maxy - miny)
    ax.set_xlim(minx - pad, maxx + pad); ax.set_ylim(miny - pad, maxy + pad)
    ax.set_title(row.ward_name, fontsize=12)
    ax.set_axis_off()
    fig.savefig(out / f"{safe(row.ward_name)}.png", dpi=100,
                bbox_inches="tight", facecolor="white")
    plt.close(fig)                          # ← the line that prevents the crash

print(f"{len(gdf)} wards rendered, {len(list(out.glob('*.png')))} files on disk")

215 wards rendered, 212 files on disk

Three wards' maps were silently replaced by another ward's — step 2 finds which.

Agg, safe filenames, comparable maps, per-region failures

1. Use the Agg backend and close every figure

import matplotlib
matplotlib.use("Agg")           # before importing pyplot
import matplotlib.pyplot as plt

Agg is a file-only renderer. It needs no display, works over SSH and in Docker, and is faster than an interactive backend. Setting it must happen before pyplot is imported.

plt.close(fig) is the line that decides whether the job finishes. Matplotlib keeps a reference to every open figure so that plt.show() can find them, so a loop that never closes accumulates all of them:

import tracemalloc
import warnings

Path("leak").mkdir(exist_ok=True)
tracemalloc.start()
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    for i, row in enumerate(gdf.head(41).itertuples()):
        fig, ax = plt.subplots(figsize=(7, 7))
        gdf.plot(ax=ax)
        fig.savefig(f"leak/{i}.png", dpi=100)
        # plt.close(fig)
        if i % 10 == 0:
            print(f"  {i:>3} figures open: {len(plt.get_fignums()):>3}  "
                  f"memory {tracemalloc.get_traced_memory()[0] / 1e6:>7.1f} MB")
tracemalloc.stop()
print("warnings:", sorted({str(w.message)[:60] for w in caught}))
plt.close("all")

    0 figures open:   1  memory     1.1 MB
   10 figures open:  11  memory    11.6 MB
   20 figures open:  21  memory    22.2 MB
   30 figures open:  31  memory    31.4 MB
   40 figures open:  41  memory    41.4 MB
warnings: ['More than 20 figures have been opened. Figures created throu']

Memory grows with every figure and never comes back. With plt.close(fig) the number stays at one. Matplotlib does warn once more than 20 figures are open — but it warns once per figure after that, in a log that is already busy, and a job writing to a file rather than a terminal never shows it to anyone.

2. Make filenames safe and unique

import re
from collections import Counter

def safe(name, max_len=80):
    s = re.sub(r"[^\w\-]+", "_", str(name)).strip("_")
    return (s or "unnamed")[:max_len]

names = [safe(n) for n in gdf["ward_name"]]
clashes = {n: c for n, c in Counter(names).items() if c > 1}
if clashes:
    print(f"filename collisions: {clashes}")

filename collisions: {'Brooklands': 2, 'St_Mary_s': 2, 'Manor': 2}

St. Mary's and St Mary's differ only in the full stop that safe() removes; the two Brooklands and the two Manors are the same name in different boroughs. Disambiguate with a stable identifier rather than a counter, so re-running produces the same names:

gdf["_file"] = gdf.apply(lambda r: f"{safe(r['ward_name'])}_{r['ward_code']}", axis=1)
print(f"{gdf['_file'].nunique()} unique filenames for {len(gdf)} wards")

215 unique filenames for 215 wards

A stable filename is what makes the job safe to re-run — a counter-based name changes if the input order changes, so yesterday's outputs and today's do not correspond.

3. Make the maps comparable

Matplotlib scales each figure independently, which produces a set of maps in which every region fills the frame and nothing can be compared. Three things need fixing:

Shared colour range — otherwise the same colour means a different value in each map:

vmin, vmax = gdf["density"].quantile([0.02, 0.98])
part.plot(column="density", cmap="YlOrRd", vmin=vmin, vmax=vmax, ax=ax)

Clipping at the 2nd and 98th percentiles keeps one extreme region from flattening the range for everyone else.

Consistent scale — one of two strategies, and they answer different questions:

# (a) same extent everywhere — comparable size, small regions become tiny
minx, miny, maxx, maxy = gdf.total_bounds
ax.set_xlim(minx, maxx); ax.set_ylim(miny, maxy)

# (b) same map scale, centred per region — comparable size, each region fills its map
CENTRE_SPAN = 12_000                     # metres across every map
cx, cy = row.geometry.centroid.coords[0]
ax.set_xlim(cx - CENTRE_SPAN/2, cx + CENTRE_SPAN/2)
ax.set_ylim(cy - CENTRE_SPAN/2, cy + CENTRE_SPAN/2)

Strategy (a) shows where each region is. Strategy (b) shows each region at a fixed scale, so a reader can compare sizes across maps while still seeing detail. Fitting each region to the frame — the default — does neither.

Context — a region alone on white gives no sense of place:

gdf.plot(ax=ax, facecolor="#f1f5f9", edgecolor="#e2e8f0", linewidth=0.3, zorder=1)
highlight.plot(ax=ax, facecolor="#ef4444", alpha=0.6, zorder=2)

4. Handle failures per region

results = []
for row in gdf.itertuples():
    try:
        path = render_one(row, gdf, out_dir)
        results.append({"name": row.ward_name, "status": "ok",
                        "path": str(path), "kb": path.stat().st_size // 1024})
    except Exception as exc:
        results.append({"name": row.ward_name, "status": "failed",
                        "error": f"{type(exc).__name__}: {exc}"})
    finally:
        plt.close("all")            # even on failure

plt.close("all") in the finally is what stops a failing region leaking a figure. Without it, a run with 40 failures leaks 40 figures on top of everything else.

The realistic causes of a per-region failure are null or empty geometry, a region with zero area, and a name that cannot be encoded — all of which are worth reporting rather than crashing on. This is the failure policy described in batch script stops at the first bad file.

5. Reuse expensive work across iterations

Anything that does not change per region should happen once. The biggest wins:

# ✅ read once, outside the loop
gdf = gpd.read_file("wards.gpkg").to_crs(27700)

# ✅ fetch one basemap image and reuse the array
import contextily as cx
img, ext = cx.bounds2img(*gdf.total_bounds, zoom=11,
                         source=cx.providers.CartoDB.Positron)

# ✅ dissolve the context layer once, not per region
context = gdf.dissolve().boundary

Fetching a basemap inside the loop means 215 sets of tile requests for the same area — slow for you and abusive to the provider. See how to add a basemap with contextily.

A map for each Greater Manchester ward, in parallel

The data used on this page

Greater Manchester's 215 electoral wards (ONS, December 2024, generalised boundaries) in British National Grid, with each ward's 2021 Census population density from Nomis. The ward names are the real ones, duplicates included. Six of the ten boroughs redrew their wards after the census and Nomis has no 2021 figures for the new ward codes, so 126 wards have no density; the maps draw them in the context colour.

import io
import re
from pathlib import Path

import geopandas as gpd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import pandas as pd
import requests

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
GM = ["Bolton", "Bury", "Manchester", "Oldham", "Rochdale", "Salford", "Stockport", "Tameside", "Trafford", "Wigan"]
r = requests.get(f"{ONS}/Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query", timeout=120, params={
    "where": "LAD24NM IN (" + ",".join(f"'{n}'" for n in GM) + ")",
    "outFields": "WD24CD,WD24NM,LAD24NM", "outSR": 27700, "f": "geojson"})
gdf = gpd.read_file(io.BytesIO(r.content)).set_crs(27700, allow_override=True)
gdf = (gdf.rename(columns={"WD24CD": "ward_code", "WD24NM": "ward_name", "LAD24NM": "borough"})
          .sort_values("ward_code").reset_index(drop=True))
pop = pd.read_csv(io.StringIO(requests.get("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv", timeout=120,
                  params={"date": "latest", "geography": ",".join(gdf["ward_code"]), "c2021_restype_3": 0,
                          "measures": 20100, "select": "geography_code,obs_value"}).text))
gdf["density"] = gdf["ward_code"].map(pop.set_index("GEOGRAPHY_CODE")["OBS_VALUE"]) / (gdf.area / 1e6)
gdf.to_file("wards.gpkg")
print(f"{len(gdf)} wards in {gdf['borough'].nunique()} boroughs; density {gdf['density'].min():,.0f} to {gdf['density'].max():,.0f} per km²")

215 wards in 10 boroughs; density 343 to 13,015 per km²

Example 1: the complete batch renderer

from pathlib import Path
from collections import Counter
import re, time
import geopandas as gpd
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe

def safe(name, max_len=80):
    s = re.sub(r"[^\w\-]+", "_", str(name)).strip("_")
    return (s or "unnamed")[:max_len]

def render_maps(gdf, out_dir, *, name_col, id_col=None, column=None,
                cmap="YlOrRd", span_m=None, dpi=150, figsize=(7, 7),
                overwrite=False, context=None, crs=27700):
    gdf = gdf.to_crs(crs)
    context = (context.to_crs(crs) if context is not None else gdf)
    out_dir = Path(out_dir); out_dir.mkdir(parents=True, exist_ok=True)

    # stable, unique filenames
    base = gdf[name_col].map(safe)
    if id_col is not None and base.duplicated().any():
        base = base + "_" + gdf[id_col].astype(str).map(safe)
    clashes = {n: c for n, c in Counter(base).items() if c > 1}
    if clashes:
        raise ValueError(f"filename collisions remain: {clashes}")
    gdf = gdf.assign(file_name=base.values)      # not "_file": itertuples() renames leading underscores

    # one shared colour range for the whole set
    vmin = vmax = None
    if column:
        vmin, vmax = gdf[column].quantile([0.02, 0.98])

    results, t0 = [], time.perf_counter()
    for row in gdf.itertuples():
        path = out_dir / f"{row.file_name}.png"
        if path.exists() and not overwrite:
            results.append({"name": getattr(row, name_col), "status": "skipped", "file": path.name})
            continue
        try:
            fig, ax = plt.subplots(figsize=figsize)
            context.plot(ax=ax, facecolor="#f1f5f9", edgecolor="#e2e8f0",
                         linewidth=0.3, zorder=1)
            one = gpd.GeoDataFrame(
                [{c: getattr(row, c) for c in [column] if column}],
                geometry=[row.geometry], crs=crs)
            if column:
                one.plot(ax=ax, column=column, cmap=cmap, vmin=vmin, vmax=vmax,
                         edgecolor="#991b1b", linewidth=1.5, zorder=2)
            else:
                one.plot(ax=ax, facecolor="#ef4444", alpha=0.6,
                         edgecolor="#991b1b", linewidth=1.5, zorder=2)

            minx, miny, maxx, maxy = row.geometry.bounds
            cx_, cy_ = (minx + maxx) / 2, (miny + maxy) / 2
            span = span_m or (max(maxx - minx, maxy - miny) * 1.7)
            ax.set_xlim(cx_ - span / 2, cx_ + span / 2)
            ax.set_ylim(cy_ - span / 2, cy_ + span / 2)

            ax.set_title(str(getattr(row, name_col)), fontsize=12, loc="left")
            ax.annotate(f"scale bar: {span/1000:.1f} km across",
                        xy=(0.01, 0.01), xycoords="axes fraction",
                        fontsize=7, color="#64748b")
            ax.set_axis_off()

            tmp = path.with_suffix(".tmp.png")
            fig.savefig(tmp, dpi=dpi, bbox_inches="tight", facecolor="white")
            tmp.replace(path)
            results.append({"name": getattr(row, name_col), "status": "ok",
                            "file": path.name, "kb": path.stat().st_size // 1024})
        except Exception as exc:
            results.append({"name": getattr(row, name_col), "status": "failed",
                            "error": f"{type(exc).__name__}: {exc}"[:120]})
        finally:
            plt.close("all")

    df = pd.DataFrame(results)
    print(df["status"].value_counts().to_dict())
    for r in results:
        if r["status"] == "failed":
            print(f"  ✗ {r['name']}: {r['error']}")
    print(f"{len(df)} regions in {time.perf_counter() - t0:.1f} s")
    return df

report = render_maps(gpd.read_file("wards.gpkg"), "out/wards",
                     name_col="ward_name", id_col="ward_code",
                     column="density", span_m=8_000, dpi=100)
report.to_csv("out/render_report.csv", index=False)

print("second run:")
again = render_maps(gpd.read_file("wards.gpkg"), "out/wards",
                    name_col="ward_name", id_col="ward_code",
                    column="density", span_m=8_000, dpi=100)

{'ok': 215}
215 regions in 19.2 s
second run:
{'skipped': 215}
215 regions in 0.0 s

The column is file_name rather than the more private-looking _file for a reason found by running this: itertuples() builds namedtuples, which cannot have fields starting with an underscore, so pandas silently renames _file to a positional name and row._file raises AttributeError.

Four properties make this a job rather than a script. Filename collisions raise before any rendering starts, so you find out in a second rather than after the whole run. Existing files are skipped unless overwrite=True, which turns a restart into a resume. Each PNG is written to a temporary path and renamed on success, so an interruption never leaves a half-written image that looks complete. And a fixed span_m gives every map the same scale, so the set is comparable — the requirement that separates batch cartography from a loop around savefig.

Example 2: rendering in parallel

Rendering is CPU-bound and each map is independent, so it should parallelise. A worker in a separate process must be importable by that process, so it goes in a module:

%%writefile mapworker.py
import re
from concurrent.futures import ProcessPoolExecutor, as_completed
from pathlib import Path

import geopandas as gpd
import matplotlib

matplotlib.use("Agg")                     # a fresh process does not inherit the parent's backend
import matplotlib.pyplot as plt
from shapely import wkb as shapely_wkb


def safe(name, max_len=80):
    s = re.sub(r"[^\w\-]+", "_", str(name)).strip("_")
    return (s or "unnamed")[:max_len]


def render_worker(args):
    """Runs in a separate process — re-reads only the context it needs."""
    wkb, name, out_path, gpkg, crs, span = args
    try:
        geom = shapely_wkb.loads(wkb)
        c = geom.centroid
        context = gpd.read_file(gpkg, bbox=(c.x - span, c.y - span, c.x + span, c.y + span)).to_crs(crs)
        fig, ax = plt.subplots(figsize=(7, 7))
        context.plot(ax=ax, facecolor="#f1f5f9", edgecolor="#e2e8f0", linewidth=0.3)
        gpd.GeoSeries([geom], crs=crs).plot(ax=ax, facecolor="#ef4444", alpha=0.6,
                                            edgecolor="#991b1b", linewidth=1.5)
        ax.set_xlim(c.x - span / 2, c.x + span / 2)
        ax.set_ylim(c.y - span / 2, c.y + span / 2)
        ax.set_title(name, fontsize=12)
        ax.set_axis_off()
        fig.savefig(out_path, dpi=100, bbox_inches="tight", facecolor="white")
        return {"name": name, "status": "ok"}
    except Exception as exc:
        return {"name": name, "status": "failed", "error": str(exc)[:120]}
    finally:
        plt.close("all")


def render_parallel(gpkg, out_dir, name_col, id_col, workers=4, span=8_000, crs=27700):
    gdf = gpd.read_file(gpkg).to_crs(crs)
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    jobs = [(row.geometry.wkb, str(getattr(row, name_col)),
             str(out_dir / f"{safe(getattr(row, name_col))}_{getattr(row, id_col)}.png"), gpkg, crs, span)
            for row in gdf.itertuples() if row.geometry is not None]
    with ProcessPoolExecutor(max_workers=workers) as pool:
        futures = [pool.submit(render_worker, job) for job in jobs]
        return [f.result() for f in as_completed(futures)]

import os
import time

from mapworker import render_parallel

for workers in (1, 4):
    start = time.perf_counter()
    results = render_parallel("wards.gpkg", f"out/parallel_{workers}", "ward_name", "ward_code", workers=workers)
    ok = sum(r["status"] == "ok" for r in results)
    print(f"{workers} worker(s): {ok} maps in {time.perf_counter() - start:5.1f} s")
print(f"({os.cpu_count()} CPUs on this machine)")

1 worker(s): 215 maps in  15.5 s
4 worker(s): 215 maps in   4.3 s
(24 CPUs on this machine)

Geometries cross the process boundary as WKB bytes, not as Shapely objects. Shapely geometries pickle, but WKB is smaller and faster to move, and it avoids version-mismatch surprises between parent and child interpreters.

Each worker re-reads only the context it needs, using a bounding-box read so it does not load the whole layer per process. Six workers loading a 2 GB layer each is a much worse problem than the serialisation it avoids — the memory reasoning is the same as in parallel batch processing. The speed-up is real but short of the worker count — 3.3× from four processes — as it usually is: each process pays for its own imports, and every map re-reads its context from the GeoPackage.

Example 3: an index page so the output is usable

215 PNGs in a folder is a delivery problem, not a delivery:

from pathlib import Path
import html

def build_index(report_df, out_dir, title="Ward maps"):
    out_dir = Path(out_dir)
    ok = report_df[report_df["status"].isin(["ok", "skipped"])].sort_values("name")
    cards = "\n".join(
        f'<figure><a href="{html.escape(r.file)}">'
        f'<img src="{html.escape(r.file)}" loading="lazy" alt="Map of {html.escape(str(r.name))}">'
        f'</a><figcaption>{html.escape(str(r.name))}</figcaption></figure>'
        for r in ok.itertuples())
    failed = report_df[report_df["status"] == "failed"]
    notes = ("<p class='warn'>Not rendered: "
             + ", ".join(html.escape(str(n)) for n in failed["name"]) + "</p>"
             if len(failed) else "")

    (out_dir / "index.html").write_text(f"""<!doctype html>
<html lang="en"><head><meta charset="utf-8">
<title>{html.escape(title)}</title>
<style>
 body {{ font: 15px/1.5 system-ui, sans-serif; margin: 2rem; color: #1e293b; }}
 .grid {{ display: grid; gap: 1rem;
          grid-template-columns: repeat(auto-fill, minmax(220px, 1fr)); }}
 figure {{ margin: 0; }}
 img {{ width: 100%; border: 1px solid #e2e8f0; border-radius: 6px; }}
 figcaption {{ font-size: 13px; color: #475569; margin-top: .3rem; }}
 .warn {{ color: #b91c1c; }}
</style></head>
<body>
<h1>{html.escape(title)}</h1>
<p>{len(ok):,} maps · all at the same scale and colour range</p>
{notes}
<div class="grid">{cards}</div>
</body></html>""", encoding="utf-8")
    return out_dir / "index.html"

index = build_index(report, "out/wards")
links = re.findall(r'<img src="([^"]+)"', index.read_text())
print(f"{index}: {len(links)} maps, {sum((index.parent / html.unescape(l)).exists() for l in links)} links resolve")

out/wards/index.html: 215 maps, 215 links resolve

The links come from the file column of the render report, not from re-deriving a name. An earlier version of this page rebuilt each filename with safe(name) — but render_maps had appended the ward code to every file as soon as it found a duplicate name, so every link on the page pointed at a file that did not exist. Passing the name the renderer actually used is the only way the two cannot drift apart.

html.escape on every interpolated value is not optional. Ward names contain apostrophes and ampersands — Greater Manchester alone has St Mary's, Ashton St Michael's and Bramhall South & Woodford — and they produce broken markup without it — and worse if the names ever come from an untrusted source.

Listing the failures on the page is what makes the index a report rather than a gallery. A reader who cannot find their ward should be able to see that it failed, not conclude the job was incomplete.

A batch job and a cartography job at once

Batch cartography sits at the intersection of two disciplines, and it fails when either half is neglected.

As a batch job, it has the ordinary requirements: bounded resources, error isolation, resumability, and a report that distinguishes outcomes. The specific resource here is matplotlib figures. plt keeps a global registry so that plt.show() can find every open figure, which means a figure is not garbage-collected when it goes out of scope — it is still referenced. Step 1 measured about 1 MB of traced memory per figure for a simple ward map; a figure with a basemap image or thousands of features holds far more, and a few hundred of those is a crash. plt.close(fig) removes it from the registry, and plt.close("all") in a finally covers the failure path too. This is the same class of bug as too many open files: a resource acquired in a loop and released never.

As cartography, it has one requirement no single map has: the outputs must be comparable to each other. Every matplotlib default works against this, and reasonably so — autoscaling to fit the data is right for one plot and wrong for a set, because it removes exactly the information a reader would use to compare. A ward of 2 km² and one of 40 km² both fill their frame, and the two maps are individually correct and collectively misleading.

The same applies to colour. vmin and vmax default to each layer's own range, so a region whose values span £20k–£25k gets the full colour ramp and looks as varied as one spanning £20k–£200k. Fixing the range across the set is what makes a colour mean the same thing everywhere, and it is the batch equivalent of using fixed classification breaks across a time series.

The third neglected element is context. A region rendered alone is a shape on white, and a reader who does not already know the area cannot place it. Drawing the whole layer as a pale background and the subject region highlighted costs one extra plot() call and is the difference between a diagram and a map.

Finally, note where the time goes. Rendering dominates — the 215 maps above took about 70 ms each on one process — and it is CPU-bound Python and C code that holds the GIL, so processes beat threads: four workers cut the run from 15 s to 4.5 s. But before reaching for parallelism, move the fixed costs out of the loop: reading the layer once instead of 215 times, fetching one basemap image instead of 215 sets of tiles, dissolving the context once. Those changes are usually larger than the speed-up from six workers, and they do not introduce the memory multiplication that parallel processing brings with it.

Backend order, closing figures and stable filenames

• matplotlib.use("Agg") must precede import matplotlib.pyplot, and must be repeated in each worker process.
• plt.close(fig) every iteration, and plt.close("all") in a finally so failures do not leak.
• A figure is roughly 15–20 MB at 150 dpi; matplotlib warns after 20 are open, which is easy to miss.
• Filenames must be stable across runs. A counter-based suffix breaks the correspondence between yesterday's outputs and today's.
• bbox_inches="tight" changes the output size per map, so images in a set have different pixel dimensions. Drop it if uniform size matters.
• facecolor="white" on savefig, or a transparent background renders black in some viewers.
• Null and empty geometries raise on .bounds. Filter them before the loop, and report them.
• dpi multiplies file size quadratically. 300 dpi is four times the bytes of 150.
• Fetch the basemap once with cx.bounds2img and reuse the array; per-map tile requests are slow and abusive.
• Deliver an index page, not a folder. A gallery with the failures listed is a report; a folder is homework.

Internal links

• How to plot multiple layers on one map in GeoPandas — the composition being repeated

• How to add legends, labels and a scale bar to a GeoPandas map — the furniture each map needs

• How to add a basemap to a GeoPandas map with contextily — fetching once, not per map

• How to save a map as an image in Python with Matplotlib — export settings in detail

• How to batch process a folder of GIS files in Python — the batch pattern this follows

• How to speed up batch GIS jobs with parallel processing — when rendering dominates

• Choropleth classification explained — shared breaks across a set

• How to package GIS deliverables into zipped bundles — shipping hundreds of images

FAQ

Why does my script run out of memory after a few hundred maps?

Every figure stays in matplotlib's global registry until closed, with all the memory it holds. Call plt.close(fig) at the end of every iteration, and plt.close("all") in a finally.

Do I need matplotlib.use("Agg")?

Yes for headless environments — servers, containers, cron — and it is faster than an interactive backend everywhere. Set it before importing pyplot, and again inside worker processes.

Why do my maps all look the same size?

Matplotlib autoscales each figure to its data. Set explicit limits, either the same extent for every map or a fixed span centred on each region, so sizes are comparable.

How do I make colours mean the same thing across maps?

Compute vmin and vmax once from the whole dataset and pass them to every plot() call. Clipping at the 2nd and 98th percentiles stops one outlier flattening the range.

How do I stop two regions overwriting each other's file?

Sanitise the names, then check for collisions before rendering anything and disambiguate with a stable id column. Raising early beats discovering it after eighty maps.

Should I render in parallel?

Only after moving fixed costs out of the loop — one layer read, one basemap fetch. Then processes rather than threads, since rendering is CPU-bound.

How should I deliver hundreds of images?

With an index page that shows thumbnails and lists anything that failed to render. A folder of PNGs makes the recipient do the work of finding out what is in it.