feat: add raw-grid block layout with bilinear lat/lon point sampling and JS parity

Author: philipredstone

FORMAT.md

@@ -366,6 +366,7 @@ Block flag bits:
 |---:|---|---|
 | `0` | `has_leaf_directories` | The block has leaf directories |
 | `1` | `has_dict` | The block carries a per-block zstd dictionary |
+| `2` | `raw_grid` | Block stores a native source-grid raster instead of a Hilbert tile pyramid; see "Raw-Grid Blocks" below |
 
 Other bits are reserved.
 
@@ -570,6 +571,103 @@ The payload is zstd-compressed vertical deltas. This codec is valid for `u8` and
 - Each following row stores `current - value_above`.
 - Arithmetic wraps modulo `2^(8 * dtype_bytes)`.
 
+### Raw-Grid Blocks
+
+Files encoded with `--no-tiles` skip the Web-Mercator pyramid. Each
+`(variable_id, time_id)` pair becomes one **raw-grid block** that stores the
+native source lat-lon grid chunked into fixed-size sub-tiles. Point queries by
+`(lat, lon)` are O(1) range requests; the file does not drive a slippy-map
+viewer without on-the-fly resampling.
+
+A raw-grid block reuses the 64-byte block header with `BlockFlagRawGrid`
+(`1 << 2`) set in `block_flags`. Field interpretation changes:
+
+- `root_directory_offset` and `root_directory_length` point at the compressed
+  raw-grid section (the equivalent of the tile directory).
+- `leaf_directories_offset` and `leaf_directories_length` are zero.
+- `tile_data_offset` and `tile_data_length` describe the concatenated chunk
+  payload region.
+- `num_addressed_tiles` and `num_directory_entries` hold the total chunk count
+  (`chunk_count_x * chunk_count_y`).
+- `num_tile_contents` holds the deduplicated chunk count.
+
+The block-table entry's `codec` is the dominant chunk codec ID for stats; each
+chunk payload still carries its own one-byte codec tag.
+
+#### Raw-Grid Section
+
+The raw-grid section is internally compressed (same compression as the file's
+catalogs and directories). After decompression:
+
+| Offset | Size | Type | Field |
+|---:|---:|---|---|
+| `0` | 1 | `u8` | `schema_version`, currently `1` |
+| `1` | 1 | `u8` | `chunk_size_log2`; chunk side in source pixels = `1 << value` |
+| `2` | 2 | bytes | Reserved, written as zero |
+| `4` | 4 | `u32` | `nx`, source grid width in pixels |
+| `8` | 4 | `u32` | `ny`, source grid height in pixels |
+| `16` | 8 | `f64` | `lat0`, latitude at source row `0` |
+| `24` | 8 | `f64` | `lon0`, longitude at source column `0` |
+| `32` | 8 | `f64` | `dy`, latitude step per row (may be negative) |
+| `40` | 8 | `f64` | `dx`, longitude step per column (may be negative) |
+| `48` | 8 | `f64` | `missing_value`, source NoData sentinel (may be NaN) |
+| `56` | 4 | `u32` | `chunk_count_x`, = `ceil(nx / chunk_size)` |
+| `60` | 4 | `u32` | `chunk_count_y`, = `ceil(ny / chunk_size)` |
+| `64` | varint × N | LEB128 `u64` each | `chunk_offsets`, one per chunk row-major |
+| ... | varint × N | LEB128 `u64` each | `chunk_lengths`, one per chunk row-major |
+
+`chunk_size_log2` must be in `[4, 12]`, allowing chunk sides of 16..4096
+source pixels. `N = chunk_count_x * chunk_count_y`. Chunk index is
+`cy * chunk_count_x + cx`, so chunks are stored in row-major order.
+
+`chunk_offsets[i]` is a byte offset relative to the block's `tile_data_offset`.
+`chunk_lengths[i]` is the chunk payload length. A `(0, 0)` pair signals
+**absent chunk**; decoders fill all pixels with NaN without fetching anything.
+
+#### Chunk Payloads
+
+Each chunk payload is a normal tile blob: one codec tag byte followed by the
+codec-specific stream. The current encoder uses `0x01` (constant) when all
+quantised values in the chunk are identical, otherwise `0x03`
+(bitshuffle + zstd). Lorenzo and delta codecs are not emitted for raw-grid
+chunks because edge chunks at the right/bottom border may be non-square.
+
+Chunk pixel count is `chunk_width(cx) * chunk_height(cy)` where:
+
+```text
+chunk_width(cx)  = min(chunk_size, nx - cx*chunk_size)
+chunk_height(cy) = min(chunk_size, ny - cy*chunk_size)
+```
+
+Pixels inside a chunk are row-major: `chunk[row*chunk_width + col]` is the
+quantised value at source coordinates `(cx*chunk_size + col, cy*chunk_size + row)`,
+which corresponds to lat/lon `(lat0 + (cy*chunk_size + row)*dy, lon0 + (cx*chunk_size + col)*dx)`.
+
+#### Point Sampling
+
+To compute the value at `(lat, lon)`:
+
+1. Compute source-grid coordinates `gx = (lon - lon0) / dx`,
+   `gy = (lat - lat0) / dy`. NaN or out-of-range inputs return NaN.
+2. Find the four neighbours `(x0, y0)`, `(x1, y0)`, `(x0, y1)`, `(x1, y1)` with
+   `x0 = floor(gx)`, `y0 = floor(gy)`, `x1 = x0+1`, `y1 = y0+1` clamped to
+   `[0, nx-1]` / `[0, ny-1]`.
+3. For each neighbour, locate the chunk `cx = x / chunk_size`,
+   `cy = y / chunk_size`, fetch and decode the chunk, then read the pixel.
+4. Bilinearly interpolate: `wx = gx - x0`, `wy = gy - y0`,
+   `v = ((1-wx)*v00 + wx*v10)*(1-wy) + ((1-wx)*v01 + wx*v11)*wy`.
+   If any neighbour is NaN, the result is NaN.
+
+For batched queries the reader unions the chunk indices touched by all points
+(including the 2×2 bilinear neighbourhood) and coalesces adjacent chunk byte
+ranges into shared HTTP-range requests, the same strategy as tile coalescing.
+
+A raw-grid file's header records `min_zoom = 0`, `max_zoom = 0`, and the
+authoritative signal that consumers should branch on is `BlockFlagRawGrid` on
+each block header. Mixing tiled and raw-grid blocks in the same file is
+permitted by the wire format but the current encoder does not produce mixed
+files.
+
 ### Lorenzo Zstd (`0x05`)
 
 The payload is zstd-compressed 2D Lorenzo residuals. This codec is valid for

