Source code: Query Overture Maps with DuckDB spatial and GeoParquet in Python

"""Query Overture Maps with DuckDB spatial and GeoParquet: Bristol's buildings and places.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/geoparquet-duckdb-overture-maps/
Part 1 of "Cloud-native vector GIS with Python". Part 2 (vector-tiles-pmtiles-maplibre) turns the
GeoParquet files this script writes to output/ into PMTiles vector tiles.

Data: Overture Maps Foundation, latest release (read from the STAC catalogue at run time).
Buildings, water: ODbL-1.0 (© OpenStreetMap contributors, Microsoft ML Buildings, ...).
Places: CDLA-Permissive-2.0 / Apache-2.0 / CC0-1.0, depending on the source of each record.
"""
import json
import statistics
import sys
import time
from pathlib import Path

import duckdb
import geopandas as gpd
import matplotlib
import pandas as pd
import pyarrow
import pyarrow.parquet as pq
import pyogrio
import requests
import shapely

matplotlib.use("Agg")
import matplotlib.pyplot as plt  # noqa: E402
from matplotlib.patches import Patch  # noqa: E402

HERE = Path(__file__).resolve().parent
DATA = HERE / "data"      # downloads and benchmark files (not committed)
OUT = HERE / "output"     # figures, console log and the GeoParquet files used in Part 2
DATA.mkdir(exist_ok=True)
OUT.mkdir(exist_ok=True)

BBOX = (-2.66, 51.42, -2.52, 51.50)   # xmin, ymin, xmax, ymax in lon/lat: central Bristol, about 9.7 x 8.9 km
CENTRE = (-2.600, 51.452, -2.578, 51.462)  # Broadmead and Cabot Circus shopping area, for the zoomed map
H3_RES = 9                              # H3 resolution 9: hexagons of about 0.1 km²
BNG = "EPSG:27700"                      # British National Grid, for areas in m²

plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})


def section(title: str) -> None:
    print(f"\n=== {title} ===")


def timed(fn, repeat: int = 1):
    """Run fn() `repeat` times; return the last result and the median wall-clock time in seconds."""
    times = []
    for _ in range(repeat):
        t0 = time.perf_counter()
        result = fn()
        times.append(time.perf_counter() - t0)
    return result, statistics.median(times)


def size_mb(path: Path) -> float:
    """File size in MB; for a Shapefile, add up all its sidecar files."""
    files = list(path.parent.glob(path.stem + ".*")) if path.suffix == ".shp" else [path]
    return sum(f.stat().st_size for f in files) / 1e6


section("versions")
print(f"python {sys.version.split()[0]} | duckdb {duckdb.__version__} | geopandas {gpd.__version__} | "
      f"pyogrio {pyogrio.__version__} (GDAL {pyogrio.__gdal_version_string__}) | pyarrow {pyarrow.__version__} | "
      f"shapely {shapely.__version__}")

# ---------------------------------------------------------------------------- 1. latest release
section("Overture release")
catalog = requests.get("https://stac.overturemaps.org/catalog.json", timeout=60).json()
RELEASE = catalog["latest"]
release_cat = requests.get(f"https://stac.overturemaps.org/{RELEASE}/catalog.json", timeout=60).json()
print(f"latest release: {RELEASE} (schema {release_cat['schema:version']})")
print("themes:", ", ".join(l["title"] for l in release_cat["links"] if l["rel"] == "child"))

con = duckdb.connect()
con.sql("INSTALL spatial; LOAD spatial; INSTALL httpfs; LOAD httpfs; INSTALL h3 FROM community; LOAD h3;")
con.sql("SET s3_region = 'us-west-2'")  # the bucket is public: no credentials needed

