#!/usr/bin/env python3
"""Build Wisp's street TILES from an OSM PBF: one WBM1 file per grid tile, gzipped, plus a manifest.

    python3 build_tiles.py --pbf ohio-latest.osm.pbf --poly ohio.poly --out DIR
    python3 build_tiles.py --self-test

WHY TILES. Ruled by Phill 2026-09-15 (decisions.md, streets-everywhere entry): the chart's streets draw
WHEREVER the player is, from our own pre-built tiles, worldwide, with privacy designed in.
build_basemap.py built ONE bundled Columbus file through Overpass; that cannot scale past a metro.

THE GRID (world_contract.md 3.3g.8 is the authority; this docstring restates only what the code does):
  - A tile is TILE_DEG x TILE_DEG degrees, keyed by integer (row, col):
        row = floor((lat + 90) / TILE_DEG),  col = floor((lon + 180) / TILE_DEG)
  - Published path:  <GRID_ID>/<row>/<col>.wbm.gz   (gzip of an unmodified WBM1 file)
  - Geometry is clipped EXACTLY to the tile box, so two neighbouring tiles meet at a shared cut vertex
    and draw no ink twice. The header bbox is that box nudged one Float32 ulp outward (3.3g.3's clip
    obligation); multiples of 0.25 are exact in Float32, so the nudge is belt, not the seam mechanism.
    There is no buffer. Which tiles the client asks for is the client's business (H3Index.swift
    `tiles(around:)`); this file only guarantees that every tile it writes is complete for its own box.
  - Every tile the region polygon TOUCHES is written, INCLUDING EMPTY ONES (a 24-byte header), so a
    missing object means "never built", never "open ocean".
  - ⚖ A tile only PARTLY inside the polygon (a BORDER tile) IS written and is marked "border": true in
    the manifest (desk ruling 2026-09-15: dropping them was a build defect that left downtown Cincinnati
    dark). The client reads a tile's header box as covered ground, which is TRUE where the outside part
    is sea and FALSE where it is land in an extract not yet built. The complete fix builds each border
    tile from EVERY extract that touches it, de-duplicating ways by OSM ID; the flag names the tiles
    that merge has to rebuild.

THE FORMAT IS FROZEN (3.3g.3) and this file reuses build_basemap.py's packer, decoder and Float32
helpers rather than re-implementing them.

DATA (c) OpenStreetMap contributors, ODbL 1.0. Publishing these tiles publishes a derived database;
the share-alike offer and licence page are part of the same change (decisions.md 2026-09-15).
"""

from __future__ import annotations

import argparse
import gzip
import hashlib
import json
import math
import os
import shutil
import struct
import sys
import time

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from build_basemap import (  # noqa: E402  (the frozen format lives in one place)
    DEFAULT_STREET_KINDS, DEFAULT_WATERWAY_KINDS, KIND_STREET, KIND_WATER, MAGIC, MAX_N,
    decode, f32, f32_out_hi, f32_out_lo, split_long,
)

TILE_DEG = 0.25
# The client's inflate ceiling (BasemapTileStore.maxTileBytes). A tile above it would be refused on
# every phone forever, so the builder refuses to write one.
MAX_TILE_BYTES = 64 * 1024 * 1024
GRID_ID = "v1/q4"          # q4 = four tiles per degree. A different TILE_DEG is a different GRID_ID.
TILES_PER_DEG = 4
assert TILES_PER_DEG * TILE_DEG == 1.0


def tile_of(lat: float, lon: float) -> tuple[int, int]:
    """Longitude wraps (+180 is -180); latitude clamps at the poles, so +90 lands in the top row."""
    rows, cols = 180 * TILES_PER_DEG, 360 * TILES_PER_DEG
    row = min(rows - 1, max(0, math.floor((lat + 90.0) * TILES_PER_DEG)))
    return (row, math.floor((lon + 180.0) * TILES_PER_DEG) % cols)


def _grid_key(lat: float, lon: float) -> tuple[int, int]:
    """tile_of WITHOUT the longitude wrap, for building: a vertex at exactly +180 stays in the last column, whose
    box (179.75..180) contains it, instead of wrapping to column 0, whose box (-180..-179.75) does not and would
    fail the round-trip check after a continent-long scan (L1 F4, 2026-09-15)."""
    rows, cols = 180 * TILES_PER_DEG, 360 * TILES_PER_DEG
    return (min(rows - 1, max(0, math.floor((lat + 90.0) * TILES_PER_DEG))),
            min(cols - 1, max(0, math.floor((lon + 180.0) * TILES_PER_DEG))))