README.md

@@ -82,6 +82,11 @@ wmtiles encode 'composite_wn_*-hd5' -o radar.wmt --max-zoom 7
 
 # CF-1.x / NetCDF4 file (regular lat-lon coords)
 wmtiles encode model.nc4 -o model.wmt
+
+# Weather-API mode: skip the Web-Mercator pyramid, store the source grid
+# chunked in source coords. Bilinear point queries via Sample / sample();
+# encodes ~50-100x faster, files shrink to ~GRIB×0.5..1.0.
+wmtiles encode forecast.grib2 -o api.wmt --no-tiles
 ```
 
 The input format is auto-detected by magic bytes (`GRIB` vs `\x89HDF`) with a
@@ -128,6 +133,23 @@ const px = await t2m.tile({ time: 12, z: 5, x: 16, y: 11 });
 // Float32Array(256*256), NaN where the encoder marked NoData
 ```
 
+For `--no-tiles` archives use the lat/lon sample API; the same range
+coalescing keeps a batch of points down to a single byte-range request when
+they fall in the same source-grid chunk neighbourhood:
+
+```ts
+const v = wmt.variable("2t_2m");
+const tempBerlin = await v.sample({ time: 0, lat: 52.52, lon: 13.40 });
+
+const cities = [
+  { lat: 52.52, lon: 13.40 },  // Berlin
+  { lat: 48.14, lon: 11.58 },  // Munich
+  { lat: 53.55, lon:  9.99 },  // Hamburg
+];
+const values = await v.samples({ time: 0, points: cities });
+// Float32Array(3) — NaN outside the source bbox.
+```
+
 For multi-tile fetches at the same `(variable, time)`, `tiles()`
 coalesces 9 viewport tiles into 1 to 2 range requests:
 
@@ -171,6 +193,18 @@ Read one tile:
 pixels, err := wmt.ReadTile("2t", 12, decode.Coord(5, 16, 11))
 ```
 
+For `--no-tiles` files use point sampling:
+
+```go
+v, err := wmt.Sample("2t", 12, 52.52, 13.40) // lat, lon
+// Float32; NaN outside the source bbox.
+
+values, err := wmt.Samples("2t", 12, []decode.SamplePoint{
+    {Lat: 52.52, Lon: 13.40},
+    {Lat: 48.14, Lon: 11.58},
+})
+```
+
 Read a viewport worth of tiles with range coalescing:
 
 ```go