# ---------------------------------------------------------------------------- 2. find the files with STAC
section("STAC index: which files cover Bristol?")
con.sql(f"""
    CREATE TABLE stac AS
    SELECT collection, num_rows, assets.aws."file:size" AS bytes,
           assets.aws.alternate.s3.href AS s3_url, bbox
    FROM read_parquet('https://stac.overturemaps.org/{RELEASE}/collections.parquet')
""")
xmin, ymin, xmax, ymax = BBOX
OVERLAPS = f"bbox.xmin <= {xmax} AND bbox.xmax >= {xmin} AND bbox.ymin <= {ymax} AND bbox.ymax >= {ymin}"
print(con.sql(f"""
    SELECT collection, count(*) AS files, sum(num_rows)::BIGINT AS features,
           round(sum(bytes) / 1e9, 1) AS total_gb,
           count(*) FILTER (WHERE {OVERLAPS}) AS files_over_bristol
    FROM stac WHERE collection IN ('building', 'place', 'water')
    GROUP BY ALL ORDER BY collection
""").df().to_string(index=False))


def files_for(collection: str) -> list[str]:
    return [r[0] for r in con.sql(f"SELECT s3_url FROM stac WHERE collection = '{collection}' AND {OVERLAPS}").fetchall()]


# ---------------------------------------------------------------------------- 3. bbox query: naive vs pruned
section("bbox query: whole theme vs STAC-pruned file list")
# The pruned query runs first, so DuckDB's cache of remote file data cannot flatter it.
n_stac, t_stac = timed(lambda: con.sql(
    f"SELECT count(*) FROM read_parquet({files_for('building')}) WHERE {OVERLAPS}").fetchone()[0])
print(f"STAC-pruned file list, bbox filter  : {n_stac:,} buildings in {t_stac:.1f} s")
glob = f"s3://overturemaps-us-west-2/release/{RELEASE}/theme=buildings/type=building/*.parquet"
n_glob, t_glob = timed(lambda: con.sql(f"SELECT count(*) FROM read_parquet('{glob}') WHERE {OVERLAPS}").fetchone()[0])
n_files = con.sql("SELECT count(*) FROM stac WHERE collection = 'building'").fetchone()[0]
print(f"all {n_files} building files, bbox filter : {n_glob:,} buildings in {t_glob:.1f} s "
      f"({t_glob / t_stac:.0f}x slower)")

# ---------------------------------------------------------------------------- 4. download to GeoParquet
section("download Bristol to local GeoParquet")
queries = {
    # `sources` lists who supplied what: an empty `property` means the whole feature,
    # '/properties/height' means only the height came from that source.
    "buildings": f"""
        SELECT id, names.primary AS name, subtype, class, height, num_floors,
               list_filter(sources, s -> s.property = '')[1].dataset AS source,
               list_filter(sources, s -> s.property = '')[1].license AS licence,
               list_filter(sources, s -> s.property = '/properties/height')[1].dataset AS height_source,
               geometry
        FROM read_parquet({files_for('building')})
        WHERE {OVERLAPS}""",
    "places": f"""
        SELECT id, names.primary AS name, basic_category, taxonomy.hierarchy[1] AS category_group,
               round(confidence, 3) AS confidence,
               list_filter(sources, s -> s.property = '')[1].dataset AS source,
               list_filter(sources, s -> s.property = '')[1].license AS licence, geometry
        FROM read_parquet({files_for('place')})
        WHERE {OVERLAPS}""",
    # Water polygons can be huge (the Severn estuary), so clip them to the box.
    "water": f"""
        SELECT id, names.primary AS name, subtype, class,
               ST_Intersection(geometry, ST_MakeEnvelope({xmin}, {ymin}, {xmax}, {ymax})) AS geometry
        FROM read_parquet({files_for('water')})
        WHERE {OVERLAPS} AND ST_GeometryType(geometry) IN ('POLYGON', 'MULTIPOLYGON')""",
}
for name, sql in queries.items():
    path = DATA / f"bristol_{name}_raw.parquet"
    _, t = timed(lambda: con.sql(f"COPY ({sql}) TO '{path.as_posix()}' (FORMAT parquet, COMPRESSION zstd)"))
    con.sql(f"CREATE OR REPLACE TABLE {name} AS SELECT * FROM '{path.as_posix()}'")
    n = con.sql(f"SELECT count(*) FROM {name}").fetchone()[0]
    print(f"{name:<9} {n:>7,} features  {size_mb(path):5.1f} MB  {t:5.1f} s")