LICENSE_TEXT = """Wisp street tiles
Derived from OpenStreetMap. (c) OpenStreetMap contributors.

This database is made available under the Open Database License 1.0:
https://opendatacommons.org/licenses/odbl/1-0/
Any rights in individual contents of the database are licensed under the Database Contents License:
https://opendatacommons.org/licenses/dbcl/1-0/
OpenStreetMap copyright and licence: https://www.openstreetmap.org/copyright

Layout: v1/q4/<row>/<col>.wbm.gz, row = floor((lat + 90) * 4), col = floor((lon + 180) * 4).
manifest.json in this directory lists every tile. The format and the build method are described on
the licence page published with these tiles.
"""


def tile_box(row: int, col: int) -> tuple[float, float, float, float]:
    """(min_lat, min_lon, max_lat, max_lon), exact multiples of TILE_DEG."""
    return (row / TILES_PER_DEG - 90.0, col / TILES_PER_DEG - 180.0,
            (row + 1) / TILES_PER_DEG - 90.0, (col + 1) / TILES_PER_DEG - 180.0)


# ---------------------------------------------------------------------------
# Splitting a polyline at grid lines
# ---------------------------------------------------------------------------

def split_at_grid(pts):
    """Cut a lat/lon polyline wherever it crosses a tile edge.

    Returns [((row, col), [pts...]), ...]. A cut vertex is computed ONCE and appended to both sides,
    so the two pieces share it bit for bit. Each sub-segment belongs to the tile containing its
    MIDPOINT, which settles a segment lying exactly on a grid line without a special case.
    """
    out = []
    cur_key = None
    cur = []

    def emit():
        nonlocal cur
        if cur_key is not None and len(cur) >= 2:
            out.append((cur_key, cur))
        cur = []

    for i in range(len(pts) - 1):
        a, b = pts[i], pts[i + 1]
        dlat, dlon = b[0] - a[0], b[1] - a[1]
        ts = [0.0, 1.0]
        for delta, base in ((dlat, a[0] + 90.0), (dlon, a[1] + 180.0)):
            if delta == 0.0:
                continue
            lo, hi = sorted((base, base + delta))
            k = math.floor(lo * TILES_PER_DEG) + 1
            while k / TILES_PER_DEG < hi:
                t = (k / TILES_PER_DEG - base) / delta
                if 0.0 < t < 1.0:
                    ts.append(t)
                k += 1
        ts = sorted(set(ts))
        for j in range(len(ts) - 1):
            t0, t1 = ts[j], ts[j + 1]
            p0 = a if t0 == 0.0 else (a[0] + dlat * t0, a[1] + dlon * t0)
            p1 = b if t1 == 1.0 else (a[0] + dlat * t1, a[1] + dlon * t1)
            # snap the coordinate that sits on a grid line to that line exactly
            p0, p1 = _snap(p0), _snap(p1)
            tm = (t0 + t1) / 2
            key = _grid_key(a[0] + dlat * tm, a[1] + dlon * tm)
            if key != cur_key or not cur or cur[-1] != p0:
                emit()
                cur_key = key
                cur = [p0]
            cur.append(p1)
    emit()
    return out


def split_antimeridian(pts):
    """Break a polyline where a segment crosses the antimeridian (a longitude jump of more than 180°), ending one
    run exactly on the seam and starting the next on its other side.

    ⛔ Without this a segment from lon 179.9 to -179.9 reads as a 359.8° step WEST, and split_at_grid cuts it into
    a piece for every column on Earth (codex, L3 2026-09-15). OSM usually splits dateline features, but nothing in
    this builder may rely on that.
    """
    runs, cur = [], [pts[0]]
    for a, b in zip(pts, pts[1:]):
        if abs(b[1] - a[1]) > 180.0:
            east = a[1] > 0.0
            seam_a, seam_b = (180.0, -180.0) if east else (-180.0, 180.0)
            b_lon = b[1] + (360.0 if east else -360.0)   # b unwrapped onto a's side of the seam
            t = (seam_a - a[1]) / (b_lon - a[1])
            lat = a[0] + (b[0] - a[0]) * t
            if cur[-1] != (lat, seam_a):
                cur.append((lat, seam_a))
            if len(cur) >= 2:
                runs.append(cur)
            cur = [(lat, seam_b)]
        cur.append(b)
    if len(cur) >= 2:
        runs.append(cur)
    return runs