@@ -400,6 +434,8 @@ wmtiles serve            <file.wmt> [--addr :8080]     bundled web viewer
 | `--tile-size-log2 N` | `8` (256 px) | tile pixel size, allowed `7..10` (128..1024) |
 | `--filter SHORTNAMES` | (none = all) | comma-separated shortNames to keep (GRIB shortName, ODIM quantity, or CF mapping) |
 | `--precision NAME=K,…` | shortName/unit lookup, then 10-bit auto-cap | quantisation precision overrides; `=0` forces full-range u16 |
+| `--no-tiles` | off | skip the Web-Mercator pyramid; store source-grid chunks for point-query (lat/lon) API use. Output is not viewable on a slippy map without on-the-fly tiling |
+| `--raw-chunk-size-log2 N` | `8` (256 px) | source-pixel side of one raw-grid chunk as log2 (4..12 → 16..4096). Only consulted with `--no-tiles` |
 
 ---
 

cmd/gen-testdata/main.go

@@ -39,6 +39,57 @@ func main() {
 	multistepPath := filepath.Join(outDir, "multistep.wmt")
 	makeMultistep(multistepPath, pixSize, refTime, fixedNow)
 	fmt.Printf("wrote %s\n", multistepPath)
+
+	rawGridPath := filepath.Join(outDir, "rawgrid.wmt")
+	makeRawGrid(rawGridPath, refTime, fixedNow)
+	fmt.Printf("wrote %s\n", rawGridPath)
+}
+
+// makeRawGrid writes a small --no-tiles file with a deterministic synthetic
+// grid. Used by the JS round-trip tests to lock in raw-grid sample values.
+func makeRawGrid(path string, refTime time.Time, now time.Time) {
+	const nx, ny = 64, 33
+	const lat0, lon0 = 30.0, 0.0
+	const dy, dx = 0.5, 0.5
+	values := make([]float32, nx*ny)
+	for y := 0; y < ny; y++ {
+		for x := 0; x < nx; x++ {
+			// pixel = lat + lon/1000, deterministic and easy to verify.
+			values[y*nx+x] = float32(lat0+float64(y)*dy) + float32(lon0+float64(x)*dx)/1000
+		}
+	}
+
+	opts := encoder.Options{
+		TilePixelSizeLog2:     7,
+		MinZoom:               0,
+		MaxZoom:               0,
+		ReferenceForecastTime: refTime,
+		TimeCatalog: format.TimeCatalog{
+			Regular: true, StartMs: refTime.UnixMilli(), IntervalMs: 0, Count: 1,
+		},
+		BBox: [4]float64{lon0, lat0, lon0 + float64(nx-1)*dx, lat0 + float64(ny-1)*dy},
+		Variables: []encoder.VariableSpec{
+			{Name: "temp", Unit: "K", Precision: 0.001},
+		},
+		CreationTime:        now,
+		SkipInternalWorkers: true,
+	}
+	enc, err := encoder.NewStreamingEncoder(opts, path)
+	if err != nil {
+		die(err)
+	}
+	if err := enc.EncodeRawGridBlock(encoder.RawGridSpec{
+		Variable: "temp", TimeStep: 0,
+		Nx: nx, Ny: ny,
+		Lat0: lat0, Lon0: lon0,
+		DY: dy, DX: dx,
+		Precision: 0.001,
+	}, values); err != nil {
+		die(err)
+	}
+	if err := enc.Finish(); err != nil {
+		die(err)
+	}
 }
 
 // 4 hourly steps × 2 variables. Pixel values encode (timeStep, variable) so a

cmd/wmtiles/cli.go

@@ -8,8 +8,6 @@ import (
 	"github.com/hstin-de/wmtiles/parser"
 )
 