# ---------------------------------------------------------------------------- 5. what is in the data?
section("sources and licences")
print(con.sql("""
    SELECT 'buildings' AS theme, source, licence, count(*) AS n FROM buildings GROUP BY ALL
    UNION ALL
    SELECT 'places', source, licence, count(*) FROM places GROUP BY ALL
    ORDER BY theme, n DESC
""").df().to_string(index=False))

section("buildings: subtype and height coverage")
print(con.sql("""
    SELECT coalesce(subtype, '(none)') AS subtype, count(*) AS n,
           round(100 * count(*) / sum(count(*)) OVER (), 1) AS pct,
           round(100 * count(height) / count(*), 1) AS pct_with_height
    FROM buildings GROUP BY ALL ORDER BY n DESC LIMIT 8
""").df().to_string(index=False))
print(con.sql("""
    SELECT CASE WHEN height IS NULL THEN '(no height)' ELSE coalesce(height_source, source) END AS height_from,
           count(*) AS n, round(median(height), 1) AS median_height_m
    FROM buildings GROUP BY ALL ORDER BY n DESC
""").df().to_string(index=False))

section("places: top-level categories")
print(con.sql("""
    SELECT coalesce(category_group, '(none)') AS category_group, count(*) AS n,
           round(median(confidence), 2) AS median_confidence
    FROM places GROUP BY ALL ORDER BY n DESC
""").df().to_string(index=False))

# ---------------------------------------------------------------------------- 6. spatial SQL
section("spatial SQL: footprint areas in British National Grid")
con.sql(f"""
    ALTER TABLE buildings ADD COLUMN area_m2 DOUBLE;
    UPDATE buildings SET area_m2 = round(ST_Area(ST_Transform(geometry, 'OGC:CRS84', '{BNG}')), 1);
""")
print(con.sql("""
    SELECT count(*) AS buildings, round(sum(area_m2) / 1e6, 2) AS footprint_km2,
           round(median(area_m2)) AS median_m2, round(max(area_m2)) AS max_m2
    FROM buildings
""").df().to_string(index=False))
print(con.sql("""
    SELECT name, subtype, round(area_m2) AS area_m2 FROM buildings
    WHERE name IS NOT NULL ORDER BY area_m2 DESC LIMIT 5
""").df().to_string(index=False))

section("spatial SQL: which places are inside a building?")
_, t_join = timed(lambda: con.sql("""
    CREATE OR REPLACE TABLE place_building AS
    SELECT p.id AS place_id, b.id AS building_id
    FROM places p JOIN buildings b ON ST_Intersects(b.geometry, p.geometry)
"""))
n_places = con.sql("SELECT count(*) FROM places").fetchone()[0]
n_inside = con.sql("SELECT count(DISTINCT place_id) FROM place_building").fetchone()[0]
print(f"point-in-polygon join of {n_places:,} places x {n_stac:,} buildings: {t_join:.2f} s")
print(f"places inside a building footprint: {n_inside:,} ({100 * n_inside / n_places:.1f}%)")
print(con.sql("""
    SELECT b.name, b.subtype, round(b.area_m2) AS area_m2, count(*) AS places
    FROM place_building pb JOIN buildings b ON b.id = pb.building_id
    GROUP BY ALL ORDER BY places DESC LIMIT 6
""").df().to_string(index=False))

