Add out-of-core OpenStreetMap example (examples/osm)

examples/osm/.gitignore

@@ -0,0 +1,10 @@
+# Rust build artifacts for the bundled native parser/sorter
+osmium-rs/target/
+
+# Local OSM data + derived caches (can be hundreds of GB — never commit)
+*.pbf
+*.f64
+*.f32
+*.idx
+*.xyzones
+osm-data/

examples/osm/README.md

@@ -0,0 +1,99 @@
+# Every OpenStreetMap node, out-of-core
+
+Render **all ~10.7 billion OpenStreetMap planet nodes** as an interactive
+density scatter in `xy`, straight off disk. This is the end-to-end proof of the
+design dossier's §27 *"mmap (native)"* canonical-store row and the §2 *"100M+ /
+out-of-core — interactive via viewport tiling, bounded RAM"* target, at true
+planet scale: **172 GB of lon/lat lives on disk, resident memory stays
+screen-bounded** (the density pyramid), never data-bounded.
+
+The dataset used here: `planet-latest.osm.pbf` decoded to
+**10,742,674,832 nodes** → two f64 columns (`osm_lon.f64`, `osm_lat.f64`,
+86 GB each).
+
+## What's in here
+
+| File | Role |
+|---|---|
+| `osmium-rs/` | Native Rust crate: `osm-nodes` (PBF → f64 columns, ~250× pyosmium) and `osm-sort` (builds the Tier-3 spatial index). |
+| `ingest.py` | Parse a planet `.pbf` into f64 columns and render the out-of-core density scatter (with timings). |
+| `viewer.py` | Interactive browser viewer — serves `xy`'s real WebGL client over plain HTTP; pan/zoom re-aggregate every viewport from disk. |
+| `_pbf_split.py` | Standalone PBF blob-framing scanner (byte-range enumeration for parallel decode); reference/legacy helper. |
+
+## Prerequisites
+
+- **Rust** (for the native parser/sorter): `cd examples/osm/osmium-rs && cargo build --release`
+- **The render client**, built once from the repo root (it is not committed — see #214): `npm ci && node js/build.mjs`
+- **A planet PBF**: download `planet-latest.osm.pbf` from a mirror (~80 GB). *(Any `.osm.pbf` extract works — a city/country extract is a fast way to try the pipeline.)*
+- Disk: ~260 GB free for the full planet (172 GB columns + 84 GB spatial index).
+
+All commands below are run **from the repo root**.
+
+## 1. Parse: PBF → canonical f64 columns
+
+```bash
+python examples/osm/ingest.py \
+    --pbf /path/to/planet-latest.osm.pbf \
+    --out /path/to/osm-data
+```
+
+`osm-nodes` decodes DenseNodes (delta + zigzag varint) across all CPU cores
+with a zlib-ng inflate backend, writing `osm_lon.f64` / `osm_lat.f64` straight
+to disk. `ingest.py` then builds the figure and prints the memory report and a
+few zoom timings — proving `canonical_bytes` (RAM-resident) stays **0** while
+`canonical_mapped_bytes` (disk) holds all 172 GB.
+
+Re-run against columns already on disk with `--reuse` (skips parsing).
+
+## 2. (Optional but recommended) Build the deep-zoom spatial index
+
+The density pyramid alone answers zoomed-out views crisply, but its finest
+level over the whole planet is ~2.4 km/cell — so city/street zoom would only
+upsample (blurry). `osm-sort` spatially sorts every point into a grid-bucketed
+layout so a viewport reads only its in-window points — deep zoom then gets
+**sharper and cheaper the further you go**.
+
+```bash
+examples/osm/osmium-rs/target/release/osm-sort \
+    /path/to/osm-data/osm_lon.f64 \
+    /path/to/osm-data/osm_lat.f64 \
+    /path/to/osm-data/osm_spatial \
+    --grid 8192 --partitions 512
+```
+
+Produces `osm_spatial.lon.f32` / `osm_spatial.lat.f32` (sorted by cell) and
+`osm_spatial.idx` (header + cumulative cell offsets). It's an external counting
+sort — sequential I/O only, peak RAM ≈ one partition. `viewer.py` attaches it
+automatically if `osm_spatial.idx` sits next to the columns.
+
+## 3. View
+
+```bash
+python examples/osm/viewer.py --out /path/to/osm-data --port 8777
+# open http://localhost:8777/
+```
+
+Drag to pan, scroll to zoom. Every viewport is re-aggregated through the same
+`channel.handle_message` kernel protocol the notebook widget uses. The tier
+served rides each update as `binning`:
+
+- `pyramid-L<l>` / `-upsampled` — zoomed-out aggregate from the pre-binned pyramid.
+- `spatial-exact` (`filter: nearest`) — deep zoom re-binned exactly from the
+  in-window points at full screen resolution (crisp streets, no interpolation blur).
+- `spatial-points` — deep enough that the window fits the direct budget: the
+  real points ship as individual marks.
+
+## Rough timings (10.7B nodes, NVMe SSD, 16 cores)
+
+| Stage | Cost | Notes |
+|---|---|---|
+| Parse (`osm-nodes`) | minutes | ~250× the pyosmium binding; multithreaded + zlib-ng |
+| Spatial sort (`osm-sort`) | ~8 min (~21 M pts/s) | 3-pass parallel external counting sort |
+| Figure build (zone maps) | **~51 s first time, then instant** | one 172 GB scan; result cached to a `.xyzones` sidecar next to each column, so every later run (e.g. viewer restart) loads it in milliseconds |
+| Pan / zoom | single-digit ms | screen-bounded per viewport, independent of N |
+
+The `.xyzones` zone-map cache is a general `xy` out-of-core feature (see
+`python/xy/columns.py`): the first `figure()` over a memmapped column persists
+its per-chunk statistics, and reopening validates the sidecar against the
+file's size/mtime before reusing it — so the 51 s domain scan is paid once, not
+on every restart.

examples/osm/_pbf_split.py

@@ -0,0 +1,93 @@
+"""Minimal PBF *framing* scanner — enumerate blob byte-ranges for parallel
+parsing. This reads only the fixed 4-byte length prefixes and the small
+`BlobHeader` protobuf (type + datasize) to walk from one blob to the next; it
+never touches the compressed OSM payload. pyosmium still performs the actual
+(trusted) coordinate decode on each byte range via `osmium.io.FileBuffer`.
+
+A .osm.pbf file is a sequence of:
+    [4-byte big-endian BlobHeader length][BlobHeader protobuf][Blob protobuf]
+The first blob is an `OSMHeader`; the rest are `OSMData`. A byte range made of
+`OSMHeader` bytes followed by any run of whole `OSMData` blobs is itself a valid
+PBF stream, so it can be parsed independently in another process.
+"""
+
+from __future__ import annotations
+
+import struct
+from dataclasses import dataclass
+
+
+def _read_varint(buf: memoryview, pos: int) -> tuple[int, int]:
+    """Decode a protobuf base-128 varint at ``pos``; return (value, new_pos)."""
+    result = 0
+    shift = 0
+    while True:
+        b = buf[pos]
+        pos += 1
+        result |= (b & 0x7F) << shift
+        if not (b & 0x80):
+            return result, pos
+        shift += 7
+
+
+def _blob_header_fields(bh: memoryview) -> tuple[str, int]:
+    """Extract (type, datasize) from a BlobHeader protobuf (fields 1 and 3)."""
+    pos = 0
+    btype = ""
+    datasize = -1
+    n = len(bh)
+    while pos < n:
+        tag, pos = _read_varint(bh, pos)
+        field, wire = tag >> 3, tag & 0x7
+        if wire == 2:  # length-delimited
+            ln, pos = _read_varint(bh, pos)
+            if field == 1:  # type (string)
+                btype = bytes(bh[pos : pos + ln]).decode("ascii")
+            pos += ln
+        elif wire == 0:  # varint
+            val, pos = _read_varint(bh, pos)
+            if field == 3:  # datasize
+                datasize = val
+        else:  # 32/64-bit — not used by BlobHeader, skip defensively
+            pos += 4 if wire == 5 else 8
+    if datasize < 0:
+        raise ValueError("BlobHeader missing datasize")
+    return btype, datasize
+
+
+@dataclass
+class Blob:
+    start: int  # byte offset of the 4-byte length prefix
+    end: int  # byte offset just past the blob body
+    btype: str
+
+
+def scan_blobs(mm: memoryview) -> list[Blob]:
+    """Walk every blob in the mapped file, returning their byte ranges."""
+    blobs: list[Blob] = []
+    pos = 0
+    total = len(mm)
+    while pos + 4 <= total:
+        (hlen,) = struct.unpack_from(">I", mm, pos)
+        hstart = pos + 4
+        btype, datasize = _blob_header_fields(mm[hstart : hstart + hlen])
+        end = hstart + hlen + datasize
+        if end > total:
+            break  # truncated tail (e.g. partial download) — stop cleanly
+        blobs.append(Blob(pos, end, btype))
+        pos = end
+    return blobs
+
+
+def partition(blobs: list[Blob], n_parts: int) -> list[tuple[int, int]]:
+    """Split the OSMData blobs into ``n_parts`` contiguous (start,end) byte
+    ranges. The OSMHeader blob (blobs[0]) is excluded — callers prepend it."""
+    data = [b for b in blobs if b.btype == "OSMData"]
+    if not data:
+        return []
+    per = -(-len(data) // n_parts)  # ceil
+    parts: list[tuple[int, int]] = []
+    for i in range(0, len(data), per):
+        group = data[i : i + per]
+        parts.append((group[0].start, group[-1].end))
+    return parts

examples/osm/ingest.py

@@ -0,0 +1,147 @@
+"""Ingest **every OpenStreetMap node** (~9 billion, from `planet-latest.osm.pbf`)
+into xy's out-of-core canonical store and render it as a density scatter.
+
+This is the end-to-end proof of the §27 "mmap (native)" canonical-store row and
+the §2 "100M+ / out-of-core — interactive via viewport tiling, bounded RAM"
+target, at true planet scale. Peak resident memory stays screen-bounded (the
+density pyramid), not data-bounded (172 GB of lon/lat on disk).
+
+Pipeline:
+  planet.pbf ──osm-nodes (native)──▶ osm_lon.f64 / osm_lat.f64 (disk f64)
+            ──▶ xy.scatter(density=True) ──▶ build_payload / density_view
+
+Usage (from the repo root):
+  python examples/osm/ingest.py \
+      --pbf /path/to/planet-latest.osm.pbf \
+      --out /path/to/osm-data \
+      [--reuse]   # skip parsing; reuse osm_lon.f64/osm_lat.f64 on disk
+
+The native parser (examples/osm/osmium-rs) is built with:
+  cd examples/osm/osmium-rs && cargo build --release
+"""
+
+from __future__ import annotations
+
+import argparse
+import os
+import resource
+import subprocess
+import sys
+import time
+
+import numpy as np
+
+# examples/osm/ → repo root is two levels up; the package lives in python/.
+_HERE = os.path.dirname(__file__)
+_REPO_ROOT = os.path.abspath(os.path.join(_HERE, "..", ".."))
+sys.path.insert(0, os.path.join(_REPO_ROOT, "python"))
+import xy  # noqa: E402
+from xy._ooc import open_f64  # noqa: E402
+
+
+def _rss_gb() -> float:
+    return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / (1024 * 1024)  # GB on Linux
+
+
+def parse_planet(pbf: str, out_dir: str) -> tuple[np.memmap, np.memmap]:
+    """Decode all node coordinates via the native `osm-nodes` parser
+    (`osmium-rs/`, ~250x faster than the Python binding), then memmap the
+    resulting f64 columns. The parser writes canonical columns straight to
+    disk; xy never holds the data in RAM."""
+    px = os.path.join(out_dir, "osm_lon.f64")
+    py = os.path.join(out_dir, "osm_lat.f64")
+    binary = os.path.join(_HERE, "osmium-rs", "target", "release", "osm-nodes")
+    if not os.path.exists(binary):
+        raise SystemExit(
+            f"native parser not built: {binary}\n"
+            "build it with:  cd examples/osm/osmium-rs && cargo build --release"
+        )
+    t0 = time.perf_counter()
+    subprocess.run([binary, pbf, px, py], check=True)
+    print(f"native parse wall-clock: {time.perf_counter() - t0:.0f}s", flush=True)
+    xcol, ycol = open_f64(px), open_f64(py)
+    print(
+        f"canonical on disk: {xcol.nbytes + ycol.nbytes:,} bytes across 2 columns "
+        f"({len(xcol):,} nodes)",
+        flush=True,
+    )
+    return xcol, ycol
+
+
+def bench(xcol: np.memmap, ycol: np.memmap) -> None:
+    n = len(xcol)
+    print(f"\n=== rendering {n:,} points (out-of-core density scatter) ===", flush=True)
+
+    t0 = time.perf_counter()
+    fig = xy.chart(xy.scatter(x=xcol, y=ycol, density=True)).figure()
+    print(
+        f"figure build (ingest + zone maps, one disk scan): {time.perf_counter() - t0:.0f}s",
+        flush=True,
+    )
+
+    rep = fig.store.memory_report()
+    print(f"  canonical_bytes (RAM-resident): {rep['canonical_bytes']:,}")
+    print(f"  canonical_mapped_bytes (disk):  {rep['canonical_mapped_bytes']:,}")
+    assert rep["canonical_bytes"] == 0, "canonical data leaked into RAM"
+
+    t0 = time.perf_counter()
+    spec, blob = fig.build_payload(2048)
+    tr = spec["traces"][0]
+    print(
+        f"first paint: {time.perf_counter() - t0:.0f}s | wire blob = {len(blob):,} B "
+        f"| tier={tr.get('tier')}",
+        flush=True,
+    )
+
+    # Zoom into a continent-scale window; should serve from the pyramid.
+    for label, (x0, x1, y0, y1) in {
+        "world": (-180.0, 180.0, -85.0, 85.0),
+        "europe": (-12.0, 32.0, 35.0, 60.0),
+        "london": (-0.5, 0.3, 51.3, 51.7),
+    }.items():
+        t0 = time.perf_counter()
+        s2, bufs2 = fig.density_view(0, x0, x1, y0, y1, 1000, 800)
+        tr2 = s2["traces"][0] if s2.get("traces") else {}
+        nbytes = sum(len(b) for b in bufs2)
+        print(
+            f"  zoom [{label}]: {(time.perf_counter() - t0) * 1e3:.0f} ms "
+            f"| {len(bufs2)} bufs / {nbytes:,} B "
+            f"| tier={tr2.get('tier')} | binning={tr2.get('binning')} | visible={tr2.get('visible')}",
+            flush=True,
+        )
+
+    full = fig.memory_report()
+    print(
+        f"\nresident_array_bytes: {full['resident_array_bytes']:,} "
+        f"(canonical RAM {full['canonical_bytes']:,} + channels {full['channel_bytes']:,} "
+        f"+ pyramid {full['pyramid_bytes']:,})"
+    )
+    print(f"canonical on disk (mapped): {full['canonical_mapped_bytes']:,} bytes")
+    print(f"peak process RSS: {_rss_gb():.2f} GB")
+    print("\nOK — every OSM node rendered with screen-bounded resident memory.")
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--pbf", help="planet-latest.osm.pbf (omit with --reuse)")
+    ap.add_argument("--out", required=True, help="output dir for the f64 columns")
+    ap.add_argument(
+        "--reuse",
+        action="store_true",
+        help="skip parsing; reuse osm_lon.f64/osm_lat.f64 already on disk",
+    )
+    args = ap.parse_args()
+    if not args.reuse and not args.pbf:
+        ap.error("--pbf is required unless --reuse is given")
+
+    if args.reuse:
+        xcol = open_f64(os.path.join(args.out, "osm_lon.f64"))
+        ycol = open_f64(os.path.join(args.out, "osm_lat.f64"))
+        print(f"reusing {len(xcol):,} nodes from disk", flush=True)
+    else:
+        xcol, ycol = parse_planet(args.pbf, args.out)
+    bench(xcol, ycol)
+
+
+if __name__ == "__main__":
+    main()

examples/osm/osmium-rs/.gitignore

@@ -0,0 +1 @@
+/target