Metadata-Version: 2.4
Name: amphi-ncdata
Version: 0.6.4
Summary: Python client library for loading and sampling Amphitrite's NetCDF datasets.
Author: Amphitrite S.A.
License: Proprietary
Classifier: Programming Language :: Python :: 3
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: Microsoft :: Windows
Requires-Python: >=3.10
Description-Content-Type: text/markdown
Requires-Dist: netCDF4>=1.6.5
Requires-Dist: numpy>=1.24
Provides-Extra: dev
Requires-Dist: build>=1.2; extra == "dev"
Requires-Dist: pytest>=8.0; extra == "dev"
Dynamic: requires-python

# amphi-ncdata

`amphi-ncdata` is a Python client library for loading and sampling Amphi NetCDF datasets.

It exposes these main entry points:

- `AmphiNcData`: the easiest high-level API for clients using Amphi-style dataset folders.
- `MultiAmphiNcData`: an `AmphiNcData`-style env spanning **several bounding boxes**, routing each query to the owning box in O(1).
- `NcData`: a lower-level API for loading arbitrary NetCDF files or basepaths.

This README is for clients using an already built wheel. Repository maintainers should use `README_DEV.md`.

All timestamps in the public API are Unix timestamps in seconds.

## Install From The Wheel

Ask for the wheel that matches:

- your operating system
- your Python version
- your machine architecture

Install it with:

```bash
python3 -m pip install /path/to/amphi_ncdata-0.1.0-<platform>.whl
```

Quick import check:

```bash
python3 -c "import amphi_ncdata; print(amphi_ncdata.__version__)"
```

If installation succeeds but import fails with a missing shared-library error, ask us for a wheel built for your environment or install the required NetCDF runtime libraries for your platform.

## Quick Start

The fastest way to use the library is with `AmphiNcData`.

Example script:

```python
from amphi_ncdata import AmphiNcData

env = AmphiNcData(
    {
        "currentsBasepaths": ["/data/data_api/FORECAST-CURRENT/"],
        "windBasepaths": ["/data/data_api/FORECAST/wind/"],
        "waveBasepaths": ["/data/data_api/FORECAST/wave/"],
        "climWindBasepaths": ["/data/data_api/HINDCAST-3Y/"],
        "climWaveBasepaths": ["/data/data_api/HINDCAST-3Y-WAVE/"],
    },
    start=1781478000,
    end=1781650800,
    south=35.0,
    west=-73.0,
    north=43.0,
    east=-30.0,
)

lat, lon = 40.0, -50.0
timestamp = 1781524800

currents = env.currents()
wind = env.wind()
waves = env.waves()

currents.set_interpolation_method("linear")
wind.set_interpolation_method("linear")
waves.set_interpolation_method("linear")

u_current, v_current = currents.at(lat, lon, timestamp)
u_wind, v_wind = wind.at(lat, lon, timestamp)
wave_height, wave_period, wave_direction = waves.at(lat, lon, timestamp)

print("Currents:", u_current, v_current)
print("Wind:", u_wind, v_wind)
print("Waves:", wave_height, wave_period, wave_direction)
```

A ready-to-copy version also lives in `examples/client_quickstart.py`.

## Main Concepts

- `AmphiNcData` is the convenience layer for Amphi folder layouts, masks, and tides.
- `NcData` is the generic layer when you want to choose the files or variables yourself.
- `load()` registers which variables you want available.
- `set_interpolation_method(method, skip_nan=...)` sets the dataset interpolation mode once. `skip_nan=True` lets `"linear"` interpolate from the corners that hold data instead of returning NaN whenever one of the four is a coastal no-data cell.
- `set_no_data_value(value, when=...)` gives the product's no-data mask a value instead of NaN, **per variable**, either `when="before_interpolation"` or `when="after_interpolation"`.
- `MultiAmphiNcData(..., inclusive=..., snap_to_grid=..., cell_size=...)` tune how boxes claim the routing grid; all three are optional and every default matches earlier behaviour.
- `set_out_of_window_behaviour(behaviour, value, dim=...)` sets what a query outside the loaded window returns, **per dimension** (`dim="time"` / `"space"` / `"both"`): `"nearest"`, `"nan"`, `"value"`, or `"fallback"` (another `NcData`). Defaults: time→`nearest`, space→`nan`.
- `sample()` returns one variable.
- `at()` returns all loaded variables for the dataset.
- `at()` / `sample()` are fast for sequential, point-by-point access (e.g. routing); `sample_batch()` / `at_batch()` are best when you already have many independent points (roughly a few dozen or more) to sample at once. Pass `out=` (a buffer from `empty_batch()`) to sample in place with no per-call allocation.
- `env.wind()` / `env.waves()` **auto-fall-back to climatology past the forecast horizon** (when climatology basepaths are set); `env.climWind()` / `env.climWave()` are the underlying climatology datasets.

## Using `AmphiNcData`

`AmphiNcData` is the recommended API for most client code.

### Using an `input.json`-style config

```python
from amphi_ncdata import AmphiNcData
import json

with open("input.json", "r", encoding="utf-8") as f:
    cfg = json.load(f)

env = AmphiNcData(
    cfg,
    start=1781478000,
    end=1781650800,
    south=35.0,
    west=-73.0,
    north=43.0,
    east=-30.0,
)
```

### Using a flat config