section(f"spatial SQL: aggregate to H3 resolution {H3_RES} hexagons")
_, t_h3 = timed(lambda: con.sql(f"""
    CREATE OR REPLACE TABLE hex AS
    WITH b AS (
        SELECT h3_latlng_to_cell(ST_Y(ST_Centroid(geometry)), ST_X(ST_Centroid(geometry)), {H3_RES}) AS h3,
               count(*) AS buildings, sum(area_m2) AS footprint_m2,
               count(*) FILTER (WHERE subtype = 'residential') AS residential
        FROM buildings GROUP BY h3),
    p AS (
        SELECT h3_latlng_to_cell(ST_Y(geometry), ST_X(geometry), {H3_RES}) AS h3,
               count(*) AS places,
               count(*) FILTER (WHERE category_group = 'food_and_drink') AS food_drink
        FROM places GROUP BY h3)
    SELECT h3_h3_to_string(h3) AS h3, coalesce(buildings, 0) AS buildings, coalesce(residential, 0) AS residential,
           coalesce(places, 0) AS places, coalesce(food_drink, 0) AS food_drink,
           round(100 * coalesce(footprint_m2, 0) / h3_cell_area(h3, 'm^2'), 1) AS built_pct,
           ST_GeomFromText(h3_cell_to_boundary_wkt(h3)) AS geometry
    FROM b FULL JOIN p USING (h3)
    -- keep whole hexagons only: edge cells would be cut by the download box
    WHERE ST_Within(ST_GeomFromText(h3_cell_to_boundary_wkt(h3)), ST_MakeEnvelope({xmin}, {ymin}, {xmax}, {ymax}))
"""))
print(f"built {con.sql('SELECT count(*) FROM hex').fetchone()[0]} hexagons in {t_h3:.2f} s")
print(con.sql("""
    SELECT h3, buildings, built_pct, places, food_drink FROM hex ORDER BY built_pct DESC LIMIT 5
""").df().to_string(index=False))
print(con.sql("""
    SELECT h3, buildings, built_pct, places, food_drink FROM hex ORDER BY food_drink DESC LIMIT 5
""").df().to_string(index=False))
print(con.sql("""
    SELECT round(corr(built_pct, places), 2) AS r_built_places,
           round(corr(residential::DOUBLE / greatest(buildings, 1), food_drink), 2) AS r_residential_share_food
    FROM hex WHERE buildings > 0
""").df().to_string(index=False))

# ---------------------------------------------------------------------------- 7. GeoPandas and GeoParquet
section("DuckDB -> GeoPandas")


def to_gdf(table: str) -> gpd.GeoDataFrame:
    """DuckDB table -> GeoDataFrame, passing geometry as WKB."""
    df = con.sql(f"SELECT * EXCLUDE geometry, ST_AsWKB(geometry) AS wkb FROM {table}").df()
    geom = shapely.from_wkb(df.pop("wkb").map(bytes))
    return gpd.GeoDataFrame(df, geometry=geom, crs="OGC:CRS84")


gdfs, t_gdf = timed(lambda: {t: to_gdf(t) for t in ("buildings", "places", "water", "hex")})
buildings, places, water, hexes = (gdfs[k] for k in ("buildings", "places", "water", "hex"))
print(f"4 tables to GeoDataFrames in {t_gdf:.2f} s; buildings: {len(buildings):,} rows, crs {buildings.crs.name}")
print(buildings.drop(columns="geometry").head(3).to_string(index=False))

gdf_raw = gpd.read_parquet(DATA / "bristol_buildings_raw.parquet")
print(f"gpd.read_parquet on DuckDB's file: {len(gdf_raw):,} rows, geometry column '{gdf_raw.geometry.name}', "
      f"crs {gdf_raw.crs.to_string() if gdf_raw.crs else None}")

section("write GeoParquet 1.1 with GeoPandas")
for name, gdf in gdfs.items():
    gdf.to_parquet(OUT / f"bristol_{name}.parquet", compression="zstd",
                   schema_version="1.1.0", write_covering_bbox=True)
    print(f"output/bristol_{name}.parquet  {len(gdf):>7,} rows  {size_mb(OUT / f'bristol_{name}.parquet'):5.2f} MB")

