Source code: Make a map in Python with matplotlib and cartopy: projections, classified choropleths, scale bars, north arrows and inset maps

"""Make publication-quality maps in Python with matplotlib and cartopy.

Natural Earth for world and context layers, ONS local authority boundaries with Census 2021 population density
for a classified choropleth of England and Wales, and the Copernicus DEM for a shaded-relief map of Greater
Manchester. Covers projections, classification and legends, scale bars, north arrows, inset locator maps,
gridlines, light and dark styling and print export.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/python-maps-matplotlib-cartopy/
"""
import io
import os
import sys
import warnings
from importlib.metadata import version
from pathlib import Path

for var in ("PROJ_LIB", "PROJ_DATA", "GDAL_DATA"):   # stray installs (e.g. PostGIS) break the wheels' own copies
    os.environ.pop(var, None)

import numpy as np
import pandas as pd
import geopandas as gpd
import requests
import rasterio
import rasterio.merge
import rasterio.features
from rasterio.warp import reproject, Resampling, calculate_default_transform
import pystac_client
import planetary_computer as pc
import mapclassify
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io import shapereader
import matplotlib
matplotlib.use("Agg")   # render to files only
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe
from matplotlib.colors import LightSource, ListedColormap, Normalize
from matplotlib.patches import Patch, Rectangle
from matplotlib_scalebar.scalebar import ScaleBar
from shapely.geometry import box, Polygon
from pyproj import Geod
from PIL import Image

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=cartopy.io.DownloadWarning)
DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)
cartopy.config["data_dir"] = str(Path("cache") / "cartopy")    # Natural Earth downloads land here

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Local_Authority_Districts_December_2022_UK_BGC_V2/FeatureServer/0/query")
NOMIS = "https://www.nomisweb.co.uk/api/v01/dataset/NM_2026_1.data.csv"   # Census 2021 TS006 population density
STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
OSGB = ccrs.OSGB(approx=False)        # British National Grid, EPSG:27700
MM = 1 / 25.4                         # millimetres to inches

# Two themes: the same map must read on a white page and on a dark screen.
THEMES = {
    "light": dict(bg="white", sea="#e8eef2", land="#d9d9d4", edge="white", text="#1a1a1a", muted="#555555",
                  grid="#8a8a8a", frame="#333333", halo="white"),
    "dark": dict(bg="#0d1012", sea="#161b1e", land="#2a3033", edge="#0d1012", text="#e7ebe9", muted="#a0aaa7",
                 grid="#6d7774", frame="#a0aaa7", halo="#0d1012"),
}
SEQ = ListedColormap(plt.get_cmap("Spectral_r")(np.linspace(0.42, 1.0, 256)))   # yellow -> red: sequential


def section(title):
    print(f"\n=== {title} ===")


section("versions")
print(f"python {sys.version.split()[0]} | cartopy {cartopy.__version__} | matplotlib {matplotlib.__version__} | "
      f"geopandas {gpd.__version__} | mapclassify {mapclassify.__version__} | "
      f"matplotlib-scalebar {version('matplotlib-scalebar')} | rasterio {rasterio.__version__}")


def ne(scale, category, name):
    return gpd.read_file(shapereader.natural_earth(resolution=scale, category=category, name=name))


# ------------------------------------------------------------------ Step 1: projections
section("Step 1: Natural Earth and four projections")
land110 = ne("110m", "physical", "land")
print(f"Natural Earth 110m land: {len(land110)} polygons, CRS {land110.crs.to_string()}")
projs = [("PlateCarree (equirectangular)", ccrs.PlateCarree()),
         ("Mercator", ccrs.Mercator(min_latitude=-80, max_latitude=84)),
         ("Equal Earth (equal area)", ccrs.EqualEarth()),
         ("Orthographic over the UK", ccrs.Orthographic(central_longitude=-2, central_latitude=54))]
