Source code: Coordinate reference systems in Python: EPSG codes, pyproj transforms and reprojecting vectors and rasters

"""Coordinate reference systems in Python: EPSG codes, pyproj transforms and reprojecting vectors and rasters.

Greater Manchester, UK: ONS borough boundaries, 12,000+ NaPTAN bus stops with published lat/lon AND
British National Grid coordinates (a free accuracy test for pyproj), and the Copernicus DEM GLO-30
from the Microsoft Planetary Computer, reprojected from WGS 84 to EPSG:27700.

Part 1 of the TechnicalGIS "Geospatial Python essentials" series:
https://www.technicalgis.com/tutorials/coordinate-reference-systems-python/
"""
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 pyproj
from pyproj import CRS, Transformer, Geod
from pyproj.transformer import TransformerGroup
import rasterio
import rasterio.merge
from rasterio.warp import calculate_default_transform, reproject, transform_bounds
from rasterio.enums import Resampling
import rioxarray  # noqa: F401  (registers the .rio accessor)
import xarray as xr
from shapely.geometry import Point, LineString
import pystac_client
import planetary_computer as pc
import matplotlib
matplotlib.use("Agg")   # render to files only
import matplotlib.pyplot as plt

plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)

ONS_LAD = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
           "Local_Authority_Districts_May_2024_Boundaries_UK_BFC/FeatureServer/0/query")
GM_BOROUGHS = ["Bolton", "Bury", "Manchester", "Oldham", "Rochdale", "Salford", "Stockport", "Tameside",
               "Trafford", "Wigan"]
NAPTAN = "https://naptan.api.dft.gov.uk/v1/access-nodes"   # ATCO area 180 = Greater Manchester
STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
TOWN_HALLS = {"Manchester": (-2.24474, 53.47938), "Wigan": (-2.63117, 53.54493)}   # lon, lat (WGS 84)


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


section("versions")
print(f"python {sys.version.split()[0]} | pyproj {pyproj.__version__} (PROJ {pyproj.proj_version_str}) | "
      f"geopandas {version('geopandas')} | rasterio {rasterio.__version__} (GDAL {rasterio.__gdal_version__}) | "
      f"rioxarray {version('rioxarray')}")

# ---------------------------------------------------------------- Step 1: what is in a CRS?
section("Step 1: inspect CRSs by EPSG code")
for code in [4326, 27700, 3857, 32630]:
    crs = CRS.from_epsg(code)
    axes = ", ".join(f"{a.abbrev} ({a.direction}, {a.unit_name})" for a in crs.axis_info)
    kind = "projected" if crs.is_projected else "geographic"
    print(f"EPSG:{code:<5} {crs.name:28s} {kind:10s} datum: {crs.datum.name}")
    print(f"           axes: {axes}")
    print(f"           area of use: {crs.area_of_use.name[:70]}")

# ---------------------------------------------------------------- Step 2: transform a point, axis order
section("Step 2: transform one point (Manchester Town Hall) and the axis-order trap")
lon, lat = TOWN_HALLS["Manchester"]
t_strict = Transformer.from_crs("EPSG:4326", "EPSG:27700")              # authority axis order: lat, lon
t_xy = Transformer.from_crs("EPSG:4326", "EPSG:27700", always_xy=True)  # always lon, lat (x, y)
print(f"strict, (lat, lon) in : E {t_strict.transform(lat, lon)[0]:,.1f}  N {t_strict.transform(lat, lon)[1]:,.1f}")
e_bad, n_bad = t_strict.transform(lon, lat)
print(f"strict, (lon, lat) in : E {e_bad:,.1f}  N {n_bad:,.1f}   <- swapped axes, silently wrong")
e, n = t_xy.transform(lon, lat)
print(f"always_xy, (lon, lat) : E {e:,.1f}  N {n:,.1f}")
back = Transformer.from_crs("EPSG:27700", "EPSG:4326", always_xy=True).transform(e, n)
print(f"round trip back to WGS 84: lon {back[0]:.6f}, lat {back[1]:.6f}")
print(f"operation used: {t_xy.get_last_used_operation().description}")

# ---------------------------------------------------------------- Step 3: which transformation? test on NaPTAN
section("Step 3: datum transformations, tested on NaPTAN bus stops")
stops_csv = DATA / "naptan_gm.csv"
if not stops_csv.exists():
    r = requests.get(NAPTAN, params={"dataFormat": "csv", "atcoAreaCodes": "180"}, timeout=300)
    r.raise_for_status()
    stops_csv.write_bytes(r.content)