geo = json.loads(pq.read_metadata(OUT / "bristol_buildings.parquet").metadata[b"geo"])
col = geo["columns"]["geometry"]
print("'geo' metadata (GeoPandas):", json.dumps({
    "version": geo["version"], "primary_column": geo["primary_column"], "encoding": col["encoding"],
    "geometry_types": col["geometry_types"], "bbox": [round(v, 4) for v in col["bbox"]],
    "covering": col.get("covering"), "crs": "PROJJSON " + col["crs"]["id"]["authority"] + ":" + str(col["crs"]["id"]["code"])}))
geo_ddb = json.loads(pq.read_metadata(DATA / "bristol_buildings_raw.parquet").metadata[b"geo"])
print("'geo' metadata (DuckDB):  ", json.dumps({k: geo_ddb[k] for k in ("version", "primary_column")}),
      "| keys of the geometry column:", sorted(geo_ddb["columns"]["geometry"]))

# ---------------------------------------------------------------------------- 8. format benchmark
section("format benchmark: buildings")
bench = buildings.drop(columns=["licence", "height_source"])
small = (-2.596, 51.447, -2.582, 51.456)   # ~1 km box: Temple Meads to Redcliffe
BENCH = DATA / "bench"
BENCH.mkdir(exist_ok=True)
formats = {
    "GeoParquet (DuckDB)": BENCH / "duckdb.parquet",
    "GeoParquet (GeoPandas)": BENCH / "geopandas.parquet",
    "FlatGeobuf": BENCH / "buildings.fgb",
    "GeoPackage": BENCH / "buildings.gpkg",
    "Shapefile": BENCH / "buildings.shp",
}
for p in BENCH.iterdir():
    p.unlink()
con.sql("CREATE OR REPLACE TABLE bench AS SELECT * EXCLUDE (licence, height_source) FROM buildings")
writers = {
    "GeoParquet (DuckDB)": lambda p: con.sql(f"COPY bench TO '{p.as_posix()}' (FORMAT parquet, COMPRESSION zstd)"),
    "GeoParquet (GeoPandas)": lambda p: bench.to_parquet(p, compression="zstd", write_covering_bbox=True,
                                                         schema_version="1.1.0"),
    "FlatGeobuf": lambda p: pyogrio.write_dataframe(bench, p, driver="FlatGeobuf"),
    "GeoPackage": lambda p: pyogrio.write_dataframe(bench, p, driver="GPKG"),
    "Shapefile": lambda p: pyogrio.write_dataframe(bench, p, driver="ESRI Shapefile"),
}
readers = {
    "GeoParquet (DuckDB)": lambda p, bbox=None: gpd.read_parquet(p) if bbox is None else con.sql(
        f"SELECT count(*) FROM '{p.as_posix()}' WHERE ST_Intersects(geometry, ST_MakeEnvelope{bbox})").fetchone(),
    "GeoParquet (GeoPandas)": lambda p, bbox=None: gpd.read_parquet(p, bbox=bbox),
}
rows = []
for fmt, path in formats.items():
    _, t_write = timed(lambda: writers[fmt](path), repeat=3)
    read = readers.get(fmt, lambda p, bbox=None: pyogrio.read_dataframe(p, bbox=bbox, use_arrow=True))
    full, t_read = timed(lambda: read(path), repeat=3)
    sub, t_bbox = timed(lambda: read(path, bbox=small), repeat=3)
    n_sub = sub[0] if isinstance(sub, tuple) else len(sub)
    rows.append({"format": fmt, "size_MB": round(size_mb(path), 1), "write_s": round(t_write, 2),
                 "read_s": round(t_read, 2), "bbox_read_s": round(t_bbox, 3), "bbox_rows": n_sub})
bench_df = pd.DataFrame(rows)
print(bench_df.to_string(index=False))
print("Shapefile field names:", list(pyogrio.read_info(formats["Shapefile"])["fields"]))
print("note: the DuckDB bbox query counts rows with ST_Intersects; the others return a GeoDataFrame")