```python
env = AmphiNcData(
    {
        "currentsBasepaths": ["/data/data_api/FORECAST-CURRENT/"],
        "tidesBasepaths": ["/data/data_api/FORECAST-TIDE/"],
        "windBasepaths": ["/data/data_api/FORECAST/wind/"],
        "waveBasepaths": ["/data/data_api/FORECAST/wave/"],
        "climWindBasepaths": ["/data/data_api/HINDCAST-3Y/"],
        "climWaveBasepaths": ["/data/data_api/HINDCAST-3Y-WAVE/"],
        "bathymetryPath": "/data/data_save/ETOPO/bathy.nc",
        "tidesMaskPath": "/data/data_api/FORECAST-TIDE/masks_tide.nc",
        "nogoPath": "/data/data_user/nogo.nc",
        "coastalZonePath": "/data/data_save/coastal_zone_mask.nc",
    },
    start=1781478000,
    end=1781650800,
    south=35.0,
    west=-73.0,
    north=43.0,
    east=-30.0,
)
```

### High-level helpers

```python
currents = env.currents()
wind = env.wind()
waves = env.waves()
tidal = env.tidal_currents()
bathy = env.bathymetry()
nogo = env.nogo()
coastal = env.coastal_zone()
tides_mask = env.tides_mask()
```

Return shapes:

- `env.currents().at(...)` -> `(ucos, vcos)`
- `env.tidal_currents().at(...)` -> `(utotal, vtotal)`
- `env.wind().at(...)` -> `(u10, v10)`
- `env.waves().at(...)` -> `(swh, mwp, mwd)`
- `env.bathymetry().at(...)` -> `bathymetry`
- `env.nogo().at(...)` -> `mask`
- `env.coastal_zone().at(...)` -> `mask`
- `env.tides_mask().at(...)` -> `mask_TIDE`

### Tidal currents: the fused `tidal_currents()`

Tides only matter near coasts (the global tide mask is ~9% active), so
`env.tidal_currents()` does **not** load the whole bounding box. With a
`tidesMaskPath` configured it loads **eager 1-degree masked tiles**: every
1-degree tile that the mask marks active is pulled into RAM up front, in the raw
NetCDF dtype, so the routing hot path never touches NetCDF I/O. Open-ocean and
land tiles allocate nothing.

`tidal_currents()` is **fused**: `at()` returns the hourly tidal current where
the mask is active (coastal) and the open-ocean current where it is not — so a
single call gives the right water velocity everywhere, no manual branch:

```python
from datetime import datetime, timezone
from amphi_ncdata import AmphiNcData

start = int(datetime(2026, 6, 23, tzinfo=timezone.utc).timestamp())
env = AmphiNcData(
    {
        "tidesBasepaths":    ["/data/data_api/FORECAST-TIDE/"],
        # The mask is what turns tides into eager 1-degree masked tiles:
        "tidesMaskPath":     "/data/data_api/FORECAST-TIDE/masks_tide.nc",
        # "tidesTileSizeDegrees": 1.0,   # optional; 1.0 is the default
    },
    start=start, end=start + 24 * 3600,
    south=48.0, west=-6.0, north=52.5, east=4.5,    # English Channel
)

tides = env.tidal_currents()        # eager 1° masked tiles + a daily open-ocean step
u, v = tides.at(lat, lon, timestamp)   # tide near the coast, current offshore — one call
```

**How the open ocean is served.** Where the tide product is not tidal (mask 0),
its values are constant over the day and equal the surface currents (verified:
`p99 = 0`, mean diff `3.7e-5` m/s over thousands of open-ocean points). So
`tidal_currents()` loads **one** daily tide step over the full bbox as the
open-ocean field — self-contained (no separate currents product needed), same
variable names, and roughly half the RAM of loading the multi-step currents
product. Coastal cells are overridden by the hourly tide tiles.

You normally never call the gate yourself, but it is exposed for inspection and
for `fused=False`:

```python
tides.mask(lat, lon)                 # compiled gate (~0.2 us): is this point tidal?
tides.mask_batch(lats, lons)         # vectorised
```

`mask()` is the nearest high-res `mask_TIDE` cell (active = non-zero); it is
spatial, so `timestamp` is accepted for symmetry but ignored. A masked-in cell
can still be land / fill (NaN) — `mask()` selects the data *source*, not its
presence.

Knobs: `tidesTileSizeDegrees` (default `1.0`) sets the tile size;
`TidalCurrents(..., eager=False)` keeps tiles lazy with an LRU cache (smaller
resident footprint, per-call I/O); `TidalCurrents(..., fused=False)` makes `at()`
tide-only again (NaN in the open ocean) if you want to gate manually.

#### What if you omit the tide paths?

| configuration | `env.tidal_currents()` behaviour |
|---|---|
| `tidesBasepaths` **and** `tidesMaskPath` set | **Eager 1° masked tiles + fused open ocean** (the recommended setup). `at()` is tide-where-coastal, current-where-not; `mask()` reflects the high-res tide mask. |
| `tidesBasepaths` set, **`tidesMaskPath` omitted** | Falls back to the **full-bbox tidal view** (lazy, no tiling, no fusion). `at()` samples tides everywhere in the box, and `mask()` returns `True` for any in-box point (no mask to gate on). Memory is the full bbox × all tide steps — fine for a small box, large for a wide one. |
| **`tidesBasepaths` omitted** | `env.tidal_currents()` raises `AmphiDiscoveryError` ("No valid NetCDF files were found …"). Only call it when tides are configured; `currents`, `wind`, and `waves` are independent and still work — `env.currents()` stays the currents-only, memory-light path. |

