Source code: Raster to vector in Python: polygonize with rasterio, rasterize boundaries and convert a GeoTIFF to a Cloud Optimized GeoTIFF

"""Raster to vector in Python: polygonise ESA WorldCover, rasterize ONS boundaries and write a Cloud Optimized GeoTIFF.

Greater Manchester, UK: water and built-up land are vectorised from ESA WorldCover 2021 (10 m) with
rasterio.features.shapes, simplified and measured; the ten ONS local authority districts are rasterized onto
the WorldCover grid for per-district land-cover shares; and the subset is converted to a COG with rio-cogeo,
validated and read over HTTP with range requests.

Part 3 of the TechnicalGIS "Geospatial Python essentials" series.
TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/raster-to-vector-python/
"""
import http.server
import logging
import os
import re
import sys
import threading
import time
import warnings
from functools import partial
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.shutil
from rasterio import features
from rasterio.windows import from_bounds
from rasterio.transform import array_bounds
import shapely
from shapely.geometry import shape, box
from shapely.validation import explain_validity
import pyproj
import pystac_client
import planetary_computer as pc
from rio_cogeo.cogeo import cog_translate, cog_validate
from rio_cogeo.profiles import cog_profiles
import matplotlib
matplotlib.use("Agg")   # render to files only
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

warnings.filterwarnings("ignore", category=UserWarning)
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 = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Local_Authority_Districts_December_2024_Boundaries_UK_BFC/FeatureServer/0/query")
GM_CODES = [f"E080000{i:02d}" for i in range(1, 11)]    # the ten Greater Manchester metropolitan boroughs
STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
WORLDCOVER = {10: ("Tree cover", "#006400"), 20: ("Shrubland", "#ffbb22"), 30: ("Grassland", "#ffff4c"),
              40: ("Cropland", "#f096ff"), 50: ("Built-up", "#fa0000"), 60: ("Bare / sparse", "#b4b4b4"),
              70: ("Snow and ice", "#f0f0f0"), 80: ("Water", "#0064c8"), 90: ("Herbaceous wetland", "#0096a0"),
              95: ("Mangroves", "#00cf75"), 100: ("Moss and lichen", "#fae6a0")}
BNG = "EPSG:27700"


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


def mb(n):
    return n / 1024 ** 2


def nverts(geoms):
    return int(shapely.get_num_coordinates(np.asarray(geoms)).sum())


section("versions")
print(f"python {sys.version.split()[0]} | rasterio {rasterio.__version__} (GDAL {rasterio.__gdal_version__}) | "
      f"rio-cogeo {version('rio-cogeo')} | geopandas {gpd.__version__} | shapely {shapely.__version__}")

# ------------------------------------------------------------------ Step 1: boundaries
section("Step 1: Greater Manchester districts (ONS LAD December 2024, BFC)")
gm_path = DATA / "gm_districts.gpkg"
if not gm_path.exists():
    where = "LAD24CD IN (" + ",".join(f"'{c}'" for c in GM_CODES) + ")"
    r = requests.get(ONS, params={"where": where, "outFields": "LAD24CD,LAD24NM", "outSR": 27700,
                                  "f": "geojson"}, timeout=120)
    r.raise_for_status()
    gpd.GeoDataFrame.from_features(r.json()["features"], crs=BNG).to_file(gm_path)
districts = gpd.read_file(gm_path).sort_values("LAD24CD").reset_index(drop=True)
districts["zone"] = np.arange(1, len(districts) + 1, dtype="uint8")       # 1..10, 0 = outside
print(f"{len(districts)} districts, CRS {districts.crs.to_string()}, "
      f"total area {districts.area.sum() / 1e6:,.1f} km2")
gm_wgs84 = districts.to_crs(4326)
bounds = gm_wgs84.total_bounds
print("bounds (lon/lat):", np.round(bounds, 4))

# ------------------------------------------------------------------ Step 2: WorldCover subset
section("Step 2: read ESA WorldCover 2021 for Greater Manchester (Planetary Computer)")
catalog = pystac_client.Client.open(STAC_URL, modifier=pc.sign_inplace)
items = [it for it in catalog.search(collections=["esa-worldcover"], bbox=list(bounds)).items() if "2021" in it.id]
item = items[0]
href = item.assets["map"].href
print("STAC item:", item.id)