# ---------------------------------------------------------------------------- 9. figures
section("figures")
b27, w27, p27, h27 = (g.to_crs(BNG) for g in (buildings, water, places, hexes))
X0, Y0, X1, Y1 = gpd.GeoSeries([shapely.box(*BBOX)], crs="OGC:CRS84").to_crs(BNG).total_bounds

# Figure 1: buildings by subtype
SUBTYPE_COLOURS = {"residential": "#e8a33c", "commercial": "#c21f3a", "education": "#5e3c99",
                   "industrial": "#6b6b6b", "outbuilding": "#a6761d", "(other)": "#1b9e77", "(none)": "#9aa3a8"}
main = list(SUBTYPE_COLOURS)[:-2]
b27["group"] = b27["subtype"].where(b27["subtype"].isin(main), "(other)").where(b27["subtype"].notna(), "(none)")
fig, ax = plt.subplots(figsize=(10, 9))
w27.plot(ax=ax, color="#9ecae1", linewidth=0)
for g, colour in SUBTYPE_COLOURS.items():
    b27[b27["group"] == g].plot(ax=ax, color=colour, linewidth=0)
ax.set_xlim(X0, X1)
ax.set_ylim(Y0, Y1)
ax.set_axis_off()
counts = b27["group"].value_counts()
ax.legend(handles=[Patch(color=c, label=f"{g} ({counts.get(g, 0):,})") for g, c in SUBTYPE_COLOURS.items()]
          + [Patch(color="#9ecae1", label="water")], loc="lower left", framealpha=0.9, title="building subtype")
ax.set_title(f"{len(b27):,} Overture buildings in central Bristol (release {RELEASE})", loc="left")
ax.text(1, -0.01, "© OpenStreetMap contributors, Microsoft, Overture Maps Foundation (ODbL)",
        transform=ax.transAxes, ha="right", va="top", fontsize=7, color="#4a5553")
fig.savefig(OUT / "01-bristol-buildings-subtype.png")
plt.close(fig)

# Figure 2: places by source/licence and by category
LIC_COLOURS = {"CDLA-Permissive-2.0": "#1543b0", "Apache-2.0": "#e8a33c", "CC0-1.0": "#1b9e77"}
src = places.groupby(["source", "licence"]).size().reset_index(name="n").sort_values("n")
cat = places["category_group"].fillna("(none)").value_counts().sort_values()
fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 4.2), gridspec_kw={"width_ratios": [1, 1.3]})
a1.barh(src["source"], src["n"], color=[LIC_COLOURS.get(l, "#9aa3a8") for l in src["licence"]])
for y, n in enumerate(src["n"]):
    a1.text(n, y, f" {n:,}", va="center", fontsize=8)
a1.set_xlim(0, src["n"].max() * 1.2)
a1.set_title("Places by source (colour = licence)", loc="left")
a1.legend(handles=[Patch(color=c, label=l) for l, c in LIC_COLOURS.items()], loc="lower right")
a2.barh(cat.index.str.replace("_", " "), cat.values, color="#5e3c99")
for y, n in enumerate(cat.values):
    a2.text(n, y, f" {n:,}", va="center", fontsize=8)
a2.set_xlim(0, cat.max() * 1.15)
a2.set_title("Places by top-level category", loc="left")
for a in (a1, a2):
    a.spines[["top", "right"]].set_visible(False)
fig.tight_layout()
fig.savefig(OUT / "02-places-sources-categories.png")
plt.close(fig)

# Figure 3: H3 hexagons
fig, axes = plt.subplots(1, 2, figsize=(12, 5.6), layout="constrained")
for ax, col_, title, cmap in [(axes[0], "built_pct", "Building footprint cover (% of hexagon)", "Spectral_r"),
                              (axes[1], "food_drink", "Food and drink places per hexagon", "magma_r")]:
    w27.plot(ax=ax, color="#9ecae1", linewidth=0)
    h27.plot(ax=ax, column=col_, cmap=cmap, linewidth=0.2, edgecolor="white", alpha=0.9, legend=True,
             legend_kwds={"shrink": 0.7})
    ax.set_xlim(X0, X1)
    ax.set_ylim(Y0, Y1)
    ax.set_axis_off()
    ax.set_title(title, loc="left")