naptan = pd.read_csv(stops_csv, low_memory=False)
stops = naptan[(naptan.Status == "active") & (naptan.StopType == "BCT")
               & naptan[["Easting", "Northing", "Longitude", "Latitude"]].notna().all(axis=1)].copy()
print(f"NaPTAN Greater Manchester: {len(naptan):,} access nodes, {len(stops):,} active on-street bus stops "
      f"with both lon/lat and easting/northing")

tg = TransformerGroup("EPSG:4326", "EPSG:27700", always_xy=True)
if not tg.best_available:   # the OSTN15 grid is not installed yet: fetch it once from cdn.proj.org
    pyproj.network.set_network_enabled(True)
    tg.download_grids(verbose=True)
    tg = TransformerGroup("EPSG:4326", "EPSG:27700", always_xy=True)
print(f"{len(tg.transformers)} candidate operations EPSG:4326 -> EPSG:27700, in PROJ's order (first = default):")
for tr in tg.transformers:
    print(f"  accuracy {tr.accuracy:5.1f} m  {tr.description.split(' + ')[0]}")

# an old-style .prj/proj-string with no datum shift: PROJ can only apply a 'ballpark' transformation
bng_no_datum = CRS.from_proj4("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 "
                              "+ellps=airy +units=m +no_defs")
candidates = {
    "OSTN15 grid (EPSG default)": tg.transformers[0],
    "7-parameter Helmert": next(t for t in tg.transformers if "WGS 84 (6)" in t.description),
    "proj string, no datum shift": Transformer.from_crs("EPSG:4326", bng_no_datum, always_xy=True),
}
lon_a, lat_a = stops.Longitude.to_numpy(), stops.Latitude.to_numpy()
pub = np.c_[stops.Easting.to_numpy(), stops.Northing.to_numpy()]
resid = {}
for name, tr in candidates.items():
    x, y = tr.transform(lon_a, lat_a)
    d = np.hypot(x - pub[:, 0], y - pub[:, 1])
    resid[name] = d
    print(f"{name:28s} vs NaPTAN easting/northing: median {np.median(d):6.2f} m, 95th pct "
          f"{np.percentile(d, 95):6.2f} m, max {d.max():7.2f} m")
rng = np.random.default_rng(0)   # NaPTAN eastings/northings are whole metres: what does rounding alone give?
r_only = np.hypot(*rng.uniform(-0.5, 0.5, (2, 100_000)))
print(f"rounding to whole metres alone: median {np.median(r_only):.2f} m, max {r_only.max():.2f} m")

# ---------------------------------------------------------------- Step 4: reproject a GeoDataFrame
section("Step 4: reproject a GeoDataFrame (ONS boroughs, GeoJSON in WGS 84 -> EPSG:27700)")
lad_file = DATA / "gm_lad24_wgs84.geojson"
if not lad_file.exists():
    where = "LAD24NM IN (" + ",".join(f"'{b}'" for b in GM_BOROUGHS) + ")"
    r = requests.get(ONS_LAD, params={"where": where, "outFields": "LAD24CD,LAD24NM,Shape__Area",
                                      "outSR": 4326, "f": "geojson"}, timeout=180)
    r.raise_for_status()
    lad_file.write_bytes(r.content)
boroughs = gpd.read_file(lad_file).sort_values("LAD24NM").reset_index(drop=True)
print(f"read {len(boroughs)} boroughs, CRS {boroughs.crs.to_string()}, bounds "
      + ", ".join(f"{v:.3f}" for v in boroughs.total_bounds))
boroughs_bng = boroughs.to_crs(27700)
print(f"to_crs(27700): CRS {boroughs_bng.crs.to_string()}, bounds "
      + ", ".join(f"{v:,.0f}" for v in boroughs_bng.total_bounds))
wrong = boroughs.set_crs(27700, allow_override=True)   # labels the degrees as metres; moves nothing
wb = wrong.total_bounds
wlon, wlat = wrong.to_crs(4326).union_all().centroid.coords[0]
print(f"set_crs(27700) by mistake: bounds " + ", ".join(f"{v:.3f}" for v in wb)
      + f" -> a {wb[2] - wb[0]:.2f} x {wb[3] - wb[1]:.2f} m speck at lon {wlon:.2f}, lat {wlat:.2f} (the BNG false origin)")
