Memory-map snapshots (MAPSNAP2) with a stored Hilbert R-tree

.gitignore

@@ -6,3 +6,4 @@
 .DS_Store
 .vscode
 *.snap
+!/tests/fixtures/*

Cargo.toml

@@ -8,9 +8,11 @@ license = "MIT"
 [dependencies]
 anyhow = "1"
 axum = { version = "0.8", default-features = false, features = ["http1", "tokio", "json", "query"] }
+bytemuck = { version = "1.25", features = ["derive", "extern_crate_alloc"] }
 clap = { version = "4.6", features = ["derive"] }
 earcutr = "0.5"
 flate2 = "1"
+memmap2 = "0.9.12"
 mimalloc = { version = "0.1", default-features = false }
 moka = { version = "0.12", features = ["future"] }
 osmpbf = "0.3.8"

compose.yaml

@@ -1,5 +1,6 @@
 # MAP=./nyc.osm.pbf docker compose up --build
-# The map lives in memory: NYC peaks near 3 GB loading, then about 1.7 GB plus cache.
+# Built from a .osm.pbf the map lives in memory (NYC peaks near 3 GB, then about
+# 1.7 GB plus cache). A snapshot (`maps build`) is mapped from disk instead.
 services:
   maps:
     build: .

examples/near.rs

@@ -8,7 +8,7 @@ fn main() {
     let p = project(lon, lat);
     let mpu = meters_per_unit(p[1]);
     for f in &map.features {
-        if !(f.kind.is_extrusion() || f.kind.is_point()) {
+        if !(f.kind().is_extrusion() || f.kind().is_point()) {
             continue;
         }
         let cx = (f.bbox[0] + f.bbox[2]) as f64 / 2.0 + map.origin[0];
@@ -23,6 +23,6 @@ fn main() {
             (f.bbox[3] - f.bbox[1]) as f64 * mpu,
         );
         let (clon, clat) = unproject([cx, cy]);
-        println!("{:?} flags={} base={} top={} roof={:?} rh={} vis={} {w:.0}x{h:.0}m @{clon:.5},{clat:.5}", f.kind, f.flags, f.base, f.height, det.roof, det.roof_height, f.vis_zoom);
+        println!("{:?} flags={} base={} top={} roof={:?} rh={} vis={} {w:.0}x{h:.0}m @{clon:.5},{clat:.5}", f.kind(), f.flags, f.base, f.height, det.roof, det.roof_height, f.vis_zoom);
     }
 }

examples/towers.rs

@@ -4,10 +4,10 @@ fn main() {
     let map = maps::map::Map::build(maps::ingest::read(std::path::Path::new(&path)).unwrap());
     let mut seen = std::collections::BTreeMap::new();
     for f in &map.features {
-        if f.kind != maps::classify::Kind::Tower {
+        if f.kind() != maps::classify::Kind::Tower {
             continue;
         }
-        let p = map.ring(f.ring_start)[0];
+        let p = map.first_point(f);
         let (lon, lat) =
             maps::geo::unproject([p[0] as f64 + map.origin[0], p[1] as f64 + map.origin[1]]);
         let e = seen

readme.md

@@ -18,8 +18,8 @@ pyramid offline.
 ## Features
 
 **WebGL client, like Google Maps**
-- `maps serve` keeps the map in memory and builds geometry tiles on
-  request: styled, triangulated ground, road ribbons, buildings and model
+- `maps serve` maps the map from a snapshot (or builds it from an
+  extract) and builds geometry tiles on request: styled, triangulated ground, road ribbons, buildings and model
   instances (`src/vtile.rs`).
 - The browser draws every frame itself with WebGL2 and a perspective
   camera, so pan, zoom, rotate and tilt are smooth at 60 fps.
@@ -97,9 +97,12 @@ cargo build --release
 # Feature counts and bounds
 ./target/release/maps info nyc.osm.pbf
 
-# Build once into a snapshot; every command loads it in seconds
+# Build once into a snapshot; every command maps it from disk at once
 ./target/release/maps build us-northeast.osm.pbf -o northeast.snap
 ./target/release/maps serve northeast.snap
+
+# Convert an older MAPSNAP1 snapshot without rebuilding
+./target/release/maps convert old.snap northeast.snap
 ```
 
 In the viewer, drag to pan and scroll to zoom. Right-drag or Ctrl+drag
@@ -158,28 +161,58 @@ Endpoints:
 | `/models/{version}/{theme}.bin` | Model templates for instances |
 | `/tiles/{version}/{style}/{bearing}/{pitch}/{z}/{x}/{y}[@2x].png` | Tiles. Immutable and cacheable forever, because `version` changes with the data. |
 | `/healthz` | Liveness probe (`ok`); `maps check` probes it from inside the image |
-| `/metrics` | Prometheus counters: requests, cache hits, renders, render time, geometry tiles and bytes sent, cache size |
+| `/metrics` | Prometheus counters: requests, cache hits, renders, render time, geometry tiles and bytes sent, and the size of the raster and geometry tile caches |
 
 The server shuts down gracefully on SIGTERM. Pushes to `main` publish
 `ghcr.io/kensac/maps` (`latest` and `sha-<short>`).
 
 ### Large extracts
 
-Building the map from an extract briefly needs about twice the memory it
-serves with. `maps build` does that once and writes a snapshot, which every
-command accepts in place of the `.osm.pbf`:
+Building the map from an extract needs about twice the memory the map
+itself takes. `maps build` does that once and writes a snapshot, which every
+command accepts in place of the `.osm.pbf`. A snapshot is memory-mapped and
+used in place, spatial index included, so a server is up at once and holds
+only the pages tiles touch, which the kernel can drop again:
 
 | US Northeast (1.8 GB PBF, 25 M features) | From PBF | From snapshot |
 | --- | ---: | ---: |
-| Load time (M1 Pro) | 42 s | 3.5 s |
-| Peak memory | 12 GB | 5.7 GB |
-
-Ingest keeps the build lean: node coordinates are written once into a
-packed 32-bit fixed-point array (about 1 cm), node references become
-32-bit indices so the IDs can be freed, raw ways are released as they are
-built, and the feature sort is in place.
-
-The snapshot is 3.8 GB. Serving it uses about 5.7 GB plus the tile cache.
+| Ready to serve | 72 s | 0.04 s |
+| Peak memory | 10 GB | 3 MB, plus 200 MB of file pages read by the checks |
+| Anonymous memory after 180 tiles in five cities (Linux) | | 135 MB, plus 1.5 GB of file pages |
+
+The snapshot is 4.3 GB, of which the index is 0.5 GB (21 bytes a feature).
+A tile that touches pages not yet read waits for the disk: in that test
+first renders took about twice as long as repeats.
+
+**Format.** `MAPSNAP2` (`src/snapshot.rs`): a fixed header (bounds, origin,
+background, tallest structure, and the offset and length of each section),
+then 64-byte-aligned sections: points, ring starts, features (56-byte
+records), roof and facade details (24-byte records), elevations, building
+groups, per-kind feature counts (so `/v1/kinds` needs no scan) and one packed
+Hilbert R-tree per zoom bucket (`src/index.rs`; 16 children a node, `f32`
+boxes, `u32` ids). Opening checks the header, every section's bounds and
+alignment, that ring starts are in order and end at the last point, and
+that the index has exactly the shape the builder writes, with every id in
+range. That reads the ring starts and index ids, about 200 MB for US
+Northeast (40 ms from the page cache). The 1.4 GB of feature records are not
+scanned: their kind and roof bytes go through checked accessors, and their
+ring, detail, elevation and group indices are bounds-checked where they are
+used, so a damaged feature draws nothing. A render that panics anyway fails
+its own request, not the server.
+
+**Older snapshots.** `MAPSNAP1` files still load, copied onto the heap with
+the index rebuilt (5.3 GB peak, 4.7 GB after loading, for US Northeast).
+`maps convert old.snap new.snap` converts one in seconds without the
+extract, copying the arrays and building only the index (1.3 GB peak for US
+Northeast).
+
+**Replacing a snapshot.** A running server reads the file in place, so
+changing it under the server crashes it (SIGBUS) or serves garbage. Never
+`cp` over a snapshot or write into it. Write the new file next to it as
+`<name>.partial`, then `mv` it over the old one on the same filesystem (a
+rename leaves the old file intact for anyone still mapping it), then
+restart. `maps build` and `maps convert` already write `<output>.partial`
+and rename it, syncing both the file and its directory.
 
 ## Performance
 
@@ -189,7 +222,7 @@ Measured on an Apple M1 Pro (10 cores), release build.
 | --- | ---: | ---: | ---: |
 | JFK airport | 2.3 MB | 59 k | 0.12 s |
 | Manhattan + Brooklyn | 16 MB | 361 k | 0.34 s |
-| New York metro | 184 MB | 3.9 M | 4.5 s, then 1.7 GB resident |
+| New York metro | 184 MB | 3.9 M | 4.5 s, then 1.7 GB resident; from a snapshot (590 MB) at once, with a 70 MB footprint after 300 tile requests |
 
 Serving the NYC extract (Midtown, 512 px retina tiles at bearing 29° and
 pitch 45°, 16 concurrent clients):
@@ -226,8 +259,8 @@ It classifies tags into physical kinds with real or typical dimensions
 - Turns coastlines into land.
 - Drops building outlines that are drawn through their parts, and lifts
   rooftop objects onto their buildings.
-- Sorts features into draw order and indexes them in one R-tree per zoom
-  bucket.
+- Sorts features into draw order and indexes them in one packed Hilbert
+  R-tree per zoom bucket (`src/index.rs`), stored in the snapshot.
 
 **Rendering** (`src/render/`) uses an orthographic camera with bearing and
 pitch. The ground is foreshortened by cos(pitch) and heights rise by

src/api.rs

@@ -198,10 +198,10 @@ impl Api {
 
     fn geometry(&self, f: &Feature) -> Value {
         let rings: Vec<&[[f32; 2]]> = self.map.rings(f).collect();
-        if f.kind.is_point() {
+        if f.kind().is_point() {
             return json!({ "type": "Point", "coordinates": self.lon_lat(rings[0][0]) });
         }
-        if !f.kind.is_area() {
+        if !f.kind().is_area() {
             let elev = self.map.elevations(f);
             let mut offset = 0;
             let lines: Vec<Value> = rings
@@ -255,8 +255,8 @@ impl Api {
 
     fn properties(&self, f: &Feature) -> Value {
         let mut p = json!({
-            "kind": kind_name(f.kind),
-            "geometry_type": geometry_type(f.kind),
+            "kind": kind_name(f.kind()),
+            "geometry_type": geometry_type(f.kind()),
         });
         let o = p.as_object_mut().unwrap();
         if f.layer != 0 {
@@ -272,13 +272,13 @@ impl Api {
                 o.insert(name.into(), json!(true));
             }
         }
-        if f.kind.is_extrusion() || f.kind.is_point() {
+        if f.kind().is_extrusion() || f.kind().is_point() {
             o.insert("height_m".into(), json!(round(f.height as f64, 2)));
             if f.base > 0.0 {
                 o.insert("min_height_m".into(), json!(round(f.base as f64, 2)));
             }
         }
-        if f.kind.is_extrusion() {
+        if f.kind().is_extrusion() {
             let d = self.map.detail(f);
             o.insert("roof_shape".into(), json!(roof_name(d.roof)));
             if d.levels > 0.0 {
@@ -416,7 +416,7 @@ async fn features(State(api): State<Shared>, Query(q): Query<FeaturesQuery>) ->
         .query(bbox_area(&api, bbox))
         .into_iter()
         .filter(|&i| {
-            let k = map.features[i as usize].kind;
+            let k = map.features[i as usize].kind();
             k != Kind::Land && kinds.as_ref().is_none_or(|ks| ks.contains(&k))
         })
         .collect();
@@ -441,7 +441,7 @@ async fn features(State(api): State<Shared>, Query(q): Query<FeaturesQuery>) ->
 
 async fn feature(State(api): State<Shared>, Path(id): Path<u32>) -> Response {
     match api.map.features.get(id as usize) {
-        Some(f) if f.kind != Kind::Land => ok(api.feature_json(id), 300),
+        Some(f) if f.kind() != Kind::Land => ok(api.feature_json(id), 300),
         _ => error(
             StatusCode::NOT_FOUND,
             "NOT_FOUND",
@@ -471,10 +471,10 @@ async fn lookup(State(api): State<Shared>, Query(q): Query<LookupQuery>) -> Resp
     let mut hits: Vec<(f32, u32)> = Vec::new();
     for i in api.query([p[0] - r, p[1] - r, p[0] + r, p[1] + r]) {
         let f = &map.features[i as usize];
-        if f.kind == Kind::Land {
+        if f.kind() == Kind::Land {
             continue;
         }
-        let hit = if f.kind.is_area() {
+        let hit = if f.kind().is_area() {
             let pt = [p[0] as f64, p[1] as f64];
             let inside = map
                 .rings(f)
@@ -538,7 +538,7 @@ async fn tallest(State(api): State<Shared>, Query(q): Query<TallestQuery>) -> Re
     let mut ids: Vec<u32> = api
         .query(bbox_area(&api, bbox))
         .into_iter()
-        .filter(|&i| map.features[i as usize].kind.is_extrusion())
+        .filter(|&i| map.features[i as usize].kind().is_extrusion())
         .collect();
     ids.sort_by(|&a, &b| {
         let (fa, fb) = (&map.features[a as usize], &map.features[b as usize]);
@@ -635,8 +635,9 @@ async fn static_map(State(api): State<Shared>, Query(q): Query<StaticQuery>) ->
             return;
         }
         let t = Instant::now();
-        let renderer = Renderer::new(map, theme.theme(), buildings);
-        let png = encode_png(&renderer.render(&vp));
+        let png = crate::server::contained(|| {
+            encode_png(&Renderer::new(map, theme.theme(), buildings).render(&vp))
+        });
         let _ = tx.send(png.map(|b| (b, t.elapsed())));
     });
     match rx.await {

src/classify.rs

@@ -375,6 +375,15 @@ pub enum RoofShape {
     Dome,
 }
 
+impl RoofShape {
+    pub fn from_u8(v: u8) -> Option<RoofShape> {
+        use RoofShape::*;
+        [Flat, Gabled, Hipped, Pyramidal, Skillion, Dome]
+            .get(v as usize)
+            .copied()
+    }
+}
+
 /// Architectural detail of an extruded solid.
 #[derive(Clone, Copy, Debug, Default, PartialEq)]
 pub struct Detail {
@@ -627,6 +636,15 @@ impl Group {
 }
 
 impl Kind {
+    /// Number of kinds; discriminants run `0..COUNT`.
+    pub const COUNT: usize = Kind::NavLight as usize + 1;
+
+    /// The kind with discriminant `v`, if there is one.
+    pub fn from_u8(v: u8) -> Option<Kind> {
+        // SAFETY: Kind is repr(u8) with contiguous discriminants 0..COUNT.
+        ((v as usize) < Kind::COUNT).then(|| unsafe { std::mem::transmute::<u8, Kind>(v) })
+    }
+
     /// Polygon geometry: ground areas and extruded solids.
     pub fn is_area(self) -> bool {
         self <= Kind::Tomb