### Out-of-window behaviour (per dimension, with climatology fallback)

A query can fall outside the loaded window in **time** (past/before the forecast)
or in **space** (outside the bounding box). Each is handled independently:

```python
wind.set_out_of_window_behaviour(behaviour, value=0.0, dim="both")
```

- `dim`: `"time"`, `"space"`, or `"both"` (default).
- `behaviour`:
  - `"nearest"` — clamp to the nearest timestep / edge cell.
  - `"nan"` — return NaN.
  - `"value"` — return the constant `value` (a float).
  - `"fallback"` (alias `"data"`) — sample `value` instead, which must be another
    `NcData` that covers the query (e.g. a climatology). The fallback dataset
    carries its own out-of-window behaviour.
- **Defaults: time → `nearest`, space → `nan`.**

Precedence when out of window: a `fallback` dimension delegates the whole sample
to its dataset (time is checked before space); otherwise a `nan`/`value` on space
applies before time; `nearest` simply clamps that dimension. This is independent
of NaN coming from the data itself (land / fill), which is always NaN.
(`set_out_of_window_behaviour("nan")` / `("value", -999)` still set both
dimensions, for backward compatibility.)

#### Climatology fallback (automatic for wind / waves)

`env.wind()` and `env.waves()` **automatically fall back to climatology past the
forecast horizon** when climatology basepaths are configured — equivalent to:

```python
wind.set_out_of_window_behaviour("fallback", env.climWind(), dim="time")
```

So a long voyage just keeps sampling, with no manual gap handling:

```python
wind = env.wind()                      # forecast + climatology time-fallback
u, v = wind.at(40.0, -50.0, ts)        # forecast in-window, climatology past the horizon
```

The climatology (`env.climWind()` -> `(u10, v10)`, `env.climWave()` ->
`(swh, mwp, mwd)`) is still its own dataset and is configured as the safety net:
`nearest` in both time and space and `linear` interpolation, so it always returns
a value. Pass no climatology basepaths, or call
`wind.set_out_of_window_behaviour("nan", dim="time")`, to opt out (e.g. to detect
the gap yourself).

This holds for a window lying **wholly** past the horizon too, not just one that
straddles it. With no forecast file for the window at all the dataset is *empty*
rather than absent, and an empty dataset is out of window everywhere — so every
query goes to climatology:

```python
# Window 15 days out; wind forecast only reaches ~today+14, so nothing covers it.
env = AmphiNcData(config, start=far_future, end=far_future + 86400, **bbox)
wind = env.wind()
wind.is_empty            # True  -- no file backs it
wind.warnings            # WARN_NO_FILES
wind.at(lat, lon, ts)    # the climatology value, not an exception
wind.time_range          # None -- it holds no timestep; guard before indexing
```

### Missing data is reported, not raised: `warnings`

A forecast archive routinely has gaps, so *ordinary* discovery problems do not
raise. They set a bitmask on the dataset and the accessor still returns something
usable:

| code | meaning |
| --- | --- |
| `WARN_NONE` (0) | nothing was wrong |
| `WARN_NO_FILES` (1) | no file at all for the window; the dataset is empty, so queries return the out-of-window value or a wired fallback |
| `WARN_PARTIAL_COVERAGE` (2) | some of the requested days resolved, some did not — including a day whose file has no time axis and was skipped |
| `WARN_LAST_AVAILABLE` (4) | nothing covered the window, so the newest file on disk was served — **the values are stale** |
| `WARN_CLIMATOLOGY_MISSING` (8) | a climatology basepath was configured but its `MM.nc` is absent, so there is no safety net |
| `WARN_HETEROGENEOUS_FILES` (16) | the window's files change format partway (grid, resolution or packing); each format is sampled on its own grid, so values are right, but the resolution changes at the switch — see [below](#windows-that-span-a-change-of-file-format) |

It is a bitmask, not a sequential code, because these co-occur — a window can be
both partially covered and served from a stale run:

```python
from amphi_ncdata import WARN_LAST_AVAILABLE, WARN_NO_FILES, warning_names

currents = env.currents()
if currents.warnings & WARN_LAST_AVAILABLE:
    print("stale:", currents.warning_messages)

env.warnings                # OR across every dataset built so far; 0 == clean
warning_names(env.warnings) # ('no_files', 'last_available')
env.warning_messages        # ('wind: no NetCDF files were found for 2026-08-22..2026-08-23', ...)
```

`MultiAmphiNcData` aggregates the same way, per box: `routed.warnings` ORs its
members and `routed.warning_messages` prefixes each line with `box <i>`.

#### Turning them back into exceptions: `on_warning`

```python
AmphiNcData(config, ..., on_warning="raise")                      # every code raises
AmphiNcData(config, ..., on_warning={"no_files": "raise"})         # just this one
AmphiNcData(config, ..., on_warning={"default": "raise",
                                     "partial_coverage": "silent"})
```

Genuine caller errors are **never** downgraded and always raise, whatever
`on_warning` says: an unknown variable name, a bounding box that does not overlap
the grid, an unconfigured mask or climatology path. The distinction is whether the
data is missing (a fact about the archive) or the request is wrong.