gm = boroughs_bng.dissolve()[["geometry"]]
boroughs_bng.to_file(OUT / "greater_manchester_boroughs_27700.gpkg", layer="boroughs", driver="GPKG")

# ---------------------------------------------------------------- Step 5: areas and distances
section("Step 5: area and distance errors in the wrong CRS")
geod = Geod(ellps="WGS84")
true_area = np.array([abs(geod.geometry_area_perimeter(g)[0]) for g in boroughs.geometry])   # geodesic, m2
areas = {}
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    areas["EPSG:4326 (degrees)"] = boroughs.area.to_numpy()
print("warning from .area in EPSG:4326:", str(caught[0].message).split(".")[0] if caught else "none")
for code in [3857, 32630, 27700]:
    areas[f"EPSG:{code}"] = boroughs.to_crs(code).area.to_numpy()
print(f"Greater Manchester, geodesic area on the WGS 84 ellipsoid: {true_area.sum() / 1e6:,.2f} km2 "
      f"(ONS Shape__Area: {boroughs.Shape__Area.sum() / 1e6:,.2f} km2)")
area_err = {}
for k, a in areas.items():
    if "degrees" in k:
        print(f"  {k:22s} total = {a.sum():.4f} square degrees (not an area)")
        continue
    err = (a.sum() / true_area.sum() - 1) * 100
    area_err[k] = err
    print(f"  {k:22s} total = {a.sum() / 1e6:9,.2f} km2  error {err:+7.3f}%")

geod_d = geod.inv(*TOWN_HALLS["Manchester"], *TOWN_HALLS["Wigan"])[2]
line = gpd.GeoSeries([LineString([TOWN_HALLS["Manchester"], TOWN_HALLS["Wigan"]])], crs=4326)
dist = {"EPSG:4326 (degrees)": line.length.iloc[0]}
for code in [3857, 32630, 27700]:
    dist[f"EPSG:{code}"] = line.to_crs(code).length.iloc[0]
print(f"Manchester Town Hall -> Wigan Town Hall, geodesic: {geod_d:,.1f} m")
dist_err = {}
for k, d in dist.items():
    if "degrees" in k:
        print(f"  {k:22s} length = {d:.5f} degrees (not a distance)")
        continue
    dist_err[k] = (d / geod_d - 1) * 100
    print(f"  {k:22s} length = {d:,.1f} m  error {dist_err[k]:+7.3f}%")

# a "1 km" buffer drawn in degrees: 0.009 degrees is about 1 km of latitude, but not of longitude
p = gpd.GeoSeries([Point(TOWN_HALLS["Manchester"])], crs=4326)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    buf_deg = p.buffer(0.009).to_crs(27700)
buf_m = p.to_crs(27700).buffer(1000)
bx0, by0, bx1, by1 = buf_deg.total_bounds
print(f"buffer(0.009) in degrees, measured in BNG: {(bx1 - bx0) / 1000:.2f} km wide x {(by1 - by0) / 1000:.2f} km tall, "
      f"area {buf_deg.area.iloc[0] / 1e6:.2f} km2 (a true 1 km buffer: {buf_m.area.iloc[0] / 1e6:.2f} km2)")
scale = CRS.from_epsg(27700)
fac = pyproj.Proj(scale).get_factors(*TOWN_HALLS["Manchester"])
print(f"BNG point scale factor at Manchester Town Hall: {fac.meridional_scale:.6f}; Web Mercator: "
      f"{pyproj.Proj(CRS.from_epsg(3857)).get_factors(*TOWN_HALLS['Manchester']).meridional_scale:.4f}")