fig = plt.figure(figsize=(170 * MM, 112 * MM), layout="constrained")
for i, (name, proj) in enumerate(projs):
    ax = fig.add_subplot(2, 2, i + 1, projection=proj)
    ax.set_global()
    ax.add_feature(cfeature.OCEAN.with_scale("110m"), facecolor="#e8eef2")
    ax.add_geometries(land110.geometry, crs=ccrs.PlateCarree(), facecolor="#d9d9d4", edgecolor="#8a8a8a", lw=0.3)
    # Tissot indicatrices: identical 500 km circles on the ground, drawn in each projection
    ax.tissot(rad_km=500, lons=np.arange(-150, 180, 60), lats=np.arange(-60, 80, 30), n_samples=40,
              facecolor="#d53e4f", alpha=0.6, edgecolor="none")
    ax.gridlines(color="#8a8a8a", lw=0.3, alpha=0.6)
    ax.set_title(name, fontsize=8.5)
fig.savefig(OUT / "01-map-projections-tissot.png", dpi=200, facecolor="white")
plt.close(fig)

# Distortion in numbers: the area of a 1 degree x 1 degree cell at different latitudes in each projection
rows = []
for lat in (0, 30, 53.5, 70):
    cell = gpd.GeoSeries([Polygon([(0, lat), (1, lat), (1, lat + 1), (0, lat + 1)])], crs=4326)
    true_area = cell.to_crs("+proj=laea +lat_0=%s +lon_0=0.5" % (lat + 0.5)).area.iloc[0]
    row = {"lat": lat, "true_km2": true_area / 1e6}
    # PlateCarree's native units are degrees, so measure it as the same projection in metres (+proj=eqc)
    for name, crs in [("PlateCarree", "+proj=eqc +R=6371000"), ("Mercator", projs[1][1].proj4_init),
                      ("Equal Earth", projs[2][1].proj4_init)]:
        row[name] = cell.to_crs(crs).area.iloc[0] / true_area
    rows.append(row)
print("area of a 1 x 1 degree cell in the projection / true area:")
print(pd.DataFrame(rows).round(3).to_string(index=False))

# ------------------------------------------------------------------ Step 2: world choropleth
section("Step 2: world choropleth of population density (Natural Earth 50m, Equal Earth)")
world = ne("50m", "cultural", "admin_0_countries")
world = world[world.POP_EST > 0].copy()
world["area_km2"] = world.to_crs("+proj=eqearth").area / 1e6
world["density"] = world.POP_EST / world.area_km2
print(f"{len(world)} countries with a population estimate; POP_YEAR values: {sorted(world.POP_YEAR.unique())[-3:]}")
q = mapclassify.Quantiles(world.density, k=6)
print("quantile breaks (people per km2):", np.round(q.bins, 1).tolist())
print(world.nlargest(3, "density")[["NAME", "density"]].round(0).to_string(index=False))
print("United Kingdom:", round(float(world.loc[world.ISO_A3 == "GBR", "density"].iloc[0]), 1), "per km2")

fig = plt.figure(figsize=(170 * MM, 95 * MM), layout="constrained")
ax = fig.add_subplot(projection=ccrs.EqualEarth())
ax.set_global()
ax.set_facecolor("#e8eef2")
cols = SEQ(np.linspace(0, 1, q.k))
world["cls"] = q.yb
for k in range(q.k):
    ax.add_geometries(world.geometry[world.cls == k], crs=ccrs.PlateCarree(), facecolor=cols[k],
                      edgecolor="white", lw=0.2)
lows = np.r_[world.density.min(), q.bins[:-1]]
labels = [f"{lo:,.0f} – {hi:,.0f}  ({n})" for lo, hi, n in zip(lows, q.bins, q.counts)]
ax.legend(handles=[Patch(fc=c, ec="none", label=l) for c, l in zip(cols, labels)], loc="lower left",
          title="People per km² (countries)", fontsize=6.5, title_fontsize=7, frameon=True, framealpha=0.85,
          facecolor="white", edgecolor="none", bbox_to_anchor=(-0.01, 0.0))
ax.gridlines(color="#8a8a8a", lw=0.3, alpha=0.6)
ax.set_title("Population density by country, quantile classes (Natural Earth 1:50m, Equal Earth projection)",
             fontsize=8.5)
fig.savefig(OUT / "02-world-population-density.png", dpi=200, facecolor="white")
plt.close(fig)

