"""Turn Overture GeoParquet into PMTiles vector tiles with GDAL, and build a static MapLibre GL JS map.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/vector-tiles-pmtiles-maplibre/
Part 2 of "Cloud-native vector GIS with Python". It reads the Bristol GeoParquet files written by Part 1
(geoparquet-duckdb-overture-maps). If they are missing, it runs the Part 1 script first.

Tiles are written by GDAL's PMTiles driver (GDAL >= 3.8), which ships inside the pyogrio wheel, so this
runs on Windows without Tippecanoe. Then:
    python serve.py                 # open http://localhost:8000/ (a server with HTTP Range support)
    node screenshot_map.mjs         # optional: Playwright screenshots of the running map

Data: Overture Maps Foundation. Buildings and water: ODbL-1.0 (© OpenStreetMap contributors, Microsoft).
Places: CDLA-Permissive-2.0 / Apache-2.0 / CC0-1.0, depending on the source of each record.
"""
import json
import math
import subprocess
import sys
import time
from pathlib import Path

import geopandas as gpd
import matplotlib
import pandas as pd
import pyogrio
import requests
import shapely
from pmtiles.reader import MmapSource, Reader, all_tiles

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

HERE = Path(__file__).resolve().parent
PART1 = HERE.parent / "geoparquet-duckdb-overture-maps"
DATA = HERE / "data"          # experiment tilesets (not committed)
OUT = HERE / "output"         # figures, console log, and the finished web map in output/web/
WEB = OUT / "web"
for d in (DATA, OUT, WEB, WEB / "lib"):
    d.mkdir(parents=True, exist_ok=True)

RELEASE = "2026-09-23.1"      # the Overture release Part 1 downloaded
MAPLIBRE = "5.24.0"           # maplibre-gl, self-hosted from output/web/lib/
PMTILES_JS = "4.5.0"          # pmtiles (JavaScript), inlined into the HTML page
ATTRIBUTION = ('© <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors, '
               '<a href="https://docs.overturemaps.org/attribution/">Overture Maps Foundation</a> '
               f'(release {RELEASE})')
# The copy of the map embedded in the tutorial loads its tiles and libraries from here on technicalgis.com
SITE_BASE = "/demos/vector-tiles-pmtiles-maplibre/"

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


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


def write_pmtiles(gdf: gpd.GeoDataFrame, path: Path, layer: str, **options) -> float:
    """Write one layer to a new PMTiles file with GDAL; return the time taken in seconds."""
    path.unlink(missing_ok=True)
    t0 = time.perf_counter()
    # encoding="UTF-8": otherwise pyogrio may encode names with the Windows code page and fail on characters
    # such as Welsh or Polish letters
    pyogrio.write_dataframe(gdf, path, driver="PMTiles", layer=layer, encoding="UTF-8",
                            **{k: str(v) for k, v in options.items()})
    return time.perf_counter() - t0


def tile_stats(path: Path) -> pd.DataFrame:
    """One row per tile: z, x, y and compressed size in bytes, read from the PMTiles directory."""
    with open(path, "rb") as f:
        reader = Reader(MmapSource(f))
        rows = [(z, x, y, len(data)) for (z, x, y), data in all_tiles(reader.get_bytes)]
    return pd.DataFrame(rows, columns=["z", "x", "y", "bytes"])


def tile_bounds(z: int, x: int, y: int) -> tuple[float, float, float, float]:
    """Lon/lat bounds of an XYZ web-mercator tile."""
    n = 2 ** z
    lat = lambda t: math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * t / n))))  # noqa: E731
    return x / n * 360 - 180, lat(y + 1), (x + 1) / n * 360 - 180, lat(y)


section("versions")
print(f"python {sys.version.split()[0]} | geopandas {gpd.__version__} | pyogrio {pyogrio.__version__} "
      f"(GDAL {pyogrio.__gdal_version_string__}) | shapely {shapely.__version__}")
drivers = pyogrio.list_drivers()
print(f"GDAL vector-tile drivers: PMTiles={drivers.get('PMTiles')}  MVT={drivers.get('MVT')}")

# ---------------------------------------------------------------------------- 1. load Part 1's GeoParquet
section("load Part 1 GeoParquet")
names = ["buildings", "places", "water", "hex"]
paths = {n: PART1 / "output" / f"bristol_{n}.parquet" for n in names}
if not all(p.exists() for p in paths.values()):
    print("Part 1 output not found: running geoparquet-duckdb-overture-maps/overture_duckdb.py first")
    subprocess.run([sys.executable, str(PART1 / "overture_duckdb.py")], check=True, cwd=PART1)