# ---------------------------------------------------------------- Step 6: reproject a raster
section("Step 6: reproject a raster (Copernicus DEM GLO-30, EPSG:4326 -> EPSG:27700)")
bbox_ll = tuple(boroughs.total_bounds + np.array([-0.02, -0.02, 0.02, 0.02]))
dem_ll_tif = DATA / "cop_dem_glo30_gm_4326.tif"
if not dem_ll_tif.exists():
    catalog = pystac_client.Client.open(STAC_URL, modifier=pc.sign_inplace)
    items = catalog.search(collections=["cop-dem-glo-30"], bbox=bbox_ll).item_collection()
    print("tiles:", ", ".join(sorted(it.id for it in items)))
    arr, transform = rasterio.merge.merge([it.assets["data"].href for it in items], bounds=bbox_ll)
    with rasterio.open(items[0].assets["data"].href) as src:
        prof = src.profile
    prof.update(driver="GTiff", width=arr.shape[2], height=arr.shape[1], transform=transform, compress="deflate",
                tiled=True, blockxsize=512, blockysize=512)
    with rasterio.open(dem_ll_tif, "w", **prof) as dst:
        dst.write(arr)

with rasterio.open(dem_ll_tif) as src:
    dem_ll = src.read(1)
    dx, dy = src.res
    print(f"source: {src.crs.to_string()}, {src.width} x {src.height} px, pixel {dx * 3600:.1f}\" x {dy * 3600:.1f}\" "
          f"({dx:.6f} x {dy:.6f} deg), nodata {src.nodata}")
    mid_lat = (src.bounds.bottom + src.bounds.top) / 2
    w_m = geod.inv(-2.0, mid_lat, -2.0 + dx, mid_lat)[2]
    h_m = geod.inv(-2.0, mid_lat, -2.0, mid_lat + dy)[2]
    print(f"on the ground at {mid_lat:.2f} N one pixel is {w_m:.1f} m wide x {h_m:.1f} m tall")

    dst_crs = "EPSG:27700"
    dst_transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height,
                                                               *src.bounds, resolution=30)
    print(f"target grid: {dst_crs}, {width} x {height} px at 30 m, origin "
          f"E {dst_transform.c:,.0f} N {dst_transform.f:,.0f}")
    warped = {}
    for rs in [Resampling.nearest, Resampling.bilinear, Resampling.cubic]:
        dst = np.full((height, width), np.nan, "float32")
        reproject(rasterio.band(src, 1), dst, dst_transform=dst_transform, dst_crs=dst_crs,
                  resampling=rs, dst_nodata=np.nan)
        warped[rs.name] = dst
    prof = src.profile | {"crs": dst_crs, "transform": dst_transform, "width": width, "height": height,
                          "nodata": np.nan, "dtype": "float32"}
with rasterio.open(OUT / "cop_dem_gm_27700.tif", "w", **prof) as dst:
    dst.write(warped["bilinear"], 1)

valid = np.isfinite(warped["bilinear"]) & np.isfinite(warped["nearest"])
for name in ["nearest", "cubic"]:
    d = (warped[name] - warped["bilinear"])[valid & np.isfinite(warped[name])]
    print(f"{name:8s} minus bilinear: mean {d.mean():+.3f} m, RMS {np.sqrt((d ** 2).mean()):.2f} m, "
          f"99th pct |diff| {np.percentile(np.abs(d), 99):.2f} m")

# the same in one line with rioxarray, then clip to the county boundary
da = rioxarray.open_rasterio(dem_ll_tif, masked=True).squeeze("band", drop=True)
da_bng = da.rio.reproject("EPSG:27700", resolution=30, resampling=Resampling.bilinear)
print(f"rioxarray: {da_bng.rio.crs}, shape {da_bng.shape}, res {da_bng.rio.resolution()}")
same_grid = da_bng.rio.transform() == dst_transform and da_bng.shape == (height, width)
diff = np.abs(da_bng.values - warped["bilinear"]) if same_grid else None
print(f"same grid as rasterio.warp: {same_grid}"
      + (f", max |difference| {np.nanmax(diff):.3f} m" if same_grid else ""))
dem_gm = da_bng.rio.clip(gm.geometry, gm.crs)
iy, ix = np.unravel_index(np.nanargmax(dem_gm.values), dem_gm.shape)
hx, hy = float(dem_gm.x[ix]), float(dem_gm.y[iy])
hlon, hlat = Transformer.from_crs(27700, 4326, always_xy=True).transform(hx, hy)
print(f"clipped to Greater Manchester: {int(np.isfinite(dem_gm.values).sum()):,} cells; elevation "
      f"{float(dem_gm.min()):.0f}-{float(dem_gm.max()):.0f} m, mean {float(dem_gm.mean()):.0f} m")
print(f"highest cell: {float(dem_gm.max()):.0f} m at E {hx:,.0f} N {hy:,.0f} (lon {hlon:.4f}, lat {hlat:.4f})")