# ------------------------------------------------------------------ Step 3: ONS boundaries + census
section("Step 3: England and Wales local authorities and Census 2021 density")
lad_path = DATA / "lad_dec2022_ew_bgc.gpkg"
if not lad_path.exists():
    r = requests.get(ONS, params={"where": "LAD22CD LIKE 'E%' OR LAD22CD LIKE 'W%'", "outFields": "LAD22CD,LAD22NM",
                                  "outSR": 27700, "f": "geojson", "resultRecordCount": 1000}, timeout=300)
    r.raise_for_status()
    gpd.GeoDataFrame.from_features(r.json()["features"], crs=27700).to_file(lad_path)
lad = gpd.read_file(lad_path)
r = requests.get(NOMIS, params={"date": "latest", "geography": "TYPE154", "measures": 20100,
                                "select": "geography_code,obs_value"}, timeout=300)
r.raise_for_status()
dens = pd.read_csv(io.StringIO(r.text)).rename(columns={"GEOGRAPHY_CODE": "LAD22CD", "OBS_VALUE": "density"})
lad = lad.merge(dens, on="LAD22CD", how="left")
print(f"{len(lad)} local authority districts; {lad.density.isna().sum()} without a census value")
print(f"density: median {lad.density.median():,.0f}, min {lad.density.min():,.1f} "
      f"({lad.loc[lad.density.idxmin(), 'LAD22NM']}), max {lad.density.max():,.0f} "
      f"({lad.loc[lad.density.idxmax(), 'LAD22NM']}) usual residents per km2")

schemes = {"Equal interval": mapclassify.EqualInterval(lad.density, k=6),
           "Quantiles": mapclassify.Quantiles(lad.density, k=6),
           "Fisher-Jenks (natural breaks)": mapclassify.FisherJenks(lad.density, k=6),
           "Manual (rounded, roughly doubling)": mapclassify.UserDefined(lad.density,
                                                                        bins=[100, 250, 500, 1000, 2500, 5000])}
print("class breaks and counts:")
for name, c in schemes.items():
    print(f"  {name:<35} breaks {np.round(c.bins).astype(int).tolist()}  counts {c.counts.tolist()}  "
          f"GADF {c.gadf:.3f}")

fig, axs = plt.subplots(1, 3, figsize=(170 * MM, 82 * MM), subplot_kw={"projection": OSGB}, layout="constrained")
for ax, name in zip(axs, list(schemes)[:3]):
    c = schemes[name]
    cc = SEQ(np.linspace(0, 1, c.k))
    ax.set_extent([130_000, 660_000, 0, 660_000], crs=OSGB)
    for k in range(c.k):
        ax.add_geometries(lad.geometry[c.yb == k], crs=OSGB, facecolor=cc[k], edgecolor="white", lw=0.1)
    lows = np.r_[lad.density.min(), c.bins[:-1]]
    ax.legend(handles=[Patch(fc=col, label=f"{lo:,.0f}–{hi:,.0f} ({n})") for col, lo, hi, n
                       in zip(cc, lows, c.bins, c.counts)], loc="upper left", fontsize=5.2, frameon=False,
              handlelength=1, handletextpad=0.4, borderaxespad=0.2)
    ax.set_title(name, fontsize=7.5)
    ax.axis("off")
fig.savefig(OUT / "03-classification-schemes.png", dpi=200, facecolor="white")
plt.close(fig)


# ------------------------------------------------------------------ Step 4: the publication map
def north_arrow(ax, x, y, theme, size=0.06):
    """Arrow to true north at axes position (x, y): rotate by the grid convergence at that point."""
    px, py = ax.transAxes.transform((x, y))
    mx, my = ax.transData.inverted().transform((px, py))
    lon, lat = ccrs.PlateCarree().transform_point(mx, my, ax.projection)
    nx, ny = ax.projection.transform_point(lon, lat + 0.1, ccrs.PlateCarree())
    angle = np.degrees(np.arctan2(nx - mx, ny - my))           # clockwise from grid north
    dx, dy = size * np.sin(np.radians(angle)), size * np.cos(np.radians(angle))
    ax.annotate("", xy=(x + dx, y + dy), xytext=(x, y), xycoords="axes fraction",
                arrowprops=dict(arrowstyle="-|>,head_width=0.35,head_length=0.7", color=theme["text"], lw=1.2))
    ax.text(x + dx * 1.35, y + dy * 1.35, "N", transform=ax.transAxes, ha="center", va="center",
            fontsize=8, weight="bold", color=theme["text"])
    return angle