layers = {n: gpd.read_parquet(p) for n, p in paths.items()}
for n, g in layers.items():
    print(f"{n:<9} {len(g):>7,} features  {paths[n].stat().st_size / 1e6:5.2f} MB GeoParquet  crs {g.crs.to_string()}")

# Keep only what the web map needs: fewer attributes mean smaller tiles (measured in section 3)
b = layers["buildings"]
buildings = gpd.GeoDataFrame({"subtype": b["subtype"], "height": b["height"].round(1)}, geometry=b.geometry)
p = layers["places"]
places = gpd.GeoDataFrame({"name": p["name"], "category": p["basic_category"], "group": p["category_group"],
                           "confidence": p["confidence"].round(2), "source": p["source"], "licence": p["licence"]},
                          geometry=p.geometry)
water = layers["water"][["subtype", "geometry"]]
hexes = layers["hex"][["h3", "buildings", "places", "food_drink", "built_pct", "geometry"]]

# ---------------------------------------------------------------------------- 2. zoom-level trade-off
section("buildings: maximum zoom vs size (MINZOOM=12, GDAL defaults otherwise)")
rows, per_zoom = [], {}
for maxzoom in [12, 13, 14, 15, 16]:
    path = DATA / f"buildings_z12-{maxzoom}.pmtiles"
    secs = write_pmtiles(buildings, path, "buildings", MINZOOM=12, MAXZOOM=maxzoom)
    stats = tile_stats(path)
    per_zoom[maxzoom] = stats
    rows.append({"zooms": f"12-{maxzoom}", "tiles": len(stats), "size_MB": round(path.stat().st_size / 1e6, 2),
                 "write_s": round(secs, 1), "largest_tile_KB": round(stats["bytes"].max() / 1e3)})
zoom_df = pd.DataFrame(rows)
print(zoom_df.to_string(index=False))

by_z = per_zoom[16].groupby("z")["bytes"].agg(
    tiles="size", total_KB=lambda s: round(s.sum() / 1e3), mean_KB=lambda s: round(s.mean() / 1e3, 1),
    max_KB=lambda s: round(s.max() / 1e3))
tile_m = {z: 40_075_016.7 * math.cos(math.radians(51.46)) / 2 ** z for z in by_z.index}
by_z["tile_width_m"] = [round(tile_m[z]) for z in by_z.index]
by_z["grid_step_m"] = [round(tile_m[z] / 4096, 2) for z in by_z.index]   # EXTENT=4096 grid units per tile
print("per zoom level in the z12-16 archive (compressed tile sizes):")
print(by_z.to_string())

# What each zoom level holds: GDAL keeps tiles under MAX_SIZE (500 KB) by coarsening or dropping features
for z in [12, 13, 14, 16]:
    g = pyogrio.read_dataframe(DATA / "buildings_z12-16.pmtiles", layer="buildings", ZOOM_LEVEL=str(z))
    print(f"z{z}: {len(g):,} feature pieces read back (a building cut by a tile edge counts once per tile)")

# ---------------------------------------------------------------------------- 3. attributes and precision
section("buildings at z14 only: attributes, EXTENT and SIMPLIFICATION")
variants = {
    "geometry only": (buildings[["geometry"]], {}),
    "subtype + height (used)": (buildings, {}),
    "all 10 Part 1 columns": (b, {}),
    "subtype + height, EXTENT=1024": (buildings, {"EXTENT": 1024}),
    "subtype + height, SIMPLIFICATION=4": (buildings, {"SIMPLIFICATION": 4}),
}
rows = []
for i, (label, (gdf, opts)) in enumerate(variants.items()):
    path = DATA / f"buildings_z14_variant{i}.pmtiles"
    secs = write_pmtiles(gdf, path, "buildings", MINZOOM=14, MAXZOOM=14, **opts)
    rows.append({"variant": label, "size_MB": round(path.stat().st_size / 1e6, 2), "write_s": round(secs, 1)})
attr_df = pd.DataFrame(rows)
print(attr_df.to_string(index=False))