# ---------------------------------------------------------------- figures
section("figures")
fig, axes = plt.subplots(1, 3, figsize=(13, 4.6), layout="constrained")
for ax, (crs_lbl, g) in zip(axes, [("EPSG:4326 (WGS 84, degrees)", boroughs),
                                   ("EPSG:3857 (Web Mercator, 'metres')", boroughs.to_crs(3857)),
                                   ("EPSG:27700 (British National Grid, metres)", boroughs_bng)]):
    g.plot(ax=ax, column="LAD24NM", cmap="Spectral", edgecolor="white", linewidth=0.6)
    if "4326" in crs_lbl:
        ax.set_aspect("equal")   # plotted naively, 1 degree east = 1 degree north
    ax.set_title(crs_lbl)
    ax.tick_params(labelsize=7)
    ax.ticklabel_format(useOffset=False, style="plain")
    for lbl in ax.get_xticklabels():
        lbl.set_rotation(30)
for _, row in boroughs_bng.iterrows():
    c = row.geometry.representative_point()
    axes[2].annotate(row.LAD24NM, (c.x, c.y), ha="center", fontsize=6.5, color="#111111")
fig.suptitle("The same ten boroughs in three CRSs (axes drawn with equal units)")
fig.savefig(OUT / "01-greater-manchester-three-crs.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(8, 3.8), layout="constrained")
bins = np.logspace(-3, 2.5, 80)
colours = {"OSTN15 grid (EPSG default)": "#3288bd", "7-parameter Helmert": "#f46d43",
           "proj string, no datum shift": "#9e0142"}
for name, d in resid.items():
    ax.hist(np.clip(d, 1e-3, None), bins=bins, color=colours[name], alpha=0.8,
            label=f"{name} (median {np.median(d):.2f} m)")
ax.set_xscale("log")
ax.set_xlabel("distance between pyproj result and NaPTAN's published easting/northing (m, log scale)")
ax.set_ylabel("bus stops")
ax.set_title(f"WGS 84 -> British National Grid for {len(stops):,} Greater Manchester bus stops")
ax.legend(frameon=False)
fig.savefig(OUT / "02-transformation-accuracy.png")
plt.close(fig)

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2), layout="constrained", gridspec_kw={"width_ratios": [1.5, 1]})
ax = axes[0]
labels = list(area_err)
yy = np.arange(len(labels))
for off, data, lbl, col in [(-0.18, area_err, "area of Greater Manchester", "#3288bd"),
                            (0.18, dist_err, "Manchester-Wigan distance", "#f46d43")]:
    vals = np.array([data[k] for k in labels])
    bars = ax.barh(yy + off, np.abs(vals), height=0.34, color=col, label=lbl)
    for b, v in zip(bars, vals):
        ax.annotate(f"{v:+.3f}%", (abs(v), b.get_y() + b.get_height() / 2), xytext=(4, 0),
                    textcoords="offset points", va="center", fontsize=8)
ax.set_yticks(yy, labels)
ax.set_xscale("log")
ax.set_xlim(0.01, 3000)
ax.invert_yaxis()
ax.set_xlabel("|error| against the geodesic value on the WGS 84 ellipsoid (%, log scale)")
ax.legend(frameon=False, loc="lower right")
ax.set_title("Measuring in the wrong CRS")
ax = axes[1]
boroughs_bng[boroughs_bng.LAD24NM == "Manchester"].boundary.plot(ax=ax, color="#999999", lw=0.6)
buf_deg.plot(ax=ax, facecolor="#fdae61", edgecolor="#d53e4f", alpha=0.7, label="buffer(0.009) in degrees")
buf_m.plot(ax=ax, facecolor="none", edgecolor="#3288bd", lw=2)
mx, my = p.to_crs(27700).iloc[0].coords[0]
ax.plot(mx, my, "k+", ms=8)
ax.set_xlim(mx - 2200, mx + 2200); ax.set_ylim(my - 1600, my + 1600)
ax.set_title("'1 km' buffer in degrees (orange)\nvs a true 1 km buffer in BNG (blue)")
ax.ticklabel_format(useOffset=False, style="plain")
ax.tick_params(labelsize=7)
fig.savefig(OUT / "03-area-distance-errors.png")
plt.close(fig)