def ew_map(theme_name, path_stem, dpi):
    t = THEMES[theme_name]
    c = schemes["Manual (rounded, roughly doubling)"]
    cc = SEQ(np.linspace(0, 1, c.k + 1))[1:] if theme_name == "dark" else SEQ(np.linspace(0, 1, c.k + 1))[:-1]
    fig = plt.figure(figsize=(120 * MM, 140 * MM), facecolor=t["bg"])
    ax = fig.add_axes([0.075, 0.075, 0.905, 0.845], projection=OSGB)
    ax.set_extent([60_000, 690_000, -10_000, 680_000], crs=OSGB)
    ax.set_facecolor(t["sea"])
    # zorder matters: context land first (1), then one collection per class on top (2)
    ax.add_geometries(context.geometry, crs=ccrs.PlateCarree(), facecolor=t["land"], edgecolor=t["edge"], lw=0.3,
                      zorder=1)
    for k in range(c.k):
        ax.add_geometries(lad.geometry[c.yb == k], crs=OSGB, facecolor=cc[k], edgecolor=t["edge"], lw=0.15,
                          zorder=2)
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, xlocs=range(-8, 4, 2), ylocs=range(50, 57),
                      color=t["grid"], lw=0.3, alpha=0.7, dms=True, x_inline=False, y_inline=False)
    gl.top_labels = gl.right_labels = False
    gl.xlabel_style = gl.ylabel_style = {"size": 6.5, "color": t["muted"]}
    lows = np.r_[0, c.bins[:-1]]
    labels = [f"{lo:,.0f} – {hi:,.0f}" if i < c.k - 1 else f"{lo:,.0f} – {lad.density.max():,.0f}"
              for i, (lo, hi) in enumerate(zip(lows, c.bins))]
    labels = [f"{l}   ({n})" for l, n in zip(labels, c.counts)]
    leg = ax.legend(handles=[Patch(fc=col, ec=t["edge"], lw=0.3, label=l) for col, l in zip(cc[::-1], labels[::-1])],
                    loc="upper right", bbox_to_anchor=(0.995, 0.995), fontsize=6.5, title_fontsize=7,
                    title="Usual residents per km²\n(number of districts)", frameon=True, framealpha=0.92,
                    facecolor=t["bg"], edgecolor=t["frame"], labelcolor=t["text"])
    leg.get_title().set_color(t["text"]); leg.get_frame().set_linewidth(0.4)
    # 1 map unit = 1 m in British National Grid; fixed_value gives a round 100 km bar
    ax.add_artist(ScaleBar(1, "m", location="lower left", fixed_value=100, fixed_units="km", box_alpha=0,
                           color=t["text"], font_properties={"size": 7}, scale_loc="top"))
    angle = north_arrow(ax, 0.07, 0.17, t)
    # Inset locator: where England and Wales sit in Europe
    ins = fig.add_axes([0.09, 0.665, 0.24, 0.24],
                       projection=ccrs.LambertAzimuthalEqualArea(central_longitude=10, central_latitude=52))
    ins.set_extent([-12, 32, 34, 68], crs=ccrs.PlateCarree())
    ins.set_facecolor(t["sea"])
    ins.add_geometries(europe.geometry, crs=ccrs.PlateCarree(), facecolor=t["land"], edgecolor=t["edge"], lw=0.2)
    ins.add_geometries(ew_outline, crs=OSGB, facecolor="#d53e4f", edgecolor="none")
    ins.add_geometries([box(60_000, -10_000, 690_000, 680_000)], crs=OSGB, facecolor="none",
                       edgecolor=t["text"], lw=0.8)
    ins.spines["geo"].set_edgecolor(t["frame"]); ins.spines["geo"].set_linewidth(0.6)
    ax.spines["geo"].set_edgecolor(t["frame"]); ax.spines["geo"].set_linewidth(0.6)
    fig.text(0.02, 0.975, "Population density, England and Wales, Census 2021", fontsize=10, weight="bold",
             color=t["text"], va="top")
    fig.text(0.02, 0.945, "Local authority districts (331), British National Grid", fontsize=7.5,
             color=t["muted"], va="top")
    fig.text(0.98, 0.006, "Source: ONS Census 2021 (TS006), ONS LAD Dec 2022 BGC, OGL v3.0. Contains OS data "
             "© Crown copyright and database right 2022.\nNatural Earth (public domain). Made with cartopy.",
             ha="right", va="bottom", fontsize=4.8, color=t["muted"])
    paths = {}
    for ext in ("png", "pdf", "svg") if theme_name == "light" else ("png",):
        p = OUT / f"{path_stem}.{ext}"
        fig.savefig(p, dpi=dpi, facecolor=t["bg"])
        paths[ext] = p
    plt.close(fig)
    return angle, paths