# Log GDAL's HTTP range requests so we can count the bytes actually downloaded (URLs are not printed).
ranges = []


class RangeLog(logging.Handler):
    def emit(self, record):
        m = re.search(r"Downloading (\d+)-(\d+)", record.getMessage())
        if m:
            ranges.append((int(m.group(1)), int(m.group(2))))


logging.getLogger("rasterio").addHandler(RangeLog())
logging.getLogger("rasterio").setLevel(logging.DEBUG)

plain_tif = DATA / "worldcover_gm_plain.tif"
with rasterio.Env(CPL_DEBUG=True, GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
    with rasterio.open(href) as src:
        remote_profile = src.profile
        win = from_bounds(*bounds, transform=src.transform).round_offsets().round_lengths()
        lc = src.read(1, window=win)
        transform = src.window_transform(win)
        crs = src.crs
logging.getLogger("rasterio").setLevel(logging.WARNING)
remote_size = int(requests.head(href, timeout=60).headers["Content-Length"])
fetched = sum(b - a + 1 for a, b in ranges)
print(f"remote COG: {remote_profile['width']} x {remote_profile['height']} px, "
      f"{remote_profile['blockxsize']} x {remote_profile['blockysize']} tiles, {remote_profile['compress']}, "
      f"{mb(remote_size):,.1f} MB on the server")
print(f"window read: {lc.shape[1]} x {lc.shape[0]} px using {len(ranges)} HTTP range requests, "
      f"{mb(fetched):.2f} MB downloaded ({100 * fetched / remote_size:.1f}% of the file)")

# Save a plain (non-cloud-optimised) GeoTIFF: GDAL's defaults = uncompressed strips, no overviews.
profile = {"driver": "GTiff", "dtype": "uint8", "count": 1, "width": lc.shape[1], "height": lc.shape[0],
           "crs": crs, "transform": transform, "nodata": 0}
with rasterio.open(plain_tif, "w", **profile) as dst:
    dst.write(lc, 1)
print(f"saved {plain_tif.name}: {mb(plain_tif.stat().st_size):.1f} MB, CRS {crs.to_string()}, "
      f"pixel {transform.a:.3e} x {-transform.e:.3e} degrees")

# A "10 m" pixel in EPSG:4326 is not 100 m2: its area depends on latitude.
geod = pyproj.Geod(ellps="WGS84")
lat_top = transform.f + np.arange(lc.shape[0]) * transform.e
row_area = np.array([abs(geod.polygon_area_perimeter([0, transform.a, transform.a, 0],
                                                      [t, t, t + transform.e, t + transform.e])[0]) for t in lat_top])
print(f"pixel size on the ground: {row_area.min():.1f} to {row_area.max():.1f} m2 "
      f"(north to south edge of the subset), not 100 m2")

# ------------------------------------------------------------------ Step 3: rasterize the districts
section("Step 3: rasterize the districts onto the WorldCover grid")
zones = features.rasterize(zip(gm_wgs84.geometry, districts["zone"]), out_shape=lc.shape, transform=transform,
                           fill=0, dtype="uint8")
zones_touch = features.rasterize(((g, 1) for g in gm_wgs84.geometry), out_shape=lc.shape, transform=transform,
                                 fill=0, dtype="uint8", all_touched=True)
inside = zones > 0
area_px = (inside * row_area[:, None]).sum()
print(f"pixels inside Greater Manchester: {inside.sum():,} (centre rule) vs {int(zones_touch.sum()):,} (all_touched)")
print(f"rasterized area {area_px / 1e6:,.1f} km2 vs vector area {districts.area.sum() / 1e6:,.1f} km2 "
      f"({100 * (area_px / districts.area.sum() - 1):+.2f}%)")

# Cross-tabulate land cover by district: one bincount over (zone, class) pairs, weighted by pixel area.
weights = np.broadcast_to(row_area[:, None], lc.shape)
key = zones.astype(np.int32) * 256 + lc
tab = np.bincount(key[inside], weights=weights[inside], minlength=11 * 256).reshape(11, 256)[1:]
classes = [c for c in WORLDCOVER if tab[:, c].sum() > 0]
share = pd.DataFrame(tab[:, classes] / tab.sum(axis=1, keepdims=True) * 100,
                     index=districts["LAD24NM"], columns=[WORLDCOVER[c][0] for c in classes])
gm_share = pd.Series(tab[:, classes].sum(axis=0) / tab.sum() * 100, index=share.columns)
print("Greater Manchester land cover (% of area):")
print(gm_share.round(1).to_string())
print("\nby district (% of district area):")
print(share[["Built-up", "Tree cover", "Grassland", "Cropland", "Water"]].round(1)
      .sort_values("Built-up", ascending=False).to_string())

# ------------------------------------------------------------------ Step 4: vectorise water
section("Step 4: vectorise water with rasterio.features.shapes")
water = ((lc == 80) & inside).astype("uint8")
results = {}
for conn in (4, 8):
    t0 = time.perf_counter()
    polys = [shape(g) for g, v in features.shapes(water, mask=water.astype(bool), transform=transform,
                                                  connectivity=conn)]
    results[conn] = (polys, time.perf_counter() - t0)
    print(f"connectivity={conn}: {len(polys):,} polygons in {results[conn][1]:.2f} s")
water_gdf = gpd.GeoDataFrame({"class": ["Water"] * len(results[4][0])}, geometry=results[4][0], crs=crs).to_crs(BNG)
water_gdf["area_m2"] = water_gdf.area
water_px_area = (water.astype(bool) * row_area[:, None]).sum()
print(f"total water area: polygons {water_gdf.area_m2.sum() / 1e6:.3f} km2, "
      f"pixel count x pixel area {water_px_area / 1e6:.3f} km2")
one_px = water_gdf.area_m2 < 1.5 * row_area.mean()
print(f"single-pixel polygons: {one_px.sum():,} of {len(water_gdf):,} ({100 * one_px.mean():.0f}%), "
      f"holding {100 * water_gdf.area_m2[one_px].sum() / water_gdf.area_m2.sum():.1f}% of the water area")
q = water_gdf.area_m2.quantile([0.5, 0.9, 0.99])
print(f"polygon area: median {q[0.5]:.0f} m2, 90th percentile {q[0.9]:.0f} m2, "
      f"99th {q[0.99]:,.0f} m2, max {water_gdf.area_m2.max() / 1e4:,.1f} ha")

# Round trip: burn the polygons back and compare with the source mask.
back = features.rasterize(((g, 1) for g in water_gdf.to_crs(crs).geometry), out_shape=lc.shape,
                          transform=transform, fill=0, dtype="uint8")
print(f"round trip (shapes -> to_crs -> back -> rasterize): {int((back != water).sum())} pixels differ "
      f"out of {water.size:,}")

# Clean up with a sieve before vectorising: remove water patches smaller than 10 pixels.
sieved = features.sieve(water, size=10, connectivity=4)
n_sieved = sum(1 for _ in features.shapes(sieved, mask=sieved.astype(bool), transform=transform))
sieved_area = (sieved.astype(bool) * row_area[:, None]).sum()
print(f"after sieve(size=10): {n_sieved:,} polygons, water area {sieved_area / 1e6:.3f} km2 "
      f"({100 * (sieved_area / water_px_area - 1):+.1f}%)")

big = water_gdf[water_gdf.area_m2 >= 50_000].copy()
big["district"] = [districts.LAD24NM[districts.intersects(g)].str.cat(sep=" / ") for g in big.geometry]
big = big.sort_values("area_m2", ascending=False)
cent = big.geometry.centroid.to_crs(4326)
print(f"{len(big)} water bodies >= 5 ha; the largest:")
for (_, rw), pt in list(zip(big.iterrows(), cent))[:8]:
    print(f"  {rw.area_m2 / 1e4:6.1f} ha  {rw.district:<20} centroid {pt.y:.4f} N, {pt.x:.4f} E")

# ------------------------------------------------------------------ Step 5: built-up and simplification
section("Step 5: vectorise built-up land and simplify")
built_raw = ((lc == 50) & inside).astype("uint8")
for conn in (8, 4):
    built = features.sieve(built_raw, size=10, connectivity=conn)
    t0 = time.perf_counter()
    polys = [shape(g) for g, v in features.shapes(built, mask=built.astype(bool), transform=transform,
                                                  connectivity=conn)]
    t_built = time.perf_counter() - t0
    bad = [p for p in polys if not p.is_valid]
    print(f"connectivity={conn}: {len(polys):,} polygons in {t_built:.1f} s, {len(bad):,} invalid"
          + (f" (e.g. {explain_validity(bad[0]).split('[')[0]})" if bad else ""))
    if bad:
        fixed = shapely.make_valid(np.array(bad))
        print(f"  shapely.make_valid() repairs them; geometry types now: "
              f"{sorted(set(shapely.get_type_id(fixed).tolist()))} (6 = MultiPolygon)")
built_gdf = gpd.GeoDataFrame({"class": ["Built-up"] * len(polys)}, geometry=polys, crs=crs).to_crs(BNG)  # 4-connected
built_gdf["area_m2"] = built_gdf.area
largest = built_gdf.area_m2.max()
print(f"{len(built_gdf):,} built-up polygons kept (sieve size=10, connectivity=4); "
      f"largest {largest / 1e6:.1f} km2 ({100 * largest / built_gdf.area_m2.sum():.0f}% of all built-up area); "
      f"{sum(len(p.interiors) for p in built_gdf.geometry):,} holes in total")

rows_s = []
for tol in [0, 2.5, 5, 10, 20, 40]:
    g = built_gdf.geometry if tol == 0 else built_gdf.geometry.simplify(tol, preserve_topology=True)
    path = DATA / f"built_simplify_{tol}.gpkg"
    path.unlink(missing_ok=True)        # overwriting a layer does not shrink an existing GeoPackage
    gpd.GeoDataFrame(geometry=g, crs=BNG).to_file(path)
    rows_s.append({"tolerance_m": tol, "vertices": nverts(g.values), "area_km2": g.area.sum() / 1e6,
                   "invalid": int((~g.is_valid).sum()), "gpkg_MB": mb(path.stat().st_size)})
simp = pd.DataFrame(rows_s)
simp["area_change_%"] = 100 * (simp.area_km2 / simp.area_km2.iloc[0] - 1)
simp["vertices_%"] = 100 * simp.vertices / simp.vertices.iloc[0]
print(simp.round(3).to_string(index=False))
for f in ("gm_built_up_simplified.gpkg", "gm_water.gpkg"):
    (OUT / f).unlink(missing_ok=True)
built_gdf.assign(geometry=built_gdf.geometry.simplify(10)).to_file(OUT / "gm_built_up_simplified.gpkg")
water_gdf.to_file(OUT / "gm_water.gpkg")

# A full land-cover coverage for a 4 km x 4 km tile (Salford Quays / Trafford Park), all classes at once.
section("Step 5b: simplify a polygon coverage without gaps (shapely.coverage_simplify)")
tile_ll = gpd.GeoSeries([box(377_000, 395_000, 381_000, 399_000)], crs=BNG).to_crs(crs).total_bounds
tw = from_bounds(*tile_ll, transform=transform).round_offsets().round_lengths()
sub = lc[tw.row_off:tw.row_off + tw.height, tw.col_off:tw.col_off + tw.width]
sub_t = rasterio.windows.transform(tw, transform)
cov = gpd.GeoDataFrame([{"code": int(v), "geometry": shape(g)} for g, v in features.shapes(sub, transform=sub_t)],
                       crs=crs).to_crs(BNG)
cov_union = shapely.union_all(cov.geometry.values).area
print(f"{len(cov):,} polygons covering {cov_union / 1e6:.3f} km2, {nverts(cov.geometry.values):,} vertices")
for name, g in [("simplify(10)", cov.geometry.simplify(10).values),
                ("coverage_simplify(10)", shapely.coverage_simplify(cov.geometry.values, 10))]:
    union = shapely.union_all(g).area
    overlap = shapely.area(g).sum() - union
    print(f"{name:<22} vertices {nverts(g):>7,}  overlaps {overlap:>8,.0f} m2  gaps {max(cov_union - union, 0):>8,.0f} m2")

# ------------------------------------------------------------------ Step 6: Cloud Optimized GeoTIFF
section("Step 6: convert the GeoTIFF to a Cloud Optimized GeoTIFF")
masked_tif = DATA / "worldcover_gm_masked.tif"
with rasterio.open(masked_tif, "w", **profile) as dst:
    dst.write(np.where(inside, lc, 0).astype("uint8"), 1)

ok, errors, warns = cog_validate(str(masked_tif), quiet=True)
print(f"plain GeoTIFF valid COG? {ok}\n  errors: {errors}\n  warnings: {warns}")

cog_tif = OUT / "worldcover_gm_cog.tif"
t0 = time.perf_counter()
cog_translate(masked_tif, cog_tif, cog_profiles.get("deflate"), overview_resampling="mode", quiet=True)
t_cog = time.perf_counter() - t0
ok, errors, warns = cog_validate(str(cog_tif), quiet=True)
print(f"rio-cogeo COG written in {t_cog:.1f} s; valid COG? {ok}; errors: {errors}; warnings: {warns}")

gdal_cog = OUT / "worldcover_gm_gdal_cog.tif"
t0 = time.perf_counter()
rasterio.shutil.copy(masked_tif, gdal_cog, driver="COG", COMPRESS="DEFLATE", OVERVIEW_RESAMPLING="MODE",
                     BLOCKSIZE=512)
t_gdal = time.perf_counter() - t0
ok2, errors2, _ = cog_validate(str(gdal_cog), quiet=True)
print(f"GDAL COG driver (rasterio.shutil.copy) in {t_gdal:.1f} s: valid COG? {ok2}; errors: {errors2}")

zstd_cog = OUT / "worldcover_gm_cog_zstd.tif"
cog_translate(masked_tif, zstd_cog, cog_profiles.get("zstd"), overview_resampling="mode", quiet=True)
lzw_cog = OUT / "worldcover_gm_cog_lzw.tif"
cog_translate(masked_tif, lzw_cog, cog_profiles.get("lzw"), overview_resampling="mode", quiet=True)
deflate_strips = DATA / "worldcover_gm_deflate_strips.tif"     # compressed but not tiled, no overviews
with rasterio.open(masked_tif) as m, rasterio.open(deflate_strips, "w", **profile, compress="deflate") as dst:
    dst.write(m.read(1), 1)

sizes = {"Plain GeoTIFF (uncompressed)": masked_tif, "GeoTIFF, DEFLATE, strips": deflate_strips, "COG, LZW": lzw_cog, "COG, DEFLATE (rio-cogeo)": cog_tif,
         "COG, DEFLATE (GDAL driver)": gdal_cog, "COG, ZSTD": zstd_cog}
size_tab = pd.Series({k: mb(p.stat().st_size) for k, p in sizes.items()})
print("file sizes (MB):")
print(size_tab.round(2).to_string())
with rasterio.open(cog_tif) as c, rasterio.open(masked_tif) as m:
    print(f"COG pixels identical to the source: {np.array_equal(c.read(1), m.read(1))}; "
          f"overview factors {c.overviews(1)}; blocks {c.block_shapes[0]}; source blocks {m.block_shapes[0]}")

# ------------------------------------------------------------------ Step 7: HTTP range reads
section("Step 7: read a window and a thumbnail over HTTP (range requests)")
served = []


class RangeHandler(http.server.SimpleHTTPRequestHandler):
    """Static file server that honours single HTTP Range requests and records the bytes sent."""

    def log_message(self, *args):
        pass

    def send_head(self):
        path = self.translate_path(self.path)
        size = os.path.getsize(path)
        rng = re.match(r"bytes=(\d+)-(\d*)", self.headers.get("Range", ""))
        f = open(path, "rb")
        start, end = 0, size - 1
        if rng:
            start = int(rng.group(1))
            end = min(int(rng.group(2)) if rng.group(2) else size - 1, size - 1)
            self.send_response(206)
            self.send_header("Content-Range", f"bytes {start}-{end}/{size}")
            f.seek(start)
        else:
            self.send_response(200)
        self.send_header("Content-Type", "image/tiff")
        self.send_header("Accept-Ranges", "bytes")
        self.send_header("Content-Length", str(end - start + 1))
        self.end_headers()
        self._remaining = end - start + 1
        if self.command == "GET":
            served.append(self._remaining)
        return f

    def copyfile(self, source, outputfile):
        outputfile.write(source.read(self._remaining))


server = http.server.ThreadingHTTPServer(("127.0.0.1", 0), partial(RangeHandler, directory=str(Path.cwd())))
threading.Thread(target=server.serve_forever, daemon=True).start()
base = f"http://127.0.0.1:{server.server_address[1]}"

window_ll = gpd.GeoSeries([box(382_000, 396_500, 386_000, 400_500)], crs=BNG).to_crs(crs).total_bounds  # city centre
http_rows = []
for label, path in [("Plain GeoTIFF", masked_tif), ("COG (DEFLATE)", cog_tif)]:
    for job in ("4 km x 4 km window, full resolution", "whole-area thumbnail, 1/16 resolution"):
        served.clear()
        url = f"/vsicurl/{base}/{path.as_posix()}?run={len(http_rows)}"   # unique URL so GDAL's cache is cold
        t0 = time.perf_counter()
        with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif"):
            with rasterio.open(url) as r:
                if job.startswith("4 km"):
                    arr = r.read(1, window=from_bounds(*window_ll, transform=r.transform))
                else:
                    arr = r.read(1, out_shape=(r.height // 16, r.width // 16))
        http_rows.append({"file": label, "job": job, "requests": len(served), "MB": mb(sum(served)),
                          "pct_of_file": 100 * sum(served) / path.stat().st_size,
                          "pixels": f"{arr.shape[1]} x {arr.shape[0]}", "seconds": time.perf_counter() - t0})
server.shutdown()
http_tab = pd.DataFrame(http_rows)
print(http_tab.round(3).to_string(index=False))

# ------------------------------------------------------------------ Figures
section("figures")
codes = sorted(WORLDCOVER)
idx_lut = np.zeros(256, dtype="uint8")       # map non-contiguous class codes to colour indices
for i, c in enumerate(codes):
    idx_lut[c] = i
cmap_i = ListedColormap([WORLDCOVER[c][1] for c in codes])
b = array_bounds(*lc.shape, transform)
extent = [b[0], b[2], b[1], b[3]]
present = [c for c in codes if c in classes and gm_share[WORLDCOVER[c][0]] >= 0.05]


def show_lc(ax, arr, ext, mask=None):
    a = np.ma.masked_where((arr == 0) if mask is None else ~mask, idx_lut[arr])
    ax.imshow(a, cmap=cmap_i, vmin=-0.5, vmax=len(codes) - 0.5, extent=ext, interpolation="nearest")


# 01 overview map
fig, ax = plt.subplots(figsize=(9, 5.4))
show_lc(ax, lc, extent, mask=inside)
gm_wgs84.boundary.plot(ax=ax, color="#222222", lw=0.6)
for _, rw in gm_wgs84.iterrows():
    p = rw.geometry.representative_point()
    ax.annotate(rw.LAD24NM, (p.x, p.y), ha="center", fontsize=7.5, color="#111111",
                bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="none", alpha=0.75))
ax.set_aspect(1 / np.cos(np.deg2rad(53.5)))
ax.set_xlabel("Longitude (°)"); ax.set_ylabel("Latitude (°)")
ax.legend(handles=[Patch(color=WORLDCOVER[c][1], label=f"{WORLDCOVER[c][0]} ({gm_share[WORLDCOVER[c][0]]:.1f}%)")
                   for c in present], loc="upper left", bbox_to_anchor=(1.01, 1), frameon=False, fontsize=8,
          title="ESA WorldCover 2021", title_fontsize=8.5)
ax.set_title("Greater Manchester land cover at 10 m, clipped by rasterizing the ten ONS districts", fontsize=10)
fig.savefig(OUT / "01-worldcover-greater-manchester.png", facecolor="white")
plt.close(fig)

# 02 zoom: raster vs polygons vs simplified (Salford Quays)
zoom_ll = gpd.GeoSeries([box(380_300, 396_950, 381_300, 397_650)], crs=BNG).to_crs(crs).total_bounds
zw = from_bounds(*zoom_ll, transform=transform).round_offsets().round_lengths()
zsub = lc[zw.row_off:zw.row_off + zw.height, zw.col_off:zw.col_off + zw.width]
zt = rasterio.windows.transform(zw, transform)
zb = array_bounds(*zsub.shape, zt)
zbox = gpd.GeoSeries([box(*zb)], crs=crs).to_crs(BNG).total_bounds
zpolys = gpd.GeoDataFrame([{"code": int(v), "geometry": shape(g)} for g, v in features.shapes(zsub, transform=zt)],
                          crs=crs).to_crs(BNG)
zcov10 = zpolys.assign(geometry=shapely.coverage_simplify(zpolys.geometry.values, 10))
fig, axs = plt.subplots(1, 3, figsize=(12, 3.9))
show_lc(axs[0], zsub, [zb[0], zb[2], zb[1], zb[3]])
axs[0].set_aspect(1 / np.cos(np.deg2rad(53.47)))
axs[0].set_title(f"1. Raster: {zsub.size:,} WorldCover pixels", fontsize=9)
for ax, gdf, title in [(axs[1], zpolys, f"2. features.shapes(): {len(zpolys):,} polygons,\n"
                                         f"{nverts(zpolys.geometry.values):,} vertices"),
                       (axs[2], zcov10, f"3. coverage_simplify(10 m):\n{nverts(zcov10.geometry.values):,} vertices")]:
    gdf.plot(ax=ax, color=[WORLDCOVER[c][1] for c in gdf.code], edgecolor="#1a1a1a", linewidth=0.4)
    ax.set_xlim(zbox[0], zbox[2]); ax.set_ylim(zbox[1], zbox[3])
    ax.set_title(title, fontsize=9)
for ax in axs:
    ax.set_xticks([]); ax.set_yticks([])
fig.suptitle("Salford Quays, 1 km x 0.7 km: from pixels to polygons", fontsize=10)
fig.savefig(OUT / "02-pixels-to-polygons.png", facecolor="white")
plt.close(fig)

# 03 simplification trade-off
fig, axs = plt.subplots(1, 2, figsize=(10, 3.8))
axs[0].plot(simp.tolerance_m, simp["vertices_%"], "o-", color="#3288bd", label="Vertices")
axs[0].plot(simp.tolerance_m, simp.gpkg_MB / simp.gpkg_MB.iloc[0] * 100, "s--", color="#d53e4f",
            label="GeoPackage file size")
for _, rw in simp.drop_duplicates("vertices").iterrows():
    axs[0].annotate(f"{rw.vertices / 1e3:,.0f}k", (rw.tolerance_m, rw["vertices_%"]), textcoords="offset points",
                    xytext=(4, 5), fontsize=7)
axs[0].set_xlabel("Simplification tolerance (m)"); axs[0].set_ylabel("% of the unsimplified layer")
axs[0].set_ylim(0, 108); axs[0].legend(frameon=False, fontsize=8); axs[0].grid(alpha=0.3)
axs[0].set_title(f"{len(built_gdf):,} built-up polygons: vertices and file size", fontsize=9)
axs[1].bar(simp.tolerance_m.astype(str), simp["area_change_%"], color="#fdae61", edgecolor="#5e4fa2")
axs[1].axhline(0, color="#333333", lw=0.8)
axs[1].set_xlabel("Simplification tolerance (m)"); axs[1].set_ylabel("Change in total area (%)")
axs[1].set_title("Total built-up area after simplify()", fontsize=9); axs[1].grid(alpha=0.3, axis="y")
fig.savefig(OUT / "03-simplification-tradeoff.png", facecolor="white")
plt.close(fig)

# 04 rasterized districts + land cover by district
fig, axs = plt.subplots(1, 2, figsize=(12.5, 4.8), gridspec_kw={"width_ratios": [1.15, 1], "wspace": 0.35})
zcm = ListedColormap(plt.get_cmap("Spectral")(np.linspace(0.05, 0.95, 10)))
axs[0].imshow(np.ma.masked_equal(zones, 0), cmap=zcm, extent=extent, interpolation="nearest")
axs[0].set_aspect(1 / np.cos(np.deg2rad(53.5)))
for _, rw in gm_wgs84.iterrows():
    p = rw.geometry.representative_point()
    axs[0].annotate(f"{rw.zone}", (p.x, p.y), ha="center", va="center", fontsize=8, weight="bold")
axs[0].set_title("features.rasterize(): district codes 1–10 on the 10 m grid", fontsize=9)
axs[0].set_xlabel("Longitude (°)"); axs[0].set_ylabel("Latitude (°)")
order = share.sort_values("Built-up").index
zone_of = districts.set_index("LAD24NM")["zone"]
left = np.zeros(len(order))
for c in present:
    vals = share.loc[order, WORLDCOVER[c][0]].values
    axs[1].barh([f"{zone_of[n]}. {n}" for n in order], vals, left=left, color=WORLDCOVER[c][1],
                edgecolor="white", linewidth=0.3, label=WORLDCOVER[c][0])
    left += vals
axs[1].set_xlim(0, 100); axs[1].set_xlabel("% of district area")
axs[1].legend(loc="upper center", bbox_to_anchor=(0.5, -0.13), ncol=4, frameon=False, fontsize=7.5)
axs[1].set_title("Land cover by district (one bincount over district × class)", fontsize=9)
fig.savefig(OUT / "04-rasterize-districts.png", facecolor="white")
plt.close(fig)

# 05 COG sizes and HTTP bytes
fig, axs = plt.subplots(1, 2, figsize=(11.5, 3.9), gridspec_kw={"width_ratios": [1, 1.15]})
cols = ["#9e0142", "#f46d43", "#fdae61", "#3288bd", "#66c2a5", "#5e4fa2"]
axs[0].barh(size_tab.index[::-1], size_tab.values[::-1], color=cols[::-1])
for i, v in enumerate(size_tab.values[::-1]):
    axs[0].text(v, i, f" {v:.1f} MB", va="center", fontsize=8)
axs[0].set_xlabel("File size (MB)"); axs[0].set_xlim(0, size_tab.max() * 1.3)
axs[0].set_title("The same Greater Manchester land-cover grid on disk", fontsize=9)
jobs = list(http_tab.job.unique())
x = np.arange(len(jobs))
for j, (label, colr) in enumerate([("Plain GeoTIFF", "#9e0142"), ("COG (DEFLATE)", "#3288bd")]):
    sel = http_tab[http_tab.file == label].set_index("job").loc[jobs]
    bars = axs[1].bar(x + (j - 0.5) * 0.38, sel.MB, 0.38, color=colr, label=label)
    for bar, v, n in zip(bars, sel.MB, sel.requests):
        axs[1].text(bar.get_x() + bar.get_width() / 2, v * 1.2, f"{v:.2f} MB\n{n} req.", ha="center", fontsize=7.5)
axs[1].set_yscale("log"); axs[1].set_ylim(0.01, 1000)
axs[1].set_xticks(x, [j.replace(", ", ",\n") for j in jobs])
axs[1].set_ylabel("Data transferred (MB, log scale)")
axs[1].legend(frameon=False, fontsize=8, loc="upper left")
axs[1].set_title("Reading over HTTP with range requests", fontsize=9)
fig.savefig(OUT / "05-cog-size-and-range-reads.png", facecolor="white")
plt.close(fig)
print("saved 5 figures to output/")