def _snap(p):
    lat, lon = p
    rl = round(lat * TILES_PER_DEG) / TILES_PER_DEG
    rn = round(lon * TILES_PER_DEG) / TILES_PER_DEG
    return (rl if abs(lat - rl) < 1e-9 else lat, rn if abs(lon - rn) < 1e-9 else lon)


# ---------------------------------------------------------------------------
# Region polygon (Geofabrik .poly)
# ---------------------------------------------------------------------------

def read_poly(path):
    """Geofabrik/osmosis .poly -> list of (is_hole, [(lon, lat), ...])."""
    rings, cur, hole = [], None, False
    with open(path) as fh:
        lines = [ln.strip() for ln in fh]
    for ln in lines[1:]:
        if not ln:
            continue
        if ln == "END":
            if cur is not None:
                rings.append((hole, cur))
                cur = None
            continue
        if cur is None:
            hole = ln.startswith("!")
            cur = []
            continue
        lon, lat = (float(x) for x in ln.split()[:2])
        cur.append((lon, lat))
    return rings


def _point_in_ring(lon, lat, ring):
    inside = False
    n = len(ring)
    for i in range(n):
        x1, y1 = ring[i]
        x2, y2 = ring[(i + 1) % n]
        if (y1 > lat) != (y2 > lat):
            x = x1 + (lat - y1) * (x2 - x1) / (y2 - y1)
            if x > lon:
                inside = not inside
    return inside


def _segments_cross(p1, p2, q1, q2):
    def orient(a, b, c):
        return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])
    d1, d2 = orient(q1, q2, p1), orient(q1, q2, p2)
    d3, d4 = orient(p1, p2, q1), orient(p1, p2, q2)
    return (d1 * d2 <= 0) and (d3 * d4 <= 0)


def tile_fully_inside(rings, box) -> bool:
    """All four corners inside the region AND no region edge crosses the tile boundary."""
    min_lat, min_lon, max_lat, max_lon = box
    corners = [(min_lon, min_lat), (max_lon, min_lat), (max_lon, max_lat), (min_lon, max_lat)]

    def inside(lon, lat):
        in_outer = any(_point_in_ring(lon, lat, r) for h, r in rings if not h)
        return in_outer and not any(_point_in_ring(lon, lat, r) for h, r in rings if h)

    if not all(inside(*c) for c in corners):
        return False
    # A ring lying wholly inside the tile (a small hole, an islet) crosses no tile edge, so the edge test below
    # cannot see it (L1 F5).
    if any(min_lon < lon < max_lon and min_lat < lat < max_lat for _, ring in rings for lon, lat in ring):
        return False
    edges = [(corners[i], corners[(i + 1) % 4]) for i in range(4)]
    for _, ring in rings:
        n = len(ring)
        for i in range(n):
            a, b = ring[i], ring[(i + 1) % n]
            if max(a[0], b[0]) < min_lon or min(a[0], b[0]) > max_lon:
                continue
            if max(a[1], b[1]) < min_lat or min(a[1], b[1]) > max_lat:
                continue
            if any(_segments_cross(a, b, e0, e1) for e0, e1 in edges):
                return False
    return True


def _candidate_tiles(rings):
    """Every tile whose box meets some ring's bounding box, built PER RING and with NO longitude wrap.

    ⛔ One bounding box over the whole region, keyed through tile_of, is wrong for a region that spans the
    antimeridian: Geofabrik's north-america.poly has rings at both -180 and +180 [measured 2026-09-15], and
    tile_of wraps +180 to column 0, so the old loop scanned a single column and wrote a sliver, not a continent.
    """
    keys = set()
    for _, ring in rings:
        lons = [p[0] for p in ring]
        lats = [p[1] for p in ring]
        r0, c0 = _grid_key(min(lats), min(lons))
        r1, c1 = _grid_key(max(lats), max(lons))
        keys.update((row, col) for row in range(r0, r1 + 1) for col in range(c0, c1 + 1))
    return sorted(keys)


def tiles_inside(rings):
    return [key for key in _candidate_tiles(rings) if tile_fully_inside(rings, tile_box(*key))]


def _in_region(rings, lon, lat) -> bool:
    in_outer = any(_point_in_ring(lon, lat, r) for h, r in rings if not h)
    return in_outer and not any(_point_in_ring(lon, lat, r) for h, r in rings if h)


def _segments_cross_properly(p1, p2, q1, q2):
    """Strict: the two segments pass through each other. Meeting at an endpoint, or running along each other,
    does not count."""
    def orient(a, b, c):
        return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])
    d1, d2 = orient(q1, q2, p1), orient(q1, q2, p2)
    d3, d4 = orient(p1, p2, q1), orient(p1, p2, q2)
    return (d1 * d2 < 0) and (d3 * d4 < 0)