### Windows that span a change of file format

An archive folder can switch product from one day to the next: a wave hindcast
built from WAVERYS (0.2°, packed int16, latitude ascending, longitude −180..180,
3-hourly) up to one day and from ERA5 (0.5°, float32, latitude descending,
longitude 0..360, hourly) after it, or a forecast whose resolution is upgraded. A
window that straddles such a switch is **split into one part per format**, and each
part is sampled on its own native grid, exactly as that product would be on its
own. Nothing changes in how you query it:

```python
from amphi_ncdata import WARN_HETEROGENEOUS_FILES

waves = env.waves()                       # window 2026-04-30 .. 2026-05-01
waves.warnings & WARN_HETEROGENEOUS_FILES  # set: the files change format
waves.warning_messages
# ('the files change format at .../2026/05/20260501.nc (grid 899x1800 (0.2 deg, lat -89.8..89.8,
#   lon -180..179.8) -> 361x720 (0.5 deg, lat 90..-90, lon 0..359.5); mwd int16 x0.01 +180
#   fill -32767 -> float32 fill nan); each format is sampled on its own grid',)

for part in waves.segments:               # one per format, in time order
    print(part.time_range, [round(step, 3) for step in part.spatial_resolution])
# (1777507200, 1777582800) [0.2, 0.2]
# (1777593600, 1777593600) [0.5, 0.5]

wave_height, wave_period, wave_direction = waves.at(lat, lon, timestamp)
```

- **Inside a part**, every value is exactly what that part's files return on their own.
- **Across the switch** — between one part's last timestep and the next part's first,
  21:00 on the 0.2° grid and 00:00 on the 0.5° one above — each side is interpolated
  on its own grid and the two are blended in time, with the same rules as any two
  timesteps (`skip_nan`, `set_no_data_value`, out-of-window behaviour, climatology
  fallback).
- **No resampling and no decoding.** Each part keeps its native grid and dtype, so
  memory is the sum of the parts, and the compiled sampler serves every query —
  `MultiAmphiNcData` keeps its compiled dispatch over such a box too.
- `latitudes`, `longitudes` and `spatial_resolution` describe the **finest** part;
  read `segments` for each part's own. A single-format dataset is its own single
  segment.
- A point at the very edge of the box can fall inside one part's grid and outside
  the other's. Across the switch, where both are sampled, it counts as outside the
  box and follows the space out-of-window behaviour.
- A variable that one of the formats lacks raises `AmphiVariableError` on `load()`,
  naming the file, as it does for any file missing it.

A file with no time axis at all cannot be placed in the window. It is skipped and
reported with `WARN_PARTIAL_COVERAGE`, instead of failing the whole window.

### `load_defaults()`

```python
env.load_defaults(
    include_tides=True,
    include_bathymetry=True,
    include_nogo=True,
    include_coastal_zone=True,
)
```

This preloads:

- currents
- wind
- waves
- optionally tidal currents and tides mask (tides eager-load their active 1° tiles up
  front — fast for a route corridor, heavier for a continental box; see the gate section)
- optionally bathymetry
- optionally no-go mask
- optionally coastal zone mask

### Overriding default variable names

If your files use different names:

```python
custom_currents = env.currents(u_var="u_current", v_var="v_current")
custom_wind = env.wind(u_var="u_wind", v_var="v_wind")
custom_waves = env.waves(height_var="hs", period_var="tp", direction_var="dir")
```

## Multiple bounding boxes: `MultiAmphiNcData`

Use `MultiAmphiNcData` when one env should cover **several** bounding boxes /
time windows at once — e.g. tiling a large region, or layering a smaller box on
top of a bigger one. All boxes share the same `config`; each carries its own
`start, end, south, west, north, east`. Every accessor mirrors `AmphiNcData`
(`currents()`, `wind()`, `waves()`, `tidal_currents()`, `bathymetry()`, …) but
returns a routed dataset whose `at`/`at_batch`/`sample`/`sample_batch` dispatch
each point to the box that owns it.

```python
from amphi_ncdata import MultiAmphiNcData   # or AmphiNcData.from_bounding_boxes(...)

env = MultiAmphiNcData(
    {"currentsBasepaths": ["/data/data_api/FORECAST-CURRENT/"]},   # shared by all boxes
    boxes=[
        dict(start=1781478000, end=1781650800, south=35.0, west=-73.0, north=43.0, east=-30.0),
        dict(start=1781478000, end=1781650800, south=43.0, west=-30.0, north=60.0, east=10.0),
    ],
    # out_of_box_value=float("nan"),   # what a point in no box returns (default NaN)
)

currents = env.currents()
currents.set_interpolation_method("linear")

u, v = currents.at(40.0, -50.0, 1781524800)         # routed to the owning box in C, O(1)
uv   = currents.at_batch(lats, lons, timestamps)    # routed + sampled in C, one call

buf = currents.empty_batch(len(lats))               # reusable (n, 2) float32 buffer
currents.at_batch(lats, lons, timestamps, out=buf)  # fill in place, returns buf
```

How dispatch works:

- A dense `(lat, lon, hour)` lookup grid (1° × 1° × 1 h) stores, per cell, the
  index of the winning box. A query computes three integer indices and reads one
  cell — the heavy bilinear + time interpolation still runs in each box's
  compiled sampler.