section("Step 4: a publication-quality map of England and Wales")
countries10 = ne("10m", "cultural", "admin_0_countries")
context = countries10[countries10.ADMIN.isin(["United Kingdom", "Ireland", "France", "Isle of Man", "Belgium",
                                              "Netherlands", "Guernsey", "Jersey"])]
europe = ne("50m", "cultural", "admin_0_countries").cx[-30:45, 30:72]
ew_outline = [lad.union_all().simplify(2000)]
angle, paths = ew_map("light", "04-england-wales-density-map", dpi=300)
print(f"grid convergence at the north arrow: true north is {angle:+.2f} degrees from grid north")
for ext, p in paths.items():
    extra = ""
    if ext == "png":
        with Image.open(p) as im:
            extra = f", {im.width} x {im.height} px"
    print(f"  {p.name}: {p.stat().st_size / 1024:,.0f} KB{extra}")
_, dark_paths = ew_map("dark", "06-england-wales-density-map-dark", dpi=200)

# Scale check: what the scale bar assumes vs the true ground distance at different eastings
for e in (100_000, 400_000, 650_000):
    lon1, lat1 = ccrs.PlateCarree().transform_point(e, 300_000, OSGB)
    lon2, lat2 = ccrs.PlateCarree().transform_point(e + 100_000, 300_000, OSGB)
    d = Geod(ellps="WGS84").inv(lon1, lat1, lon2, lat2)[2]
    print(f"  100 km of grid at easting {e:,} m = {d / 1000:.2f} km on the ground")

# ------------------------------------------------------------------ Step 5: raster theme
section("Step 5: shaded relief of Greater Manchester (Copernicus DEM GLO-30)")
gm = lad[lad.LAD22CD.str.match(r"E0800000\d|E08000010")].copy()
print(f"Greater Manchester districts: {len(gm)}")
xmin, ymin, xmax, ymax = gm.total_bounds + np.array([-6000, -6000, 6000, 6000])
ll = gpd.GeoSeries([box(xmin, ymin, xmax, ymax)], crs=27700).to_crs(4326).total_bounds
dem_path = DATA / "gm_dem_bng_30m.tif"
if not dem_path.exists():
    cat = pystac_client.Client.open(STAC_URL, modifier=pc.sign_inplace)
    items = list(cat.search(collections=["cop-dem-glo-30"], bbox=list(ll)).items())
    print("DEM tiles:", sorted(i.id for i in items))
    srcs = [rasterio.open(i.assets["data"].href) for i in items]
    mosaic, mt = rasterio.merge.merge(srcs, bounds=tuple(ll))
    src_crs = srcs[0].crs
    for s in srcs:
        s.close()
    width = int((xmax - xmin) / 30); height = int((ymax - ymin) / 30)
    dst_t = rasterio.transform.from_origin(xmin, ymax, 30, 30)
    out = np.full((height, width), np.nan, dtype="float32")
    reproject(mosaic[0].astype("float32"), out, src_transform=mt, src_crs=src_crs, dst_transform=dst_t,
              dst_crs="EPSG:27700", resampling=Resampling.bilinear)
    with rasterio.open(dem_path, "w", driver="GTiff", width=width, height=height, count=1, dtype="float32",
                       crs="EPSG:27700", transform=dst_t, nodata=np.nan, compress="deflate") as dst:
        dst.write(out, 1)