def tile_touches(rings, box) -> bool:
    """The tile and the region share POSITIVE AREA: the tile's centre is inside the region, a region vertex lies
    strictly inside the tile, or a region edge properly crosses a tile edge. A tile wholly inside a hole touches
    nothing.

    ⛔ A tile that only shares an edge or a corner with the region does NOT touch it (codex + complement, L3
    2026-09-15). Written, such a tile would be an empty full-box header OUTSIDE the extract, suppressing the
    client's "no map data" state exactly where there is none.
    """
    min_lat, min_lon, max_lat, max_lon = box
    if _in_region(rings, (min_lon + max_lon) / 2, (min_lat + max_lat) / 2):
        return True
    corners = [(min_lon, min_lat), (max_lon, min_lat), (max_lon, max_lat), (min_lon, max_lat)]
    edges = [(corners[i], corners[(i + 1) % 4]) for i in range(4)]
    for _, ring in rings:
        n = len(ring)
        for i in range(n):
            a, b = ring[i], ring[(i + 1) % n]
            if min_lon < a[0] < max_lon and min_lat < a[1] < max_lat:
                return True
            if max(a[0], b[0]) < min_lon or min(a[0], b[0]) > max_lon:
                continue
            if max(a[1], b[1]) < min_lat or min(a[1], b[1]) > max_lat:
                continue
            if any(_segments_cross_properly(a, b, e0, e1) for e0, e1 in edges):
                return True
    return False


def tiles_touching(rings) -> dict:
    """Every tile the region touches -> True when it lies fully inside the region, False for a BORDER tile."""
    out = {}
    for key in _candidate_tiles(rings):
        box = tile_box(*key)
        if tile_touches(rings, box):
            out[key] = tile_fully_inside(rings, box)
    return out


# ---------------------------------------------------------------------------
# Per-tile body spool (bodies are appended to disk, never held for a continent)
# ---------------------------------------------------------------------------

class Spool:
    FLUSH_BYTES = 64 * 1024 * 1024

    def __init__(self, work_dir, wanted):
        self.dir = work_dir
        self.wanted = wanted
        self.buf: dict[tuple[int, int], bytearray] = {}
        self.count: dict[tuple[int, int], int] = {}
        self.vertices: dict[tuple[int, int], int] = {}
        self.pending = 0
        os.makedirs(work_dir, exist_ok=True)

    def add(self, key, kind, pts):
        if key not in self.wanted:
            return
        b = self.buf.setdefault(key, bytearray())
        for piece in split_long(pts):
            b += struct.pack("<BH", kind, len(piece))
            b += struct.pack(f"<{2 * len(piece)}f", *(v for p in piece for v in p))
            self.count[key] = self.count.get(key, 0) + 1
            self.vertices[key] = self.vertices.get(key, 0) + len(piece)
            self.pending += 3 + 8 * len(piece)
        if self.pending >= self.FLUSH_BYTES:
            self.flush()

    MIN_FREE_BYTES = 0

    def flush(self):
        free = shutil.disk_usage(self.dir).free
        if free < self.MIN_FREE_BYTES:
            raise SystemExit(f"✗ ABORT: {free / 1e9:.1f} GB free is under --min-free-gb "
                             f"{self.MIN_FREE_BYTES / 1e9:.0f}; spool left at {self.dir}")
        for key, b in self.buf.items():
            with open(self._path(key), "ab") as fh:
                fh.write(b)
        self.buf.clear()
        self.pending = 0

    def _path(self, key):
        return os.path.join(self.dir, f"{key[0]}_{key[1]}.body")

    def body(self, key):
        p = self._path(key)
        if not os.path.exists(p):
            return b""
        with open(p, "rb") as fh:
            return fh.read()


# ---------------------------------------------------------------------------
# Build
# ---------------------------------------------------------------------------

