How to Batch-Generate a Map Image for Every Region in Python
Problem statement
Someone needs one PNG per ward. There are 340 wards. 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 until the process dies around ward 180, because every figure is still open. Two wards are called St. Mary's and St Mary's, so one file overwrites the other. 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["income"].min(), gdf["income"].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=150,
bbox_inches="tight", facecolor="white")
plt.close(fig) # ← the line that prevents the crash
| Requirement | Line |
|---|---|
| no display, faster rendering | matplotlib.use("Agg") |
| bounded memory | plt.close(fig) every iteration |
| filenames that cannot collide | safe() plus a collision check |
| comparable outputs | shared vmin/vmax, shared extent or shared zoom |
| one bad region does not stop the run | try/except per item |
Step-by-step solution
1. Use the Agg backend and close every figure
import matplotlib
matplotlib.use("Agg") # before importing pyplot
import matplotlib.pyplot as pltAgg 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
tracemalloc.start()
for i, row in enumerate(gdf.head(50).itertuples()):
fig, ax = plt.subplots(figsize=(7, 7))
gdf.plot(ax=ax)
fig.savefig(f"/tmp/{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") 0 figures open: 1 memory 4.2 MB
10 figures open: 11 memory 184.9 MB
20 figures open: 21 memory 361.4 MB
40 figures open: 41 memory 712.8 MB
Roughly 18 MB per figure, growing without limit. With plt.close(fig) the number stays at one. Matplotlib emits a warning after 20 open figures, and it is easy to miss in a busy log.
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: {'St_Mary_s': 2, 'Central': 3}
Three wards called "Central" would produce one file. 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)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["income"].quantile([0.02, 0.98])
part.plot(column="income", 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 failureplt.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().boundaryFetching a basemap inside the loop means 340 sets of tile requests for the same area — slow for you and abusive to the provider. See how to add a basemap with contextily.
Code examples
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=base.values)
# 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}.png"
if path.exists() and not overwrite:
results.append({"name": getattr(row, name_col), "status": "skipped"})
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",
"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="income", span_m=12_000)
report.to_csv("out/render_report.csv", index=False){'ok': 336, 'skipped': 2, 'failed': 2}
✗ Offshore: ValueError: cannot render an empty geometry
✗ Docks Ward: AttributeError: 'NoneType' object has no attribute 'bounds'
2 regions in 84.2 s
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 eighty. 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 parallelises well:
from concurrent.futures import ProcessPoolExecutor, as_completed
from pathlib import Path
import geopandas as gpd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
def render_worker(args):
"""Runs in a separate process — must re-import and re-read what it needs."""
wkb, name, value, bounds, vmin, vmax, out_path, gpkg, crs, span = args
import geopandas as gpd
from shapely import wkb as shapely_wkb
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
try:
geom = shapely_wkb.loads(wkb)
context = gpd.read_file(gpkg, bbox=(
geom.centroid.x - span, geom.centroid.y - span,
geom.centroid.x + span, geom.centroid.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)
cx_, cy_ = geom.centroid.x, geom.centroid.y
ax.set_xlim(cx_ - span / 2, cx_ + span / 2)
ax.set_ylim(cy_ - span / 2, cy_ + span / 2)
ax.set_title(name, fontsize=12)
ax.set_axis_off()
fig.savefig(out_path, dpi=150, 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, workers=6, span=12_000, crs=27700):
gdf = gpd.read_file(gpkg).to_crs(crs)
out_dir = Path(out_dir); out_dir.mkdir(parents=True, exist_ok=True)
vmin = vmax = None
jobs = [(row.geometry.wkb, str(getattr(row, name_col)), None,
row.geometry.bounds, vmin, vmax,
str(out_dir / f"{safe(getattr(row, name_col))}.png"),
gpkg, crs, span)
for row in gdf.itertuples() if row.geometry is not None]
done = 0
results = []
with ProcessPoolExecutor(max_workers=workers) as pool:
futures = {pool.submit(render_worker, j): j[1] for j in jobs}
for fut in as_completed(futures):
results.append(fut.result())
done += 1
print(f" {done}/{len(jobs)}", end="\r")
print()
return resultsGeometries 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.
matplotlib.use("Agg") appears in the worker as well as the parent, because a fresh process does not inherit the parent's backend selection.
Example 3: an index page so the output is usable
340 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"] == "ok"].sort_values("name")
cards = "\n".join(
f'<figure><a href="{html.escape(safe(r.name))}.png">'
f'<img src="{html.escape(safe(r.name))}.png" 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"
build_index(report, "out/wards")html.escape on every interpolated value is not optional. Ward names contain apostrophes and ampersands, and a name like St Mary's & Docks produces 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.
Explanation
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. At roughly 18 MB each, three hundred of them 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 — usually 100–300 ms per figure at 150 dpi — and it is CPU-bound in C code that releases the GIL poorly, so processes beat threads. But before reaching for parallelism, move the fixed costs out of the loop: reading the layer once instead of 340 times, fetching one basemap image instead of 340 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.
Edge cases or notes
matplotlib.use("Agg")must precedeimport matplotlib.pyplot, and must be repeated in each worker process.plt.close(fig)every iteration, andplt.close("all")in afinallyso 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"onsavefig, or a transparent background renders black in some viewers.- Null and empty geometries raise on
.bounds. Filter them before the loop, and report them. dpimultiplies file size quadratically. 300 dpi is four times the bytes of 150.- Fetch the basemap once with
cx.bounds2imgand 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 340 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, at roughly 18 MB each. 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 340 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.