with rasterio.open(dem_path) as src:
    dem = src.read(1)
    dext = [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top]
print(f"DEM {dem.shape[1]} x {dem.shape[0]} px at 30 m; elevation {np.nanmin(dem):.0f} to {np.nanmax(dem):.0f} m")
inside_gm = rasterio.features.geometry_mask(gm.geometry, dem.shape,
                                            rasterio.transform.from_origin(dext[0], dext[3], 30, 30), invert=True)
print(f"inside Greater Manchester: {np.nanmin(dem[inside_gm]):.0f} to {np.nanmax(dem[inside_gm]):.0f} m, "
      f"mean {np.nanmean(dem[inside_gm]):.0f} m")

t = THEMES["light"]
fig = plt.figure(figsize=(170 * MM, 130 * MM), facecolor="white")
ax = fig.add_axes([0.065, 0.05, 0.775, 0.86], projection=OSGB)
ax.set_extent(dext, crs=OSGB)
norm = Normalize(0, 600)
ls = LightSource(azdeg=315, altdeg=45)
rgb = ls.shade(np.nan_to_num(dem, nan=0), cmap=plt.get_cmap("Spectral_r"), norm=norm, vert_exag=2,
               blend_mode="soft", dx=30, dy=30)
ax.imshow(rgb, extent=dext, transform=OSGB, origin="upper", interpolation="bilinear")
ax.add_geometries(gm.geometry, crs=OSGB, facecolor="none", edgecolor="#1a1a1a", lw=0.7)
for _, rw in gm.iterrows():
    p = rw.geometry.representative_point()
    ax.text(p.x, p.y, rw.LAD22NM, transform=OSGB, ha="center", va="center", fontsize=7, color="#111111",
            path_effects=[pe.withStroke(linewidth=2.2, foreground="white")])
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, color="#555555", lw=0.3, alpha=0.6,
                  xlocs=np.arange(-3, -1.5, 0.2), ylocs=np.arange(53.2, 53.9, 0.1), x_inline=False, y_inline=False)
gl.top_labels = gl.right_labels = False
gl.xlabel_style = gl.ylabel_style = {"size": 6.5}
ax.add_artist(ScaleBar(1, "m", location="lower right", length_fraction=0.2, box_alpha=0.75, font_properties={"size": 7}))
north_arrow(ax, 0.94, 0.86, t, size=0.07)
cax = fig.add_axes([0.86, 0.25, 0.02, 0.45])
cb = fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap="Spectral_r"), cax=cax, extend="max")
cb.set_label("Elevation (m)", fontsize=7.5); cb.ax.tick_params(labelsize=6.5)
ins = fig.add_axes([0.855, 0.74, 0.13, 0.2], projection=OSGB)
ins.set_extent([-30_000, 680_000, -20_000, 1_230_000], crs=OSGB)
ins.add_geometries(countries10[countries10.ADMIN.isin(["United Kingdom", "Ireland", "Isle of Man"])].geometry,
                   crs=ccrs.PlateCarree(), facecolor="#d9d9d4", edgecolor="#8a8a8a", lw=0.2)
ins.add_geometries([box(dext[0], dext[2], dext[1], dext[3])],
                   crs=OSGB, facecolor="#d53e4f", edgecolor="#d53e4f", lw=0.8)
ins.set_facecolor("#e8eef2")
fig.text(0.04, 0.975, "Greater Manchester: terrain and the ten metropolitan boroughs", fontsize=10, weight="bold",
         va="top")
ax.text(0.01, 0.99, "Copernicus DEM GLO-30 © DLR, Airbus (ESA, EU); ONS LAD Dec 2022, OGL v3.0, "
        "contains OS data © Crown copyright 2022", transform=ax.transAxes, fontsize=5, color="#333333", va="top",
        bbox=dict(fc="white", ec="none", alpha=0.7, pad=1))
fig.savefig(OUT / "05-greater-manchester-terrain.png", dpi=200, facecolor="white")
plt.close(fig)
print("saved figures:", sorted(p.name for p in OUT.glob("0*.png")))