# ---------------------------------------------------------------------------- 4. the tilesets for the web map
section("final tilesets for the web map")
final = {  # layer: (data, minzoom, maxzoom). MapLibre overzooms the last level, so z14 tiles still draw at z19
    "hex": (hexes, 10, 13),
    "water": (water, 10, 14),
    "buildings": (buildings, 14, 14),
    "places": (places, 14, 14),
}
rows = []
for name, (gdf, z0, z1) in final.items():
    path = WEB / f"{name}.pmtiles"
    secs = write_pmtiles(gdf, path, name, MINZOOM=z0, MAXZOOM=z1, NAME=f"Bristol {name}",
                         DESCRIPTION=f"Overture Maps {RELEASE}, central Bristol")
    rows.append({"file": path.name, "features": len(gdf), "zooms": f"{z0}-{z1}", "tiles": len(tile_stats(path)),
                 "size_MB": round(path.stat().st_size / 1e6, 2), "write_s": round(secs, 1)})
final_df = pd.DataFrame(rows)
print(final_df.to_string(index=False))
print(f"total: {final_df['size_MB'].sum():.2f} MB in {final_df['tiles'].sum()} tiles "
      f"(the four GeoParquet inputs are {sum(p.stat().st_size for p in paths.values()) / 1e6:.1f} MB)")
meta = pyogrio.read_info(WEB / "buildings.pmtiles", layer="buildings")["dataset_metadata"]
print("buildings.pmtiles metadata:", {k: meta[k] for k in ["name", "minzoom", "maxzoom", "bounds", "format"]})

# ---------------------------------------------------------------------------- 5. the MapLibre page
section("MapLibre GL JS page")
for url, dest in [(f"https://unpkg.com/maplibre-gl@{MAPLIBRE}/dist/maplibre-gl.js", "maplibre-gl.js"),
                  (f"https://unpkg.com/maplibre-gl@{MAPLIBRE}/dist/maplibre-gl.css", "maplibre-gl.css"),
                  (f"https://unpkg.com/maplibre-gl@{MAPLIBRE}/LICENSE.txt", "LICENSE-maplibre-gl.txt"),
                  (f"https://unpkg.com/pmtiles@{PMTILES_JS}/dist/pmtiles.js", "pmtiles.js")]:
    target = WEB / "lib" / dest
    if not target.exists():
        response = requests.get(url, timeout=60)
        response.raise_for_status()
        target.write_bytes(response.content)
pmtiles_js = (WEB / "lib" / "pmtiles.js").read_text(encoding="utf-8")

xmin, ymin, xmax, ymax = layers["hex"].total_bounds
GROUP_COLOURS = {"food_and_drink": "#c21f3a", "shopping": "#e8a33c", "services_and_business": "#1543b0",
                 "health_care": "#1b9e77", "arts_and_entertainment": "#7b3294", "education": "#5e3c99"}
place_colour = ["match", ["get", "group"], *[v for kv in GROUP_COLOURS.items() for v in kv], "#6b6b6b"]


def style(base: str) -> dict:
    """MapLibre style JSON: one vector source per PMTiles file, no third-party basemap."""
    sources = {n: {"type": "vector", "url": f"pmtiles://{base}{n}.pmtiles", "minzoom": z0, "maxzoom": z1,
                   "attribution": ATTRIBUTION} for n, (_, z0, z1) in final.items()}
    return {
        "version": 8,
        "sources": sources,
        "layers": [
            {"id": "background", "type": "background", "paint": {"background-color": "#f3f1ec"}},
            {"id": "water", "type": "fill", "source": "water", "source-layer": "water",
             "paint": {"fill-color": "#9ecae1"}},
            {"id": "hex", "type": "fill", "source": "hex", "source-layer": "hex", "maxzoom": 14.5,
             "paint": {"fill-color": ["interpolate", ["linear"], ["get", "places"],
                                      0, "#fff7ec", 25, "#fdd49e", 75, "#fc8d59", 200, "#d7301f", 500, "#7f0000"],
                       "fill-outline-color": "#ffffff",
                       "fill-opacity": ["interpolate", ["linear"], ["zoom"], 13.5, 0.85, 14.5, 0]}},
            {"id": "buildings", "type": "fill", "source": "buildings", "source-layer": "buildings", "minzoom": 14,
             "paint": {"fill-color": ["match", ["get", "subtype"], "residential", "#e8c48f", "commercial", "#e49aa6",
                                      "education", "#b7a5d6", "medical", "#9fd3c0", "#c9cdca"],
                       "fill-outline-color": "#8a908d"}},
            {"id": "places", "type": "circle", "source": "places", "source-layer": "places", "minzoom": 15,
             "paint": {"circle-color": place_colour, "circle-stroke-color": "#ffffff", "circle-stroke-width": 1,
                       "circle-radius": ["interpolate", ["linear"], ["zoom"], 15, 3, 18, 6]}},
        ],
    }