-// Each field is one subcommand for kong; the Run() method on the subcommand
-// struct receives the parsed flags and drives the real work.
 type CLI struct {
 	Encode          encodeCmd          `cmd:"" help:"encode GRIB2 or HDF5 sources into a fresh .wmt"`
 	Extend          extendCmd          `cmd:"" help:"append blocks for new (variable, time) pairs to an existing .wmt"`
@@ -24,9 +22,11 @@ type CLI struct {
 type encodeCmd struct {
 	Output            string `short:"o" required:"" placeholder:"PATH" help:"output .wmt path"`
 	Format            string `enum:"auto,grib2,hdf5" default:"auto" help:"override input-format auto-detection"`
-	MinZoom           uint8  `default:"0" help:"minimum zoom level"`
-	MaxZoom           uint8  `default:"5" help:"maximum zoom level"`
+	MinZoom           uint8  `default:"0" help:"minimum zoom level (ignored when --no-tiles is set)"`
+	MaxZoom           uint8  `default:"5" help:"maximum zoom level (ignored when --no-tiles is set)"`
 	TileSizeLog2      uint8  `name:"tile-size-log2" default:"8" help:"tile size as log2 of pixel count (7..10 -> 128..1024)"`
+	NoTiles           bool   `name:"no-tiles" help:"skip the Web-Mercator pyramid; store source-grid chunks for point-query (lat/lon) API use. Output is not viewable on a slippy map without on-the-fly tiling"`
+	RawChunkSizeLog2  uint8  `name:"raw-chunk-size-log2" default:"8" help:"source-pixel side of one raw-grid chunk as log2 (4..12 -> 16..4096). Only consulted with --no-tiles"`
 	Filter            string `help:"comma-separated source variable shortNames to keep (default: all)"`
 	Precision         string `placeholder:"NAME=K,..." help:"per-variable quantisation precision overrides (default: lookup table + auto-cap)"`
 	DisableDeltaCodec bool   `help:"force bitshuffle-only encoding (faster, larger files)"`
@@ -44,9 +44,12 @@ func (c *encodeCmd) Validate() error {
 	if c.TileSizeLog2 < 7 || c.TileSizeLog2 > 10 {
 		return fmt.Errorf("tile-size-log2 must be 7..10, got %d", c.TileSizeLog2)
 	}
-	if c.MaxZoom < c.MinZoom {
+	if !c.NoTiles && c.MaxZoom < c.MinZoom {
 		return fmt.Errorf("max-zoom %d < min-zoom %d", c.MaxZoom, c.MinZoom)
 	}
+	if c.NoTiles && (c.RawChunkSizeLog2 < 4 || c.RawChunkSizeLog2 > 12) {
+		return fmt.Errorf("raw-chunk-size-log2 must be 4..12, got %d", c.RawChunkSizeLog2)
+	}
 	return nil
 }
 

cmd/wmtiles/encode_cmd.go

@@ -187,8 +187,13 @@ func runEncode(c *encodeCmd) error {
 
 	ui.Section("Settings")
 	ui.KV("source format", formatLabel)
-	ui.KVf("zoom range", "%d..%d", c.MinZoom, c.MaxZoom)
-	ui.KVf("tile size", "%d px", 1<<c.TileSizeLog2)
+	if c.NoTiles {
+		ui.KV("mode", "raw-grid (no tile pyramid)")
+		ui.KVf("raw chunk size", "%d px", 1<<c.RawChunkSizeLog2)
+	} else {
+		ui.KVf("zoom range", "%d..%d", c.MinZoom, c.MaxZoom)
+		ui.KVf("tile size", "%d px", 1<<c.TileSizeLog2)
+	}
 	if c.Filter == "" {
 		ui.KV("filter", "all source variables")
 	} else {
@@ -232,14 +237,16 @@ func runEncode(c *encodeCmd) error {
 	}
 
 	opts := encode.Options{
-		TileSize:          1 << c.TileSizeLog2,
-		MinZoom:           c.MinZoom,
-		MaxZoom:           c.MaxZoom,
-		Precision:         overrides,
-		Metadata:          map[string]any{"sourceFormat": formatLabel, "sourceCount": len(inputs)},
-		DisableDeltaCodec: c.DisableDeltaCodec,
-		ZstdLevel:         c.ZstdLevel,
-		EnableTileDict:    c.TileDict,
+		TileSize:             1 << c.TileSizeLog2,
+		MinZoom:              c.MinZoom,
+		MaxZoom:              c.MaxZoom,
+		NoTiles:              c.NoTiles,
+		RawGridChunkSizeLog2: c.RawChunkSizeLog2,
+		Precision:            overrides,
+		Metadata:             map[string]any{"sourceFormat": formatLabel, "sourceCount": len(inputs)},
+		DisableDeltaCodec:    c.DisableDeltaCodec,
+		ZstdLevel:            c.ZstdLevel,
+		EnableTileDict:       c.TileDict,
 		OnInputScanned: func(name string, records int) {
 			if scanPhase == nil {
 				scanStart.Store(time.Now().UnixNano())