- **Both `.at()` and `.at_batch()` run compiled**, matching a plain dataset:
  `at`/`sample` route each point in C, and `at_batch`/`sample_batch` route **and**
  sample every point in C (no Python bucketing). They accept the same `out=`
  buffer and `empty_batch()` helper (see *Reusing an output buffer with `out=`*
  under **Using `NcData`** below), so batching is fast from small sizes up.
- **Overlap → first listed box wins.** When boxes overlap in both space and
  time, the earliest box in `boxes` that covers the point is chosen (so you can
  layer a high-res box on top of a coarse one by listing it first).
- A point in **no** box returns `out_of_box_value` (default `NaN`).
- **Cells cut by a box edge are resolved per query** against the real box
  bounds, so routing is exact wherever the edges fall — including a point in a
  shared edge cell, which goes to the box that actually holds it rather than to
  whichever claimed the cell. `inclusive` and `snap_to_grid` below remove that
  per-query step for layouts that do not need it.
- If every box shares the same `[start, end]`, the time axis collapses and the
  grid stays tiny (~130 KB); otherwise it grows with the total time span
  (~22 MB for a 7-day window).

### `inclusive`: say which edges each box owns

Box bounds are closed by default, so `[0°, 10°]` and `[10°, 20°]` **both** contain
latitude 10.0 — they overlap on their shared edge, and the cell holding it has to
be resolved per query.

`inclusive` lets you declare, per box, which of the six edges are inclusive:

```python
MultiAmphiNcData(config, boxes, inclusive="half_open")   # default for every box
```

`"half_open"` keeps `south`/`west` inclusive and makes `north`/`east` exclusive,
so abutting boxes **partition** space instead of overlapping. Nothing is left to
resolve, and the box edges stay exactly where you put them.

What changes, for boxes `[0°, 10°]` and `[10°, 20°]` listed in that order:

| query latitude | `closed` (default) | `half_open` |
| --- | --- | --- |
| 9.999 | box 0 | box 0 |
| **10.0** | **box 0** — both hold it, first listed wins | **box 1** — the box that starts there |
| 10.001 | box 1 | box 1 |
| 19.999 | box 1 | box 1 |
| **20.0** | box 1 — it owns its north edge | **out of box** — nothing above claims it |

Only the shared edge moves; everything else routes identically.

Accepted forms — a string of the edges that *are* inclusive, a `{edge: bool}`
dict patching the default, a set of names, or an alias:

```python
inclusive="south west start"     # the named edges are inclusive; the rest exclusive
inclusive={"north": False}       # patch one edge, keep the rest
inclusive={"south", "west"}      # set form
inclusive="closed"               # all six (the default — unchanged behaviour)
inclusive="half_open"            # south/west inclusive, north/east exclusive
inclusive="half_open_time"       # ...and `end` exclusive too
```

Any box may override the env-wide default with its own `inclusive` key. That
matters at the **outer rim**: with `north` exclusive nothing claims the topmost
box's north edge, so a query at exactly that latitude is out of box. Close the
last box in each direction to get it back:

```python
boxes = [dict(..., south=0, north=10), dict(..., south=10, north=20, inclusive="closed")]
```

#### Time edges

`start` and `end` work the same way, and matter when your windows abut rather
than sharing one instant. Three consecutive one-hour boxes:

```python
boxes = [dict(start=i * 3600, end=(i + 1) * 3600, south=0, west=0, north=10, east=10)
         for i in range(3)]
```

| timestamp | `closed` (default) | `end` exclusive |
| --- | --- | --- |
| 0 | box 0 | box 0 |
| 3599 | box 0 | box 0 |
| **3600** | **box 0** — both windows hold it, first listed wins | **box 1** — the window that starts there |
| 3601 | box 1 | box 1 |
| 7200 | box 1 | box 2 |
| 10800 | box 2 | box 2 |
| 12000 (past every window) | box 2 | box 2 |

Two things to note. Unlike latitude, a timestamp past the **last** window is not
out of box: nothing holds that instant, so routing falls back to a box that holds
the *position* and lets that box's own `set_out_of_window_behaviour()` decide —
which is what keeps a climatology fallback reachable past the forecast horizon.
And an exclusive `end` is exact to the second: `[start, end)` is `[start, end-1]`,
because timestamps are whole seconds.

Ask for it with `"half_open_time"` (spatial half-open **plus** an exclusive
`end`), or by naming the edges:

```python
MultiAmphiNcData(config, boxes, inclusive="half_open_time")
MultiAmphiNcData(config, boxes, inclusive={"end": False})   # only the time edge
```

`"half_open"` on its own is deliberately **spatial only** — most box sets pin one
instant (`start == end`), and an exclusive `end` would make those windows empty,
which raises rather than routing to nothing.

### `snap_to_grid`: a cheaper dispatch, for a coarser edge

Routing buckets space into 1° cells, and box edges rarely land on cell edges. By
default a cell cut by an edge is resolved per query against the real box bounds —
exact, and a point outside every box reliably returns `out_of_box_value`.

> For an abutting tiling, reach for `inclusive="half_open"` above first: it also
> leaves nothing to resolve, but keeps exact edges and loads no extra data.
> `snap_to_grid` earns its keep when box edges are *ragged* — not aligned to each
> other or to the cell grid — where no inclusivity rule can make them partition.