LEGEND = "".join(f'<span><i style="background:{c}"></i>{g.replace("_", " ")}</span>' for g, c in GROUP_COLOURS.items())
PAGE = """<!doctype html>
<html lang="en-GB">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<meta name="robots" content="noindex, follow">
<title>Bristol buildings and places: PMTiles and MapLibre GL JS</title>
<link rel="stylesheet" href='__BASE__lib/maplibre-gl.css'>
<script src='__BASE__lib/maplibre-gl.js'></script>
<!-- pmtiles __PMTILES_VERSION__ (BSD-3-Clause, https://github.com/protomaps/PMTiles), inlined -->
<script>__PMTILES_JS__</script>
<style>
  html, body, #map { margin: 0; height: 100%; font: 13px/1.4 system-ui, sans-serif; }
  .panel { position: absolute; top: 10px; left: 10px; max-width: 260px; background: rgba(255,255,255,.93);
           padding: 8px 10px; border-radius: 6px; box-shadow: 0 1px 4px rgba(0,0,0,.25); }
  .panel b { display: block; margin-bottom: 2px; }
  .panel span { display: inline-flex; align-items: center; margin: 2px 8px 0 0; font-size: 12px; }
  .panel i { width: 10px; height: 10px; border-radius: 50%; margin-right: 4px; display: inline-block; }
  #zoom { color: #555; font-size: 12px; }
  @media (max-width: 480px) { .panel { max-width: 180px; font-size: 12px; } .panel span { display: none; } }
</style>
</head>
<body>
<div id="map"></div>
<div class="panel"><b>Central Bristol from Overture Maps</b>
  <div id="zoom"></div>
  <div>Zoom 10&ndash;14: places per H3 hexagon. Zoom 14+: buildings. Zoom 15+: places (click one).</div>
  __LEGEND__
</div>
<script>
const protocol = new pmtiles.Protocol();
maplibregl.addProtocol("pmtiles", protocol.tile);
const params = new URLSearchParams(location.search);
const map = new maplibregl.Map({
  container: "map",
  style: __STYLE__,
  center: [+(params.get("lon") || -2.5925), +(params.get("lat") || 51.4585)],
  zoom: +(params.get("z") || 12.4),
  minZoom: 10, maxZoom: 19,
  maxBounds: [[__XMIN__ - 0.05, __YMIN__ - 0.03], [__XMAX__ + 0.05, __YMAX__ + 0.03]],
  attributionControl: { compact: window.innerWidth < 640 },
});
map.addControl(new maplibregl.NavigationControl(), "top-right");
map.addControl(new maplibregl.ScaleControl({ unit: "metric" }), "bottom-left");
// MapLibre opens the compact attribution on load; start it collapsed so it doesn't cover the scale bar on phones.
const collapseAttrib = () => document.querySelectorAll(".maplibregl-ctrl-attrib.maplibregl-compact").forEach((a) => {
  a.classList.remove("maplibregl-compact-show"); a.removeAttribute("open");
});
map.once("load", collapseAttrib); map.once("idle", collapseAttrib);
const showZoom = () => document.getElementById("zoom").textContent = "zoom " + map.getZoom().toFixed(1);
map.on("zoom", showZoom); showZoom();
map.on("click", "places", (e) => {
  const p = e.features[0].properties;
  new maplibregl.Popup().setLngLat(e.lngLat).setHTML(
    `<b>${p.name || "(no name)"}</b><br>${(p.category || "").replaceAll("_", " ")}<br>` +
    `source: ${p.source} (${p.licence})<br>confidence: ${p.confidence}`).addTo(map);
});
map.on("mouseenter", "places", () => map.getCanvas().style.cursor = "pointer");
map.on("mouseleave", "places", () => map.getCanvas().style.cursor = "");
window.map = map;
</script>
</body>
</html>
"""


def page(base: str) -> str:
    html = PAGE
    for key, value in {"__BASE__": base, "__PMTILES_VERSION__": PMTILES_JS, "__LEGEND__": LEGEND,
                       "__STYLE__": json.dumps(style(base), ensure_ascii=False), "__XMIN__": f"{xmin:.4f}",
                       "__YMIN__": f"{ymin:.4f}", "__XMAX__": f"{xmax:.4f}", "__YMAX__": f"{ymax:.4f}",
                       "__PMTILES_JS__": pmtiles_js}.items():
        html = html.replace(key, value)
    return html