with rasterio.open(dem_ll_tif) as src:
    ext_ll = [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top]
ext_bng = [dst_transform.c, dst_transform.c + width * 30, dst_transform.f - height * 30, dst_transform.f]
fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained")
dem_plot = np.where(dem_ll == -32767, np.nan, dem_ll)
im = axes[0].imshow(dem_plot, extent=ext_ll, cmap="Spectral_r", vmin=0, vmax=550)
boroughs.boundary.plot(ax=axes[0], color="#222222", lw=0.5)
axes[0].set_title(f"Copernicus DEM as delivered: EPSG:4326, {dx * 3600:.1f}\" x {dy * 3600:.1f}\" pixels")
axes[0].set_aspect("equal")
axes[1].imshow(warped["bilinear"], extent=ext_bng, cmap="Spectral_r", vmin=0, vmax=550)
boroughs_bng.boundary.plot(ax=axes[1], color="#222222", lw=0.5)
axes[1].set_title("Reprojected to EPSG:27700 at 30 m (bilinear)")
axes[1].ticklabel_format(useOffset=False, style="plain")
axes[1].plot(hx, hy, marker="^", color="black", ms=7)
axes[1].annotate(f"{float(dem_gm.max()):.0f} m", (hx, hy), xytext=(6, -10), textcoords="offset points", fontsize=8)
fig.colorbar(im, ax=axes, shrink=0.8, label="elevation (m)")
fig.savefig(OUT / "04-dem-reprojected.png")
plt.close(fig)


# zoom on 600 m near the highest cell: the source pixel grid drawn in BNG, then the two resampled 30 m grids
x0, y1 = (hx // 30) * 30 - 450, (hy // 30) * 30 + 300
x1, y0 = x0 + 600, y1 - 600
fig, axes = plt.subplots(1, 3, figsize=(13, 4.4), layout="constrained")
with rasterio.open(dem_ll_tif) as src:
    lon0, lat0, lon1, lat1 = transform_bounds("EPSG:27700", src.crs, x0 - 60, y0 - 60, x1 + 60, y1 + 60)
    win = rasterio.windows.from_bounds(lon0, lat0, lon1, lat1, src.transform).round_offsets().round_lengths()
    z = src.read(1, window=win)
    wt = src.window_transform(win)
cols_e = wt.c + wt.a * np.arange(z.shape[1] + 1)
rows_e = wt.f + wt.e * np.arange(z.shape[0] + 1)
LON, LAT = np.meshgrid(cols_e, rows_e)
XE, YN = t_xy.transform(LON, LAT)
vmin, vmax = np.nanpercentile(warped["bilinear"][int((dst_transform.f - y1) / 30):int((dst_transform.f - y0) / 30),
                                                 int((x0 - dst_transform.c) / 30):int((x1 - dst_transform.c) / 30)], [1, 99])
axes[0].pcolormesh(XE, YN, z, cmap="Spectral_r", vmin=vmin, vmax=vmax, edgecolors="#ffffff", linewidth=0.3)
axes[0].set_title(f"Source pixels ({dx * 3600:.1f}\" x {dy * 3600:.1f}\" = {w_m:.1f} x {h_m:.1f} m)\ndrawn in BNG")
for ax, name in zip(axes[1:], ["nearest", "bilinear"]):
    im = ax.imshow(warped[name], extent=ext_bng, cmap="Spectral_r", vmin=vmin, vmax=vmax)
    rr = np.arange(x0, x1 + 1, 30)
    ax.vlines(rr, y0, y1, color="#ffffff", lw=0.3); ax.hlines(np.arange(y0, y1 + 1, 30), x0, x1, color="#ffffff", lw=0.3)
    ax.set_title(f"{name} resampling to 30 x 30 m BNG pixels")
for ax in axes:
    ax.set_xlim(x0, x1); ax.set_ylim(y0, y1); ax.set_aspect("equal")
    ax.ticklabel_format(useOffset=False, style="plain"); ax.tick_params(labelsize=7)
fig.colorbar(im, ax=axes, shrink=0.8, label="elevation (m)")
fig.suptitle("600 m of Pennine moorland near the highest cell: nearest copies whole source pixels, bilinear interpolates")
fig.savefig(OUT / "05-resampling.png")
plt.close(fig)
print("saved:", ", ".join(sorted(p.name for p in OUT.glob("*.png"))))