`snap_to_grid=True` widens every box outward until it swallows whole cells, then
gives each cell to the first box listed that owns it. No cell is ever contested,
so routing becomes a single array read with nothing to resolve:

```python
multi = MultiAmphiNcData(config, boxes, snap_to_grid=True)
multi.requested_boxes[0]   # what you asked for:  north=10.0
multi.boxes[0]             # what you got:        north=11.0
```

The widening reaches the **data**, not just the routing — each box loads the
wider area, otherwise it would own cells it cannot serve.

What it costs: each box loads up to one cell of extra data on every widened side,
and returns data for up to a cell beyond what you asked for instead of
`out_of_box_value`. The finer the tiling, the larger that overhead is relative to
the boxes themselves.

Use it when box edges are ragged and you never treat out-of-box as a signal.
Leave it off when `out_of_box_value` means something to you, or when RAM is the
reason you are using multiple boxes at all.

### `cell_size`: the resolution of the dispatch grid

Routing buckets space and time into cells; `cell_size` sets how big one is, as
`(lat_deg, lon_deg, seconds)`. The default is `(1.0, 1.0, 3600)`.

```python
MultiAmphiNcData(config, boxes, cell_size=(1.0, 1.0, 86400))
```

This is about the **geometry of your box layout**, not the resolution of your
data. In particular the time cell has nothing to do with the product's timestep —
tides at 1 h, wind and wave at 3 h then 6 h, currents at 24 h are all handled
inside each box's sampler. The grid's time axis exists only to separate boxes
that have *different windows*; if every box shares one `[start, end]` it collapses
away entirely and the seconds you pass are irrelevant.

Two things follow.

**Space: pick a cell that divides your tile pitch.** A cell straddling a box edge
has to be resolved per query, and the band of ground where that happens is one
cell wide. Edges at 2.5° on a 1° grid put a full row in that band; a 0.5° grid
makes every edge a cell edge. Against that, the grid is
`(180 / lat_deg) × (360 / lon_deg) × T` cells of 2 bytes, so halving the cell
quadruples it. Both sizes must divide 180° and 360° evenly, or the constructor
raises.

**Time: match your window granularity, not your files.** One box per day wants
`seconds=86400` — at 3600 the grid carries 24× the layers to describe the same
routing. Windows that *abut* want the opposite: a cell holding a window boundary
is claimed by both neighbours, so coarse cells there create contention rather
than removing it. Finer is safer when windows are adjacent, coarser when they are
a day apart.

Routing itself is invariant: every cell size returns the same box for the same
point. Only memory, build time, and how much resolving happens at the seams
change. If you would rather not think about it, `inclusive="half_open"` removes
the seam question at any resolution.

```python
multi.cell_size          # (1.0, 1.0, 3600) -- what is actually in use
```

`env.tidal_currents()`, `env.wind()`, etc. work the same way and keep their
per-box behaviour (tide tiling, climatology fallback). `MultiAmphiNcData` is a
context manager and forwards `close()` to every box's env.

## Using `NcData`

Use `NcData` when you want more manual control.

### Load from basepaths

Use this when the library should discover the right NetCDF files for you.

```python
from amphi_ncdata import NcData

currents = NcData.from_basepaths(
    basepaths=[
        "/data/data_api/FORECAST-CURRENT/",
        "/mnt/store/data_api/HIRES_v3/FORECAST/",
    ],
    start=1781478000,
    end=1781650800,
    south=35.0,
    west=-73.0,
    north=43.0,
    east=-30.0,
    allow_last_available=True,
)

currents.load("ucos")
currents.load("vcos")
currents.set_interpolation_method("linear")

print(currents.at(40.0, -50.0, 1781524800))
```

### Load from explicit filepaths

Use this when you already know the exact files you want.

```python
from amphi_ncdata import NcData

waves = NcData.from_filepaths(
    filepaths=[
        "/tmp/20260614.nc",
        "/tmp/20260615.nc",
    ],
    start=1781391600,
    end=1781478000,
    south=35.0,
    west=-73.0,
    north=43.0,
    east=-30.0,
)

waves.load("swh")
waves.load("mwp")
waves.load("mwd")
waves.set_interpolation_method("linear")

print(waves.at(40.0, -50.0, 1781434800))
```

### Sampling

```python
currents.set_interpolation_method("linear")
value = currents.sample("ucos", 40.0, -50.0, 1781524800)
```

### Batch sampling

For large point sets, batch APIs avoid Python per-sample overhead:

```python
lats = [39.5, 40.0, 40.5]
lons = [-51.0, -50.0, -49.0]
timestamps = [1781521200, 1781524800, 1781528400]

u_values = currents.sample_batch("ucos", lats, lons, timestamps)
uv_values = currents.at_batch(lats, lons, timestamps)
```

Notes:

- `sample_batch()` returns one NumPy array
- `at_batch()` returns one NumPy array shaped like `(sample_count, variable_count)` when multiple variables are loaded
- use batch APIs when you already hold many independent points; for sequential, point-by-point access (e.g. routing), `at()` / `sample()` are the right call

#### Reusing an output buffer with `out=`

`at_batch()` and `sample_batch()` accept an optional `out=` buffer: results are
written into it in place and the same array is returned — no per-call allocation.
Build the correctly-shaped, correctly-typed buffer once with `empty_batch()` and
reuse it across calls (e.g. re-sampling the same points on every step of a loop):