def build(args) -> int:
    import osmium  # flow venv; imported here so --self-test needs no pyosmium

    rings = read_poly(args.poly)
    touching = tiles_touching(rings)
    wanted = set(touching)
    if not wanted:
        print("✗ no tile touches the region polygon — nothing to build", file=sys.stderr)
        return 2
    border = sum(1 for fully_inside in touching.values() if not fully_inside)
    print(f"REGION {args.poly}: {len(wanted)} tiles touch it, {border} of them BORDER tiles "
          f"({GRID_ID}, {TILE_DEG}°)")

    keep_streets = {s for s in DEFAULT_STREET_KINDS.split(",") if s}
    keep_water = {s for s in DEFAULT_WATERWAY_KINDS.split(",") if s}

    work = os.path.join(args.out, ".spool")
    shutil.rmtree(work, ignore_errors=True)
    Spool.MIN_FREE_BYTES = int(args.min_free_gb * 1e9)
    spool = Spool(work, wanted)

    t0 = time.time()
    ways = kept = bad_loc = 0
    fp = osmium.FileProcessor(args.pbf)
    if args.index != "none":
        fp = fp.with_locations(args.index)
    # --index none: the ways already carry their node locations (`osmium add-locations-to-ways`), so no index
    # is built here; a way without them is counted in bad_loc and dropped, never drawn at 0,0.
    fp = fp.with_filter(osmium.filter.EntityFilter(osmium.osm.WAY))
    for w in fp:
        ways += 1
        tags = w.tags
        hw = tags.get("highway")
        if hw is not None:
            if hw not in keep_streets:
                continue
            kind = KIND_STREET
        else:
            ww = tags.get("waterway")
            if not ((ww in keep_water) or tags.get("natural") == "water"
                    or tags.get("landuse") in ("reservoir", "basin")):
                continue
            kind = KIND_WATER
        pts = []
        for n in w.nodes:
            if not n.location.valid():
                bad_loc += 1
                pts = None
                break
            pts.append((f32(n.lat), f32(n.lon)))
        if not pts or len(pts) < 2:
            continue
        kept += 1
        for run in split_antimeridian(pts):
            for key, piece in split_at_grid(run):
                spool.add(key, kind, piece)
        if ways % 2_000_000 == 0:
            print(f"  … {ways:,} ways scanned, {kept:,} kept, {time.time() - t0:.0f}s")
    spool.flush()
    print(f"SCAN {ways:,} ways, {kept:,} kept, {bad_loc:,} dropped for a missing node location, "
          f"{time.time() - t0:.0f}s")
    # ⛔ A PBF without node locations under --index none drops EVERY way and would still write a complete set of
    # valid empty tiles: a continent published as open ground (L1 F2). Refused before any tile is written.
    # ⚠ Judged on MISSING LOCATIONS, never on kept == 0: a region with no qualifying ways at all is legitimately
    # empty and must still be written, or open ground reads as "never built" (codex, L3 2026-09-15).
    if bad_loc * 1000 > kept + bad_loc:
        print(f"✗ ABORT before writing any tile: {kept:,} ways kept, {bad_loc:,} missing a node location "
              f"(limit 0.1%). Under --index none the PBF must come from `osmium add-locations-to-ways`.",
              file=sys.stderr)
        return 6

    manifest = {"grid": GRID_ID, "tile_degrees": TILE_DEG, "source": os.path.basename(args.pbf),
                "region": os.path.basename(args.poly), "built": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
                "licence": "ODbL 1.0 -- (c) OpenStreetMap contributors", "tiles": {}}
    total_raw = total_gz = empty = 0
    for i, key in enumerate(sorted(wanted)):
        if i % 1000 == 0 and shutil.disk_usage(args.out).free < Spool.MIN_FREE_BYTES:
            print(f"✗ ABORT at tile {i:,}: under --min-free-gb {args.min_free_gb:.0f} while writing (L1 F6)",
                  file=sys.stderr)
            return 7
        min_lat, min_lon, max_lat, max_lon = tile_box(*key)
        header = (f32_out_lo(min_lat), f32_out_lo(min_lon), f32_out_hi(max_lat), f32_out_hi(max_lon))
        blob = struct.pack("<4sffffI", MAGIC, *header, spool.count.get(key, 0)) + spool.body(key)
        if len(blob) > MAX_TILE_BYTES:
            print(f"✗ tile {key} is {len(blob):,} B, over the client's {MAX_TILE_BYTES:,} B inflate ceiling; "
                  f"it would be refused on every phone. Needs a finer grid there.", file=sys.stderr)
            return 4

        # ROUND TRIP through the independent decoder: every vertex inside the header box.
        (lo_lat, lo_lon, hi_lat, hi_lon), dec = decode(blob)
        outside = sum(1 for _, pts in dec for la, lo in pts
                      if not (lo_lat <= la <= hi_lat and lo_lon <= lo <= hi_lon))
        if outside or len(dec) != spool.count.get(key, 0):
            print(f"✗ tile {key}: {outside} vertices outside its header, {len(dec)} decoded "
                  f"vs {spool.count.get(key, 0)} written", file=sys.stderr)
            return 3

        gz = gzip.compress(blob, compresslevel=9, mtime=0)
        dest = os.path.join(args.out, GRID_ID, str(key[0]), f"{key[1]}.wbm.gz")
        os.makedirs(os.path.dirname(dest), exist_ok=True)
        with open(dest + ".part", "wb") as fh:
            fh.write(gz)
        os.replace(dest + ".part", dest)
        manifest["tiles"][f"{key[0]}/{key[1]}"] = {
            "sha256": hashlib.sha256(blob).hexdigest(), "bytes": len(blob), "gzip_bytes": len(gz),
            "polylines": spool.count.get(key, 0), "vertices": spool.vertices.get(key, 0)}
        if not touching[key]:
            manifest["tiles"][f"{key[0]}/{key[1]}"]["border"] = True
        total_raw += len(blob)
        total_gz += len(gz)
        empty += spool.count.get(key, 0) == 0

    shutil.rmtree(work, ignore_errors=True)
    with open(os.path.join(args.out, GRID_ID, "manifest.json"), "w") as fh:
        json.dump(manifest, fh, indent=1, sort_keys=True)
        fh.write("\n")
    # ODbL 4.2/4.3: the notice sits WITH the database, where a user fetching tiles would look (L2).
    with open(os.path.join(args.out, GRID_ID, "LICENSE.txt"), "w") as fh:
        fh.write(LICENSE_TEXT)
    biggest = max(manifest["tiles"].items(), key=lambda kv: kv[1]["bytes"])
    print(f"WROTE {len(wanted)} tiles ({empty} empty) under {os.path.join(args.out, GRID_ID)}")
    print(f"  WBM {total_raw:,} B · gzip {total_gz:,} B · largest {biggest[0]} "
          f"{biggest[1]['bytes']:,} B ({biggest[1]['gzip_bytes']:,} gzip)")
    print(f"  total {time.time() - t0:.0f}s")
    return 0