(WEB / "index.html").write_text(page(""), encoding="utf-8")
(OUT / "map.html").write_text(page(SITE_BASE), encoding="utf-8")   # the copy embedded on technicalgis.com
web_files = sorted(f for f in WEB.rglob("*") if f.is_file())
print(f"output/web/: {len(web_files)} files, {sum(f.stat().st_size for f in web_files) / 1e6:.2f} MB")
for f in web_files:
    print(f"  {f.relative_to(WEB).as_posix():<28} {f.stat().st_size / 1e3:8.1f} KB")
print("serve it with:  python serve.py   then open http://localhost:8000/")

# ---------------------------------------------------------------------------- 6. figures
section("figures")

# Figure 1: zoom ladder
fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 3.8))
a1.bar(zoom_df["zooms"], zoom_df["size_MB"], color="#1543b0")
for i, (s, t) in enumerate(zip(zoom_df["size_MB"], zoom_df["write_s"])):
    a1.text(i, s, f"{s:g} MB\n{t:g} s", ha="center", va="bottom", fontsize=8)
a1.set_ylim(0, zoom_df["size_MB"].max() * 1.25)
a1.set_xlabel("zoom levels in the archive")
a1.set_title("Buildings archive: size and GDAL write time", loc="left")
a2.bar(by_z.index.astype(str), by_z["tiles"], color="#e8a33c")
a2.set_yscale("log")
for i, (n, kb) in enumerate(zip(by_z["tiles"], by_z["max_KB"])):
    a2.text(i, n, f"{n} tiles\nmax {kb} KB", ha="center", va="bottom", fontsize=8)
a2.set_ylim(1, by_z["tiles"].max() * 8)
a2.set_xlabel("zoom level")
a2.set_title("Tiles per zoom level (log scale) and the largest tile", loc="left")
for a in (a1, a2):
    a.spines[["top", "right"]].set_visible(False)
fig.tight_layout()
fig.savefig(OUT / "01-zoom-size-tradeoff.png")
plt.close(fig)

# Figure 2: tile sizes at z14, drawn over the buildings
z14 = tile_stats(WEB / "buildings.pmtiles")
fig, ax = plt.subplots(figsize=(9, 7.4))
water.to_crs(3857).plot(ax=ax, color="#9ecae1", linewidth=0)
norm = LogNorm(vmin=max(z14["bytes"].min() / 1e3, 0.1), vmax=z14["bytes"].max() / 1e3)
cmap = plt.get_cmap("YlOrRd")
for t in z14.itertuples():
    x0, y0, x1, y1 = gpd.GeoSeries([shapely.box(*tile_bounds(t.z, t.x, t.y))], crs=4326).to_crs(3857).total_bounds
    ax.add_patch(Rectangle((x0, y0), x1 - x0, y1 - y0, facecolor=cmap(norm(t.bytes / 1e3)), edgecolor="#555555",
                           linewidth=0.5, alpha=0.6, zorder=1))
    ax.text((x0 + x1) / 2, (y0 + y1) / 2, f"{t.bytes / 1e3:.0f}" if t.bytes >= 1e3 else f"{t.bytes / 1e3:.1f}",
            ha="center", va="center", fontsize=7, zorder=3,
            bbox={"boxstyle": "round,pad=0.15", "facecolor": "white", "edgecolor": "none", "alpha": 0.8})
buildings.to_crs(3857).plot(ax=ax, color="#1d2526", linewidth=0, alpha=0.55, zorder=2)
ax.set_aspect("equal")
ax.autoscale_view()
ax.set_axis_off()
fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax, shrink=0.6, label="compressed tile size (KB)")
ax.set_title(f"{len(z14)} zoom-14 building tiles, {z14['bytes'].sum() / 1e6:.2f} MB in total "
             f"(median {z14['bytes'].median() / 1e3:.0f} KB)", 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 / "02-z14-tile-sizes.png")
plt.close(fig)

# Figure 3: attributes and precision
fig, ax = plt.subplots(figsize=(8, 2.8))
ax.barh(attr_df["variant"], attr_df["size_MB"], color=["#9aa3a8", "#1543b0", "#c21f3a", "#1b9e77", "#e8a33c"])
for y, v in enumerate(attr_df["size_MB"]):
    ax.text(v, y, f" {v:g} MB", va="center", fontsize=8)
ax.invert_yaxis()
ax.set_xlim(0, attr_df["size_MB"].max() * 1.2)
ax.set_title(f"{len(buildings):,} buildings at zoom 14: what attributes and options cost", loc="left")
ax.spines[["top", "right"]].set_visible(False)
fig.savefig(OUT / "03-attributes-precision.png")
plt.close(fig)

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