```python
import numpy as np

lats = [39.5, 40.0, 40.5]
lons = [-51.0, -50.0, -49.0]
timestamps = [1781521200, 1781524800, 1781528400]

buf = currents.empty_batch(len(lats))                 # (n, 2) float32 for a 2-variable dataset
currents.at_batch(lats, lons, timestamps, out=buf)    # returns buf, filled in place

sbuf = np.empty(len(lats), dtype=np.float32)          # (n,) float32 for a single variable
currents.sample_batch("ucos", lats, lons, timestamps, out=sbuf)
```

`out` must be a C-contiguous `float32` array of shape `(n, variable_count)` for
`at_batch()` — or `(n,)` when the dataset has a single loaded variable, and always
`(n,)` for `sample_batch()`. `empty_batch(n)` returns exactly the shape `at_batch()`
needs. Passing `out=` implies a NumPy-array return (the `as_numpy` flag is ignored).

**When to batch:** the batch APIs run the same compiled sampler as `at()` /
`sample()`, but in one call. They carry a small fixed per-call cost, so for just a
few points a plain `at()` loop is still faster; batching pays off from roughly a
few dozen points upward and scales to millions.

Supported interpolation values for `set_interpolation_method()`:

- `"nearest"`
- `"linear"`

### No-data near the coast: `skip_nan` and `set_no_data_value()`

Coastal and land cells carry the product's `_FillValue`, which decodes to NaN.
Bilinear interpolation reads four corners, so **one** no-data corner is enough to
make the whole sample NaN — which is why a point near the coast can return a
value under `"nearest"` and nothing under `"linear"`.

There are two independent controls, and they compose.

**1. Interpolate from the corners that do hold data:**

```python
currents.set_interpolation_method("linear", skip_nan=True)
```

The no-data corners drop out and the survivors are reweighted: two valid corners
and two no-data corners interpolate from the two valid ones. All four no-data
still gives NaN — there is nothing to interpolate from. The same rule extends to
the time blend, so a timestep that is wholly no-data no longer poisons the one
next to it. Water away from the coast is unaffected: a stencil with no no-data in
it interpolates exactly as before.

`skip_nan` is remembered across calls, so pass it only when changing it. It has
no effect under `"nearest"`, which reads a single cell.

**2. Give the no-data mask a value:**

```python
waves.set_no_data_value(0.0)                                   # every variable
waves.set_no_data_value((0.0, 0.0, float("nan")))              # swh, mwp, mwd
waves.set_no_data_value({"mwd": float("nan")}, when="after_interpolation")
```

The value is per variable — a scalar applies to all of them, a tuple/list is
positional over `loaded_variables` (`waves` is height, period, direction; `wind`
is u, v), and a dict sets only the variables it names. Each entry is a float, or
NaN to leave that variable as it was.

`when` decides where the substitution happens, and the two are **not**
equivalent:

| `when` | What happens | Use it when |
| --- | --- | --- |
| `"before_interpolation"` | No-data cells enter the stencil as `value`, and all four corners interpolate normally. | The mask means a real physical value (nil current inside a harbour), so it should bleed into nearby water the way the data would. |
| `"after_interpolation"` | Interpolate as usual, then replace the result if it came out NaN. Applies **whatever `skip_nan` is set to**. | You need a number everywhere — it is the backstop for what `skip_nan` cannot rescue. |

`after_interpolation` never touches valid water, and it does not override
`set_out_of_window_behaviour()`: it replaces *no-data*, not an out-of-window
query, which stays with the out-of-window policy.

Two edges worth knowing. `skip_nan` reweights by interpolation weight, so a query
sitting *exactly* on a no-data grid node still gives NaN — the valid corners
carry zero weight there. And `skip_nan` never invents data. For a value in those
cases, combine the two:

```python
currents.set_interpolation_method("linear", skip_nan=True)
currents.set_no_data_value(0.0, when="after_interpolation")
```

`no_data_values` reports what is currently set, as `{variable: (value, when)}`.
Both settings are available on `NcData`, on every `AmphiNcData` accessor
(including `tidal_currents()`), and on `MultiAmphiNcData` routed datasets, where
they apply to every box.

### Lazy loading and cache

```python
cache = {"size": 128 * 1024 * 1024, "nelems": 20000, "preemption": 0.8}

tidal = NcData.from_basepaths(
    basepaths=["/data/data_api/FORECAST-TIDE/"],
    start=1781478000,
    end=1781650800,
    south=35.0,
    west=-73.0,
    north=43.0,
    east=-30.0,
    allow_last_available=True,
)

tidal.load("utotal", lazy=True, cache=cache)
tidal.load("vtotal", lazy=True, cache=cache)
```

Loaded variables keep their **raw on-disk dtype** in RAM (e.g. int16 for currents/tides/waves);
`scale_factor`, `add_offset`, and `_FillValue` are applied inside the Cython sampler at lookup
time. This halves resident memory for packed products versus decoding to float32 up front, with
no loss of precision. (The `decode=` argument to `load()` is retained for API compatibility but
is now a no-op — storage is always the raw dtype.)

### Inspect loaded state

```python
print(currents.resolved_files)
print(currents.loaded_variables)
print(currents.time_range)
print(currents.bounds)
print(currents.contains_timestamp(1781524800))
```

### Context manager