# ---------------------------------------------------------------------------
# Self test
# ---------------------------------------------------------------------------

def self_test() -> int:
    ok = True

    def check(name, good, detail=""):
        nonlocal ok
        ok = ok and good
        print(f"  {'ok  ' if good else 'FAIL'} {name} {detail}")

    print("SELF TEST -- grid keys:")
    check("Columbus downtown", tile_of(39.9612, -83.0007) == (519, 387), str(tile_of(39.9612, -83.0007)))
    check("a vertex on a grid line belongs to the tile above/right", tile_of(40.0, -83.0) == (520, 388))
    check("tile_box inverts tile_of", tile_box(519, 387) == (39.75, -83.25, 40.0, -83.0))
    check("south-west corner of the world", tile_of(-90.0, -180.0) == (0, 0))

    print("SELF TEST -- splitting at grid lines:")
    inside = split_at_grid([(39.8, -83.1), (39.9, -83.05)])
    check("a line inside one tile stays one piece", len(inside) == 1 and inside[0][0] == (519, 387))
    across = split_at_grid([(39.9, -83.1), (40.1, -83.1)])
    check("crossing a lat line makes two pieces", [k for k, _ in across] == [(519, 387), (520, 387)])
    check("the pieces share the cut vertex exactly", across[0][1][-1] == across[1][1][0] == (40.0, -83.1),
          str(across[0][1][-1]))
    diag = split_at_grid([(39.9, -83.1), (40.1, -82.9)])
    check("a diagonal through a corner reaches the far tile", diag[-1][0] == (520, 388),
          str([k for k, _ in diag]))
    joined = all(diag[i][1][-1] == diag[i + 1][1][0] for i in range(len(diag) - 1))
    check("every consecutive piece shares a vertex", joined)
    for key, pts in diag:
        lo_lat, lo_lon, hi_lat, hi_lon = tile_box(*key)
        check(f"piece in {key} lies inside its closed box",
              all(lo_lat <= a <= hi_lat and lo_lon <= b <= hi_lon for a, b in pts))
    along = split_at_grid([(40.0, -83.2), (40.0, -83.1)])
    check("a segment ON a grid line lands in exactly one tile", len(along) == 1, str(along))
    back = split_at_grid([(39.9, -83.1), (40.1, -83.1), (39.9, -83.1)])
    check("out and back re-enters the first tile as a new piece",
          [k for k, _ in back] == [(519, 387), (520, 387), (519, 387)])

    print("SELF TEST -- Float32-packed cut pieces stay inside their EXACT tile box:")
    import random
    rng = random.Random(7)
    escaped = 0
    for _ in range(2000):
        a = (39.0 + rng.random() * 2, -84.0 + rng.random() * 2)
        b = (a[0] + (rng.random() - 0.5) * 0.8, a[1] + (rng.random() - 0.5) * 0.8)
        for key, piece in split_at_grid([(f32(a[0]), f32(a[1])), (f32(b[0]), f32(b[1]))]):
            lo_lat, lo_lon, hi_lat, hi_lon = tile_box(*key)
            for la, lo in piece:
                la, lo = f32(la), f32(lo)
                if not (lo_lat <= la <= hi_lat and lo_lon <= lo <= hi_lon):
                    escaped += 1
    check("no packed vertex of 2,000 random segments leaves its exact box", escaped == 0, f"{escaped} escaped")
    check("f32_out_lo(39.75) < 39.75", f32_out_lo(39.75) < 39.75)
    print("SELF TEST -- the world's edges:")
    check("longitude +180 wraps to column 0", tile_of(10.0, 180.0) == (400, 0), str(tile_of(10.0, 180.0)))
    check("latitude +90 clamps into the top row", tile_of(90.0, 0.0) == (719, 720), str(tile_of(90.0, 0.0)))

    print("SELF TEST -- region containment:")
    square = [(False, [(-84.0, 39.0), (-82.0, 39.0), (-82.0, 41.0), (-84.0, 41.0)])]
    check("an interior tile is inside", tile_fully_inside(square, tile_box(519, 387)))
    check("a tile straddling the edge is not", not tile_fully_inside(square, tile_box(523, 383)),
          str(tile_box(523, 383)))
    notch = [(False, [(-84.0, 39.0), (-82.0, 39.0), (-82.0, 41.0), (-83.1, 41.0), (-83.1, 39.9),
                      (-83.15, 39.9), (-83.15, 41.0), (-84.0, 41.0)])]
    check("a thin notch cutting through a tile (corners all inside) is caught",
          not tile_fully_inside(notch, tile_box(519, 387)))
    found = set(tiles_inside(square))
    interior = {(r, c) for r in range(517, 523) for c in range(385, 391)}
    check("tiles_inside keeps every strictly interior tile", interior <= found, f"{len(found)} found")
    check("tiles_inside keeps nothing outside the square", all(
        -84.0 <= tile_box(*k)[1] and tile_box(*k)[3] <= -82.0 and 39.0 <= tile_box(*k)[0]
        and tile_box(*k)[2] <= 41.0 for k in found))

    print("SELF TEST -- region touching (border tiles):")
    # Edges OFF the quarter-degree grid, so a tile can truly straddle one rather than share it.
    offgrid = [(False, [(-84.1, 39.1), (-81.9, 39.1), (-81.9, 40.9), (-84.1, 40.9)])]
    check("an interior tile touches", tile_touches(offgrid, tile_box(519, 387)))
    check("a tile straddling the western edge touches", tile_touches(offgrid, tile_box(519, 383)),
          str(tile_box(519, 383)))
    check("a tile wholly outside does not", not tile_touches(offgrid, tile_box(519, 380)))
    speck = [(False, [(-83.2, 39.8), (-83.1, 39.8), (-83.1, 39.9)])]
    check("a region smaller than one tile, with no tile corner inside it, still touches",
          tile_touches(speck, tile_box(519, 387)))
    holed = [(False, [(-86.0, 37.0), (-80.0, 37.0), (-80.0, 43.0), (-86.0, 43.0)]),
             (True, [(-84.1, 39.1), (-81.9, 39.1), (-81.9, 40.9), (-84.1, 40.9)])]
    check("a tile wholly inside a hole does not", not tile_touches(holed, tile_box(519, 387)))
    t = tiles_touching(offgrid)
    fully = set(tiles_inside(offgrid))
    check("tiles_touching keeps every fully-inside tile, marked inside",
          bool(fully) and all(t.get(k) is True for k in fully), f"{len(fully)} fully inside")
    check("the straddling tile is kept and marked BORDER", t.get((519, 383)) is False)
    check("no kept tile lies wholly outside the square", all(
        tile_box(*k)[1] <= -81.9 and tile_box(*k)[3] >= -84.1 and tile_box(*k)[0] <= 40.9
        and tile_box(*k)[2] >= 39.1 for k in t), f"{len(t)} kept")

    print("SELF TEST -- a region spanning the antimeridian:")
    dateline = [(False, [(-180.0, 51.1), (-178.9, 51.1), (-178.9, 51.9), (-180.0, 51.9)]),
                (False, [(178.9, 51.1), (180.0, 51.1), (180.0, 51.9), (178.9, 51.9)])]
    cols = {c for _, c in tiles_touching(dateline)}
    check("both sides of the dateline are scanned, not one wrapped column",
          {0, 1, 2, 3, 1436, 1437, 1438, 1439} <= cols, str(sorted(cols)))
    check("tiles_inside reaches past the first column too", any(c > 0 for _, c in tiles_inside(dateline)))
    on_line = split_at_grid([(51.3, 180.0), (51.4, 180.0)])
    check("a segment lying ON lon +180 keys to the last column, not a wrapped column 0",
          [k for k, _ in on_line] == [(565, 1439)], str([k for k, _ in on_line]))
    check("...and every vertex of it lies inside that tile's box", all(
        tile_box(*k)[1] <= lo <= tile_box(*k)[3] and tile_box(*k)[0] <= la <= tile_box(*k)[2]
        for k, pts in on_line for la, lo in pts))
    approach = split_at_grid([(51.3, 179.9), (51.35, 180.0)])
    check("a segment ending at +180 stays in the last column", all(k[1] == 1439 for k, _ in approach),
          str([k for k, _ in approach]))

    print("SELF TEST -- a hole wholly inside one tile:")
    pocket = [(False, [(-86.0, 37.0), (-80.0, 37.0), (-80.0, 43.0), (-86.0, 43.0)]),
              (True, [(-83.2, 39.8), (-83.1, 39.8), (-83.1, 39.9)])]
    check("the tile holding the hole is not fully inside", not tile_fully_inside(pocket, tile_box(519, 387)))
    check("...so tiles_touching keeps it and marks it BORDER", tiles_touching(pocket).get((519, 387)) is False)
    check("a neighbouring tile away from the hole is still fully inside", tile_fully_inside(pocket, tile_box(519, 389)))

    print("SELF TEST -- a region whose edges lie ON grid lines (positive area only):")
    grid_sq = [(False, [(-84.0, 39.0), (-82.0, 39.0), (-82.0, 41.0), (-84.0, 41.0)])]
    expect = {(r, c) for r in range(516, 524) for c in range(384, 392)}
    got = set(tiles_touching(grid_sq))
    check("exactly the 64 tiles with positive-area overlap, none sharing only an edge or corner",
          got == expect, f"{len(got)} got; extra {sorted(got - expect)[:4]}; missing {sorted(expect - got)[:4]}")

    print("SELF TEST -- a way crossing the antimeridian:")
    west = split_antimeridian([(10.0, 179.9), (10.0, -179.9)])
    check("eastbound across the seam becomes two runs meeting at +180 / -180",
          west == [[(10.0, 179.9), (10.0, 180.0)], [(10.0, -180.0), (10.0, -179.9)]], str(west))
    east = split_antimeridian([(10.0, -179.9), (12.0, 179.9)])
    check("westbound across the seam splits too, with the seam latitude interpolated",
          len(east) == 2 and east[0][-1] == (11.0, -180.0) and east[1][0] == (11.0, 180.0), str(east))
    pieces = [k for run in west for k, _ in split_at_grid(run)]
    check("...and the pieces land only in the seam's two columns, never across the world",
          {c for _, c in pieces} == {0, 1439}, str(pieces))
    check("a way that never crosses is one run, untouched",
          split_antimeridian([(10.0, 170.0), (10.5, 179.0)]) == [[(10.0, 170.0), (10.5, 179.0)]])

    print("\n" + ("ALL SELF TESTS PASSED" if ok else "*** SELF TEST FAILED ***"))
    return 0 if ok else 5


def main() -> int:
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--pbf")
    ap.add_argument("--poly", help="Geofabrik .poly for the region; every tile it touches is written, partial ones flagged border")
    ap.add_argument("--out", help="output root; tiles land under OUT/" + GRID_ID)
    ap.add_argument("--index", default="flex_mem",
                    help="pyosmium node-location index: flex_mem for a state; sparse_file_array,FILE on an "
                         "osmium tags-filter'ed continent (never dense_*: sized by the max node ID, see "
                         "world_contract.md 3.3g.8)")
    ap.add_argument("--min-free-gb", type=float, default=45.0,
                    help="abort at a spool flush if the output disk has less free space than this")
    ap.add_argument("--self-test", action="store_true")
    args = ap.parse_args()
    if args.self_test:
        return self_test()
    if not (args.pbf and args.poly and args.out):
        ap.error("--pbf, --poly and --out are required")
    return build(args)


if __name__ == "__main__":
    sys.exit(main())