fig.suptitle(f"{len(h27)} H3 resolution-{H3_RES} hexagons over central Bristol, aggregated in DuckDB", x=0.02,
             ha="left", fontsize=10)
fig.savefig(OUT / "03-h3-hexagons.png")
plt.close(fig)

# Figure 4: zoom on the shopping centre: places inside vs outside buildings
cx0, cy0, cx1, cy1 = gpd.GeoSeries([shapely.box(*CENTRE)], crs="OGC:CRS84").to_crs(BNG).total_bounds
per_b = con.sql("SELECT building_id AS id, count(*) AS places FROM place_building GROUP BY 1").df()
bz = b27.cx[cx0:cx1, cy0:cy1].merge(per_b, on="id", how="left").fillna({"places": 0})
pz = p27.cx[cx0:cx1, cy0:cy1].copy()
inside_ids = set(con.sql("SELECT place_id FROM place_building").df()["place_id"])
pz["inside"] = pz["id"].isin(inside_ids)
fig, ax = plt.subplots(figsize=(10, 6.2))
w27.plot(ax=ax, color="#9ecae1", linewidth=0)
bz[bz["places"] == 0].plot(ax=ax, color="#d9dcdb", edgecolor="#b5bab8", linewidth=0.3)
bz[bz["places"] > 0].plot(ax=ax, column="places", cmap="YlOrRd", edgecolor="#7a7f7d", linewidth=0.3,
                          legend=True, legend_kwds={"label": "places inside building", "shrink": 0.7},
                          vmax=max(bz["places"].quantile(0.99), 5))
pz[pz["inside"]].plot(ax=ax, color="#0d1315", markersize=3)
pz[~pz["inside"]].plot(ax=ax, color="#1543b0", markersize=5, marker="x", linewidth=0.8)
ax.set_xlim(cx0, cx1)
ax.set_ylim(cy0, cy1)
ax.set_axis_off()
ax.legend(handles=[plt.Line2D([], [], color="#0d1315", marker="o", ls="", ms=3, label="place inside a building"),
                   plt.Line2D([], [], color="#1543b0", marker="x", ls="", ms=5, label="place outside every building")],
          loc="lower left", framealpha=0.9)
ax.set_title(f"Broadmead and Cabot Circus: {int(pz['inside'].sum()):,} of {len(pz):,} places fall inside a footprint",
             loc="left")
fig.savefig(OUT / "04-places-in-buildings.png")
plt.close(fig)

# Figure 5: format benchmark
fig, axes = plt.subplots(1, 4, figsize=(13, 3.6), sharey=True)
colours = ["#1543b0", "#5b8bd9", "#1b9e77", "#e8a33c", "#c21f3a"]
for ax, col_, title in zip(axes, ["size_MB", "write_s", "read_s", "bbox_read_s"],
                           ["File size (MB)", "Write (s)", "Read all (s)", "Read 1 km box (s)"]):
    ax.barh(bench_df["format"], bench_df[col_], color=colours)
    for y, v in enumerate(bench_df[col_]):
        ax.text(v, y, f" {v:g}", va="center", fontsize=8)
    ax.set_xlim(0, bench_df[col_].max() * 1.35)
    ax.set_title(title, loc="left")
    ax.spines[["top", "right"]].set_visible(False)
axes[0].invert_yaxis()
fig.suptitle(f"{len(bench):,} building polygons, {len(bench.columns) - 1} attribute columns (median of 3 runs)", x=0.01, ha="left", fontsize=10)
fig.tight_layout()
fig.savefig(OUT / "05-format-benchmark.png")
plt.close(fig)

print("wrote", ", ".join(sorted(p.name for p in OUT.glob("*.png"))))