```python
from amphi_ncdata import NcData

with NcData.from_basepaths(
    basepaths=["/data/data_api/FORECAST-CURRENT/"],
    start=1781478000,
    end=1781650800,
    south=35.0,
    west=-73.0,
    north=43.0,
    east=-30.0,
) as data:
    data.load("ucos")
    data.set_interpolation_method("linear")
    print(data.sample("ucos", 40.0, -50.0, 1781524800))
```

### Out-of-window behaviour

A query is "out of window" when its timestamp is outside the loaded time range,
or its `lat`/`lon` is outside the loaded bounding box. Choose what `at()` /
`sample()` (and the batch variants) return in that case:

```python
wind.set_out_of_window_behaviour("nearest")   # default: clamp to the nearest
                                               # available timestep / edge cell
wind.set_out_of_window_behaviour("nan")        # return NaN (handy to detect gaps)
wind.set_out_of_window_behaviour("value", 0.0) # return a fixed fill value
```

This is independent of NaN coming from the data itself (land / fill gaps), which
is always returned as NaN regardless of the mode. The setting applies equally to
out-of-range time and out-of-box space.

## Timestamps And Coordinates

Accepted timestamp format:

- Unix timestamp integer in seconds

Sampling coordinate format:

- `lat` and `lon` are passed as separate numeric arguments

## Errors

The package raises:

- `AmphiDiscoveryError`: file discovery, open, or setup problems
- `AmphiVariableError`: variable loading or sampling problems
- `AmphiNativeError`: lower-level native errors

**Missing data is not an error.** A window with no files, partial coverage, or only
a stale run on disk sets a code on `dataset.warnings` and still returns a usable
dataset — see [warnings](#missing-data-is-reported-not-raised-warnings). So this
does *not* raise; it returns an empty dataset whose queries yield NaN:

```python
from amphi_ncdata import NcData, WARN_NO_FILES

wind = NcData.from_basepaths(
    ["/tmp/does-not-exist"],
    start=0, end=1, south=0, west=0, north=1, east=1,
)
wind.is_empty                       # True
wind.warnings == WARN_NO_FILES      # True

# Pass on_warning="raise" for the old behaviour:
NcData.from_basepaths(["/tmp/does-not-exist"], start=0, end=1,
                      south=0, west=0, north=1, east=1,
                      on_warning="raise")   # -> AmphiDiscoveryError
```

What *does* raise, regardless of `on_warning`:

```python
from amphi_ncdata import AmphiVariableError

try:
    env.wind(u_var="not_a_variable")     # request is wrong, not the archive
except AmphiVariableError as exc:
    print(exc)
```

One caveat: variable names are validated against the file, so on an **empty**
dataset there is nothing to validate against — a typo'd name is accepted there and
only surfaces once the window has files. Check `warnings` before concluding that a
run of NaN means a bad variable name.

## Performance

- `at()` and `sample()` are backed by a compiled sampler, so single-point lookups are fast. This is the intended path for routing algorithms (Dijkstra / A* / isochrone / DP) that sample one point at a time, where each step depends on the previous one.
- Eager-loaded variables (the default) stay in memory, so sampling never touches disk. Use `lazy=True` only for variables too large to keep resident; lazy sampling reads from disk and is much slower per point.
- `at_batch()` / `sample_batch()` run the same compiled sampler over many points in one call. They carry a small fixed per-call cost, so use them once you have roughly a few dozen points or more (below that, an `at()` loop wins). Pass `out=` (from `empty_batch()`) to avoid per-call allocation when sampling repeatedly.
- An eager-loaded dataset closes the files it read as soon as it has loaded them, and reopens them on its own if you `load()` again. An open file would keep every chunk it decompressed cached — the WAVERYS hindcast stores each variable as one 26 MB chunk — and that is how a 237-box wave env used to hold 4.2 GiB of memory for 32 MiB of data. `MultiAmphiNcData` keeps the files its boxes share open while it builds them (within 512 MiB), so each is decompressed once, then closes them all. Lazily loaded variables and tides keep their files, since they read them per query.
- The climatology fallback behind `env.wind()` / `env.waves()` runs compiled too: past the horizon, `at()`, `at_batch()` and the `MultiAmphiNcData` dispatch hand the point to the climatology's own compiled sampler without leaving C, so a route that crosses the horizon costs no more per point than one that does not. A fallback you wire yourself with `set_out_of_window_behaviour("fallback", other, dim="time")` gets the same, provided `other` is eager-loaded, holds the same variables in the same order, and has no `"fallback"` on space. Otherwise it is still correct but resolved in Python, which costs the batch calls their compiled loop and a multi-box env its compiled dispatch.
- A window that [spans a change of file format](#windows-that-span-a-change-of-file-format) samples at the same speed as one that does not: each part keeps its own compiled sampler, and choosing the part costs a comparison or two. Only points in the hours across the switch itself read two grids, which costs a batch call there about 15% more.

## Practical Notes

- A wheel is specific to an OS, Python version, and architecture.
- If you change Python versions, you may need a new wheel.
- `AmphiNcData(...)` accepts missing path keys, but a dataset helper such as `env.wind()` or `env.waves()` still raises immediately if that dataset cannot be resolved.
- `AmphiNcData` is the best starting point unless you specifically need file-level control.
- This package focuses on point sampling, not bulk NumPy-style array workflows.
