diff --git a/CHANGELOG.md b/CHANGELOG.md index 0405545..7d075f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,8 +12,94 @@ and was not. Denominators are computed once per variable and shared across all ranges defined on it. +- **`group` support in `pcount`** — the stat now honours the `group` field, reporting the + percentage of valid samples *within each group* and returning a group dimension. The + denominator is the number of timesteps belonging to the group, not the length of the + whole record, so a variable with no gaps reads 100 in every group. + ### Bug Fixes +- **`statdir` with more than one function raised `MergeError`** — every stat names its + output after the variable it reduced, so `funcs: [mean, max]` produced two datasets + that both contained `hs` and the merge of the two sector results failed with + *conflicting values for variable 'hs'*. The stat was therefore usable only with the + single function it exists to avoid. Delegated outputs are now suffixed with the + function that produced them (`hs_mean`, `hs_max`), including when only one function is + requested, so a name never depends on the length of `funcs`. The pipeline no longer + appends `_statdir` on top of a name that already carries the name of the delegated + stat, which would have given `hs_mean_statdir_direc`; the written names are + `hs_mean_direc` and `hs_max_direc` — identical to those of the equivalent per-function + calls with `nsector` set, and they resolve through `attributes.yml` like any other + output. Single-function `statdir` output previously written as `hs_statdir` will be + written as `hs_mean_direc`; when appending to an existing store the old variable is not + removed automatically. +- **`nsector` on a `statdir` call binned the data twice** — `nsector` is a pipeline-level + field, so it was consumed by the pipeline's own sectoriser and never reached the stat, + which then applied its own default of 4 sectors inside each of the outer ones. The + documented eight-sector `statdir` example failed with *conflicting sizes for dimension + 'direction'*. The sector settings are now forwarded to stats that sectorise internally + instead of being applied around them, so `nsector: 8` yields the eight sectors asked + for. +- **`statdir` accepts `group` again** — it forwarded `group` to its delegated stats + through `**kwargs` but declared no `group` parameter of its own, so the new + signature-based guard classed it as unable to group and rejected the call. `group` is + now an explicit parameter passed through to every delegated stat, and a call naming a + stat that cannot honour it (`funcs: [mean, rpv]` with `group: month`) fails up front, + naming the offender, rather than silently returning an ungrouped result for it. +- **`output.chunks` was silently ignored** — the pipeline read the output chunking as + `self.config.output.__dict__.get("chunks", {})`, but `OutputConfig` declared no + `chunks` field, and Pydantic had already dropped the key while validating the config. + A user who set `chunks:` under `output:` got the default chunking of the finalised + dataset — typically a single chunk per variable — with no warning, however large the + grid. `chunks` is now a declared `dict[str, int]` field, read by normal attribute + access, and reaches `finalise`, which applies it to the dimensions present in the + output and clears the stale chunk encoding inherited from the source. +- **Grouped variables were written with no attributes at all** — the pipeline appends the + group to the output variable name, so a `group: month` call on `mean` writes + `hs_mean_month`. Attribute resolution split that name on `_` and read the last token as + the statistic, found `month` in neither the `stats` nor the `data_vars` map, and gave + up: the variable was delivered with no `units`, no `standard_name` and no `long_name`, + while its ungrouped sibling `hs_mean` was fully described. The group suffix — and the + `_direc` suffix a directionally sectorised call appends after it — are now stripped + before the statistic and the parent variable are looked up, so `hs_mean_month` and + `hs_mean_month_direc` both resolve through `hs_mean` and carry `units: m`. A grouped + variable additionally has its `long_name` prefixed with the matching adjective, taken + from an explicit mapping rather than built by string arithmetic (`month` → *monthly*, + `season` → *seasonal*, `year` → *annual*; the previous, unreachable attempt at this + would have produced "seasonly"). Only a *trailing* group token counts as a group, so a + parent variable whose own name contains a group word is not mangled. The override that + keeps a probability or a percentage dimensionless still applies on top, so + `hs_gt_1p5_range_probability_month` remains `units: "1"`. Values are unchanged, but + grouped fields already written carry no metadata and will not gain any until they are + rewritten or re-consolidated. +- **`count` no longer inherits the units of its parent variable** — `hs_count` was written + with `units: m` although it is a number of samples. `count` is now declared in the + `stat_units` section of `attributes.yml` as `"1"`, the CF/UDUNITS spelling of a + dimensionless quantity, alongside `dist`, whose joint histograms are raw counts and had + exactly the same problem. Values are unchanged; only metadata is affected. +- **Grouped `pcount` output was the ungrouped value** — `pcount` accepted `group` and + silently discarded it, yet the pipeline still appended the group to the output variable + name, so a `group: month` config wrote the whole-record coverage under a name like + `hs_pcount_month`, with no `month` dimension. Now that the stat honours `group` the name + matches the data. Any grouped `pcount` field produced before this change is wrong and + should be regenerated. +- **Requesting `group` on a stat that cannot honour it is now an error** — `rpv`, `hmo` + and `distribution3_timestep` absorbed `group` into `**kwargs` and ignored it while the + pipeline appended it to the output variable name, producing an ungrouped result wearing + a grouped name. Such a call now raises before any data is loaded, naming the stat and + listing the stats that do support grouping. Support is determined by inspecting the + stat's signature for an explicit `group` parameter, so it cannot drift out of sync with + the code. This is a behaviour change callers may notice: a config that previously ran to + completion with a meaningless `group` on one of these stats will now fail, and the fix + is to drop the `group` field from that call. Grouping is deliberately not implemented + for the three — a per-group extreme-value fit uses a fraction of the samples and is not + statistically defensible, and grouping is not well defined for a Welch spectrum or for + a histogram accumulated over windows that do not align with time groups. +- **`rpv` timestep derivation no longer assumes nanosecond time** — the sampling interval + was computed as `float(np.diff(time)[0]) / 3.6e12`, which raises `TypeError` under + NumPy 2 for any `datetime64` resolution other than nanoseconds (and would silently + mis-scale if it did not). It now goes through `pd.Timedelta`, matching the fix already + applied in `exceedance`. - **`uorb` no longer returns a value outside the water column** — the evaluation point is a height above the seabed, `z_bed`, and with `reference: surface` it is derived as `depth - z`, so `z: 15` in 3 m of water asked for a point 12 m *below* the bed. The @@ -59,12 +145,29 @@ ### Breaking Changes +- **Unknown fields under `output` are now rejected** — `OutputConfig` inherited Pydantic's + default of dropping keys it does not know, which is how `output.chunks` came to be + written in configs, accepted by validation and never applied. Any misspelling had the + same failure mode: `chunk:` for `chunks:`, or `consolidated:` for `consolidate:`, left + the option unset while the run reported success. The model now forbids extras, so such + a config fails validation naming the offending key. A config carrying a stray key under + `output` that used to run will now fail, and the fix is to correct or remove the key — + the setting it was meant to express was never in force. - **Fractional exceedance thresholds are named with `p`** — `exceedance`/`nonexceedance` output variables now write the decimal point as `p`, so a threshold of 1.5 gives `hs_1p5` rather than `hs_1.5`, keeping the name usable as an identifier and by attribute access. Integral thresholds are unchanged (`2.0` still gives `hs_2`). Configs using fractional thresholds will write new variable names; when appending to an existing store the old dotted variables are not removed automatically. +- **Default `range_probability` labels are named with `p`** — a `data_ranges` entry with + no `label` built its output name from the raw bounds, so a fractional bound produced + `hs_1.5_to_max`, a name that is not a valid identifier and is awkward in path-like + contexts. The decimal point is now written as `p` (`hs_1p5_to_max`), matching + `exceedance`; the formatting is a single helper shared by both ops so the two cannot + drift apart. Integral bounds are unchanged (`hs_1_to_2`), and entries that set `label` + explicitly are unaffected. Configs relying on the default label with a fractional bound + will write new variable names; when appending to an existing store the old dotted + variables are not removed automatically. --- diff --git a/docs/configuration/output.md b/docs/configuration/output.md index f5a7e54..f0e6d53 100644 --- a/docs/configuration/output.md +++ b/docs/configuration/output.md @@ -7,12 +7,30 @@ Controls where and how results are written. | Field | Type | Default | Description | |---|---|---|---| | `outfile` | string | — | Output file path. Extension determines format: `.nc` for NetCDF4, `.zarr` for Zarr. Supports remote paths (`gs://`, `s3://`). **Required.** | +| `chunks` | dict | `{}` | Chunk sizes applied to the finalised dataset before writing: `{dim: size}`. Only dimensions present in the output are used, so a chunk on a dimension the statistics reduced away (e.g. `time`) is ignored. See [Output chunking](#output-chunking) below. | | `global_attrs` | dict | `{}` | Global dataset attributes to add to or override the defaults. See [Global attributes](#global-attributes) below. | | `append` | bool | `false` | Add variables to an existing Zarr store rather than overwriting it. See [Parallel Zarr writes](#parallel-zarr-writes) below. | | `consolidate` | bool | `false` | Run `zarr.consolidate_metadata()` after writing. See [Parallel Zarr writes](#parallel-zarr-writes) below. | | `mask` | dict | `null` | Optional spatial mask applied to all output variables before writing. See [Masking](#masking) below. | | `updir` | string | `null` | Optional remote directory to upload the written output to after the run. See [Uploading to a remote directory](#uploading-to-a-remote-directory) below. | +Unknown fields are rejected: a misspelled option would otherwise be dropped in +silence and the output written as if it had never been set. + +## Output chunking + +Sets the chunking of the written store. Chunks are applied to the finalised dataset, +after the statistics have reduced the time dimension away, so the sizes refer to the +dimensions of the *output* grid rather than of the source. + +```yaml +output: + outfile: out.zarr + chunks: + latitude: 500 + longitude: 500 +``` + ## Masking An optional spatial mask can be applied to all output variables before writing. The mask is derived from a variable in the **source** dataset and broadcast automatically across any extra dimensions (e.g. time, quantile, direction) in the output. diff --git a/docs/ops/aggregations.md b/docs/ops/aggregations.md index 97d77f9..453b275 100644 --- a/docs/ops/aggregations.md +++ b/docs/ops/aggregations.md @@ -106,8 +106,13 @@ Percentage of non-NaN values (0–100). Indicates data coverage. - func: pcount dim: time data_vars: [hs, tp] + group: month # optional ``` +With `group` set the percentage is computed within each group: the denominator is the +number of timesteps belonging to that group, not the length of the whole record, so a +variable with no gaps reads 100 in every group. + --- ## API reference diff --git a/docs/ops/directional.md b/docs/ops/directional.md index d58cf58..dd4fc1e 100644 --- a/docs/ops/directional.md +++ b/docs/ops/directional.md @@ -14,14 +14,16 @@ Bins the dataset by direction, applies each listed function to each sector, and | `dir_var` | str | `"dpm"` | Variable name used for directional binning. Must be in the dataset. | | `nsector` | int | `4` | Number of equally-spaced sectors. `8` → 45° sectors, `16` → 22.5°, etc. | -Additional keyword arguments (e.g. `group`) are forwarded to each stat function. +Additional keyword arguments (e.g. `group`) are forwarded to each stat function. They are shared by every function in `funcs`, so statistics needing different arguments have to be requested in separate calls. When `group` is set, every function in `funcs` must support grouping. -Output variables are named as `{variable}_{func}` and gain a `direction` coordinate with sector-centre values (0°, 45°, …). +The stat names its outputs `{variable}_{func}` — a single function is suffixed too, so a name never depends on how many functions were asked for — and the pipeline then appends the group (when set) and `_direc`. A `funcs: [mean, max]` call over `hs` writes `hs_mean_direc` and `hs_max_direc`, exactly the names the equivalent `nsector` calls produce. Every variable gains a `direction` coordinate holding the sector-centre values (0°, 45°, …). + +The directional variable must be present in the data the call sees, so include it in `data_vars` when that list is explicit. It is itself reduced by each function, giving `dpm_mean_direc` and so on. ```yaml - func: statdir dim: time - data_vars: [hs, tp] + data_vars: [hs, tp, dpm] dir_var: dpm nsector: 8 funcs: [mean, max] @@ -47,14 +49,15 @@ Data points outside a sector are set to NaN before the stat function is called, Setting `nsector` on a regular call (e.g. `func: mean, nsector: 8`) is handled at the pipeline level and applies a single function. `statdir` lets you apply **multiple** functions in one pass while keeping the same directional binning, which is more efficient when you need several statistics per sector. +The two routes are otherwise interchangeable: they bin the data identically and write identical variable names, so `statdir` adds no capability of its own — it exists to share one masking of the data and one pass over it between several statistics. Prefer `nsector` on individual calls unless that matters, or when the statistics need different keyword arguments (which `statdir` cannot express, since its extra arguments go to every function). + ```yaml # Equivalent to two separate calls with nsector, but done in one pass: - func: statdir dim: time - data_vars: [hs] + data_vars: [hs, dpm] nsector: 16 - funcs: [mean, max, quantile] - q: [0.90, 0.99] + funcs: [mean, max] ``` --- diff --git a/docs/ops/distribution.md b/docs/ops/distribution.md index d5aeb68..5b0baca 100644 --- a/docs/ops/distribution.md +++ b/docs/ops/distribution.md @@ -82,7 +82,7 @@ Output variable: `dist2` with dimensions `(var1, var2)`. Memory-efficient variant of `distribution3` that accumulates histogram counts in time chunks. Suitable for very long timeseries where loading the full dataset at once would exceed memory limits. -All parameters are the same as `distribution3`, plus: +All parameters are the same as `distribution3` except `group`, which is not supported, plus: | Parameter | Type | Default | Description | |---|---|---|---| @@ -98,13 +98,17 @@ All parameters are the same as `distribution3`, plus: bins1: {start: 0, step: 0.5} bins2: {start: 0, step: 1.0} bins3: {start: 0, stop: 360, step: 45} - group: month freq: 10d # load 10 days at a time ``` !!! tip Prefer `distribution3` for datasets that fit comfortably in memory — it is faster because it avoids repeated I/O. Use `distribution3_timestep` only for multi-decade datasets. +!!! warning + `distribution3_timestep` does not support `group`: the accumulation runs over + `freq`-sized windows of the record, which do not align with time groups. A call + setting `group` raises an error — use `distribution3` with `group` instead. + --- ## API reference diff --git a/docs/ops/probability.md b/docs/ops/probability.md index 0d16888..7cf4364 100644 --- a/docs/ops/probability.md +++ b/docs/ops/probability.md @@ -23,7 +23,7 @@ Each item in `data_ranges` is a dict with the following fields: | `stop` | float \| null | — | Upper bound. `null` means no upper bound. | | `left` | `"closed"` \| `"open"` | `"closed"` | Whether the lower bound is inclusive (`>=`) or exclusive (`>`). | | `right` | `"closed"` \| `"open"` | `"closed"` | Whether the upper bound is inclusive (`<=`) or exclusive (`<`). | -| `label` | str | `{var}_{start}_to_{stop}` | Output variable name. | +| `label` | str | `{var}_{start}_to_{stop}` | Output variable name. An absent bound is written `min`/`max`, and a decimal point is written `p` so the name stays a valid identifier (`start: 1.5`, `stop: null` → `hs_1p5_to_max`). | ### Output diff --git a/docs/ops/rpv.md b/docs/ops/rpv.md index 6eea1fd..c0ad85c 100644 --- a/docs/ops/rpv.md +++ b/docs/ops/rpv.md @@ -10,10 +10,14 @@ Estimates the value expected to be exceeded once every N years using peaks-over- | `percentile` | float | `95` | Percentile used to define the POT threshold. | | `distribution` | str | `"gumbel_r"` | Any `scipy.stats` continuous distribution name. | | `duration` | float | `24` | Minimum peak separation in hours. | -| `group` | str | `null` | Time grouping (rarely used with RPV). | The output has a `period` dimension with the return period values (in years). +!!! warning + `rpv` does not support `group`. Fitting an extreme-value distribution within a single + time group uses only a fraction of the samples (a twelfth of them for `group: month`), + which is not statistically defensible, so a call setting `group` raises an error. + ```yaml - func: rpv dim: time diff --git a/gridstats/attributes.yml b/gridstats/attributes.yml index 938eefb..3bbe213 100644 --- a/gridstats/attributes.yml +++ b/gridstats/attributes.yml @@ -219,8 +219,14 @@ stats: # Units of statistics that do not preserve the units of their parent variable. # Resolved against every '_'-separated token of the output variable name, so # suffixed names such as `hs_gt_1p5_range_probability_month` resolve as well. +# "1" is the CF/UDUNITS spelling of a dimensionless quantity. stat_units: + # Counts of samples: a number, not a quantity in the parent's units. + count: "1" + + dist: "1" + pcount: "%" probability: "1" diff --git a/gridstats/config.py b/gridstats/config.py index 102a5d8..2a7202c 100644 --- a/gridstats/config.py +++ b/gridstats/config.py @@ -75,9 +75,19 @@ class ThresholdMaskConfig(BaseModel): class OutputConfig(BaseModel): - """Configuration for pipeline output.""" + """Configuration for pipeline output. + + Unknown fields are rejected: silently ignoring a misspelled option here + means the output is written with a setting the user believes is in force. + """ + + model_config = ConfigDict(extra="forbid") outfile: str + #: Chunk sizes applied to the finalised dataset before it is written, + #: ``{dim: size}``. Only dimensions present in the output are used; stale + #: chunk encoding inherited from the source is cleared. + chunks: dict[str, int] = {} #: Optional remote directory to upload the written output to after the run. #: The output is copied under its basename (e.g. ``outfile: ./scratch/hs.zarr`` #: with ``updir: gs://bucket/stats`` -> ``gs://bucket/stats/hs.zarr``). Any diff --git a/gridstats/ops/_naming.py b/gridstats/ops/_naming.py new file mode 100644 index 0000000..827e529 --- /dev/null +++ b/gridstats/ops/_naming.py @@ -0,0 +1,27 @@ +"""Shared helpers for building output variable names. + +Several ops derive an output variable name from a numeric bound — a threshold in +:mod:`gridstats.ops.exceedance`, the ends of a range in +:mod:`gridstats.ops.probability`. They must agree on how a number is spelled, so +the formatting lives here rather than being duplicated per module. +""" +from __future__ import annotations + +__all__ = ["format_number"] + + +def format_number(value: float) -> str: + """Format a numeric bound for use in an output variable name. + + The value is rendered with ``'%g'`` (so ``2.0`` loses its trailing zero) and + the decimal point is replaced by ``'p'`` (for "point"), keeping the name a + valid identifier that can be reached by attribute access and safe to use in + path-like contexts: ``hs_1p5`` rather than ``hs_1.5``. + + Args: + value: The number to render. + + Returns: + The formatted string: ``1.5`` -> ``'1p5'``, ``2.0`` -> ``'2'``. + """ + return f"{value:g}".replace(".", "p") diff --git a/gridstats/ops/aggregations.py b/gridstats/ops/aggregations.py index 8009d07..ae8177f 100644 --- a/gridstats/ops/aggregations.py +++ b/gridstats/ops/aggregations.py @@ -147,7 +147,9 @@ def quantile( @register_stat("pcount") -def pcount(data: xr.Dataset, *, dim: str = "time", **kwargs) -> xr.Dataset: +def pcount( + data: xr.Dataset, *, dim: str = "time", group: str | None = None, **kwargs +) -> xr.Dataset: """Percentage of non-NaN values along a dimension. Values are in [0, 100]. Useful for reporting data coverage. @@ -155,11 +157,24 @@ def pcount(data: xr.Dataset, *, dim: str = "time", **kwargs) -> xr.Dataset: Args: data: Input dataset. dim: Dimension to reduce along. + group: Time component for grouped climatology: `'month'`, `'season'`, or + `'year'`. When set, the percentage is computed *within each group* — + the denominator is the number of timesteps belonging to that group, + not the length of the whole record — and the output gains a + dimension named after the group. Returns: - Dataset with values in [0, 100], carrying ``units: '%'``. + Dataset with values in [0, 100], carrying ``units: '%'``. Gains a + `group` dimension when `group` is set. """ - dsout = 100 * data.count(dim) / data[dim].size + # The denominator is the number of timesteps available to each group, so a + # variable with no gaps reads 100 % in every group. + if group is None: + nsamples = data[dim].size + else: + nsamples = _groupby(xr.ones_like(data[dim], dtype="int64"), group).sum(dim) + + dsout = 100 * _groupby(data, group).count(dim) / nsamples # A percentage of valid samples, not the parent quantity, so the units and # descriptive names inherited from it through the arithmetic do not apply. diff --git a/gridstats/ops/directional.py b/gridstats/ops/directional.py index a9f4c66..1802e75 100644 --- a/gridstats/ops/directional.py +++ b/gridstats/ops/directional.py @@ -160,6 +160,7 @@ def statdir( dir_var: str = "dpm", nsector: int = 4, dim: str = "time", + group: str | None = None, **kwargs: Any, ) -> xr.Dataset: """Apply multiple stat functions over directional sectors. @@ -168,18 +169,50 @@ def statdir( function in ``funcs`` to each sector, then concatenates results along a new 'direction' dimension. + Each delegated statistic returns its result under the name of the variable + it was computed from, so the outputs of two functions would collide on + merge. They are therefore suffixed with the function that produced them: + ``funcs: [mean, max]`` over ``hs`` yields ``hs_mean`` and ``hs_max``. The + suffix is applied for a single function too, so a name never depends on how + many functions were requested. + Args: - data: Input dataset. Must contain ``dir_var``. + data: Input dataset. Must contain ``dir_var`` — when the call selects + ``data_vars`` explicitly, the directional variable has to be among + them. funcs: Names of registered stat functions to apply per sector. dir_var: Name of the directional variable used for binning. nsector: Number of equally-spaced directional sectors (default 4). dim: Time dimension name. - **kwargs: Forwarded to each stat function. + group: Time component for a grouped climatology (``'month'``, + ``'season'``, ``'year'``), passed through to every delegated stat. + All of them must support grouping. + **kwargs: Forwarded to each stat function. They are shared by every + function in *funcs*; per-function keyword arguments are not + supported. Returns: - Dataset with a 'direction' dimension containing sector-centre values. + Dataset with one variable per ``{variable}_{func}`` combination, each + gaining a 'direction' dimension holding the sector-centre values. + + Raises: + ValueError: If *funcs* is empty, if ``dir_var`` is missing from the + dataset, or if *group* is set and one of *funcs* cannot honour it. + + Note: + This stat adds no capability that the pipeline's own ``nsector`` field + does not already provide: ``func: mean`` with ``nsector: 8`` produces + exactly the same variable, under exactly the same output name, as + ``func: statdir`` with ``funcs: [mean]`` and ``nsector: 8``. Its reason + to exist is efficiency — several statistics share one directional + masking of the data and one pass over it, instead of one full pass per + call. Prefer ``nsector`` on individual calls unless that matters, or + when the statistics need different keyword arguments. """ - from gridstats.registry import get_stat + from gridstats.registry import get_stat, list_group_stats, supports_group + + if not funcs: + raise ValueError("statdir requires at least one function in 'funcs'.") if dir_var not in data: raise ValueError( @@ -187,6 +220,17 @@ def statdir( f"Available: {list(data.data_vars)}" ) + # Validate the delegated stats up front so a config error is not reported + # from inside the sector loop, after work has already been done. + if group is not None: + unsupported = [f for f in funcs if not supports_group(f)] + if unsupported: + raise ValueError( + f"statdir was given group='{group}', but these functions do not " + f"support grouping: {', '.join(unsupported)}. Stats that support " + f"grouping: {', '.join(list_group_stats())}." + ) + dirs = data[dir_var] dsector = 360.0 / nsector sector_centres = np.linspace(0, 360 - dsector, nsector) @@ -205,7 +249,13 @@ def statdir( func_results = [] for func_name in funcs: fn = get_stat(func_name) - func_results.append(fn(sector_ds, dim=dim, **kwargs)) + result = fn(sector_ds, dim=dim, group=group, **kwargs) + # Every stat names its output after the variable it reduced, so the + # results of two functions would collide here. Disambiguate them by + # the function that produced them before merging. + func_results.append( + result.rename({v: f"{v}_{func_name}" for v in result.data_vars}) + ) sector_results.append(xr.merge(func_results)) dsout = xr.concat(sector_results, dim="direction").assign_coords( diff --git a/gridstats/ops/distribution.py b/gridstats/ops/distribution.py index 24011d1..926b55a 100644 --- a/gridstats/ops/distribution.py +++ b/gridstats/ops/distribution.py @@ -228,7 +228,6 @@ def distribution3_timestep( isdir2: bool = False, isdir3: bool = True, freq: str = "30d", - group: str | None = None, **kwargs: Any, ) -> xr.Dataset: """3-D joint histogram accumulated over time chunks to limit memory use. @@ -253,7 +252,12 @@ def distribution3_timestep( isdir2: Whether `var2` is directional. isdir3: Whether `var3` is directional (default `True`). freq: Pandas-compatible frequency string for time-chunking (e.g. `'30d'`, `'1ME'`). - group: Time component to group by (e.g. `'month'`). + + Note: + Grouping is not supported: the accumulation runs over `freq`-sized + windows of the record, which do not align with time groups. Requesting + `group` on this stat raises rather than silently returning the ungrouped + histogram. Use `distribution3` with `group` instead. Returns: Dataset with accumulated joint distribution counts, same structure as diff --git a/gridstats/ops/exceedance.py b/gridstats/ops/exceedance.py index a5ae1b2..84d923e 100644 --- a/gridstats/ops/exceedance.py +++ b/gridstats/ops/exceedance.py @@ -7,6 +7,7 @@ import pandas as pd import xarray as xr +from gridstats.ops._naming import format_number from gridstats.registry import register_stat logger = logging.getLogger(__name__) @@ -18,16 +19,6 @@ def _duration_to_hours(duration: str) -> float: return td.total_seconds() / 3600 -def _format_threshold(threshold: float) -> str: - """Format a threshold for use in an output variable name. - - The decimal point is replaced by ``'p'`` (for "point") so the name stays a - valid identifier and can be reached by attribute access: ``hs_1p5`` rather - than ``hs_1.5``. Integral thresholds are unaffected (``2.0`` -> ``'2'``). - """ - return f"{threshold:g}".replace(".", "p") - - def _runs_above_duration( condition: np.ndarray, dt_hours: float, min_hours: float ) -> np.ndarray: @@ -178,7 +169,7 @@ def exceedance( condition = cmp(threshold) & (data <= maxval) dsout = _probability_of_occurrence(data, condition, duration, dim, group) return dsout.rename( - {v: f"{v}_{_format_threshold(threshold)}" for v in dsout.data_vars} + {v: f"{v}_{format_number(threshold)}" for v in dsout.data_vars} ) @@ -215,5 +206,5 @@ def nonexceedance( condition = cmp(threshold) dsout = _probability_of_occurrence(data, condition, duration, dim, group) return dsout.rename( - {v: f"{v}_{_format_threshold(threshold)}" for v in dsout.data_vars} + {v: f"{v}_{format_number(threshold)}" for v in dsout.data_vars} ) diff --git a/gridstats/ops/frequency_domain.py b/gridstats/ops/frequency_domain.py index e663177..1706a0b 100644 --- a/gridstats/ops/frequency_domain.py +++ b/gridstats/ops/frequency_domain.py @@ -316,7 +316,6 @@ def hmo( segsec: float = 256, bands: dict[str, list[float | None]] = BANDS, fs: float | None = None, - group: str | None = None, **kwargs, ) -> xr.Dataset: """Compute frequency-domain significant wave height (Hmo) per frequency band. @@ -334,7 +333,12 @@ def hmo( fs: Sampling frequency in Hz. If ``None`` (default) it is inferred from the mean sampling interval of ``dim`` (rounded to 3 decimal places), which tolerates small timing jitter. Set explicitly to override. - group: Time component to group by (not typically meaningful for Hmo). + + Note: + Grouping is not supported — a Welch spectrum is not well defined over + the discontinuous record a time grouper produces — so requesting + ``group`` on this stat raises rather than silently returning the + ungrouped result. Returns: Dataset with variables named 'hs_{band_label}', each with a 'band' diff --git a/gridstats/ops/probability.py b/gridstats/ops/probability.py index c867578..7e8dca3 100644 --- a/gridstats/ops/probability.py +++ b/gridstats/ops/probability.py @@ -7,6 +7,7 @@ import dask.array as da import xarray as xr +from gridstats.ops._naming import format_number from gridstats.registry import register_stat logger = logging.getLogger(__name__) @@ -69,7 +70,10 @@ def range_probability( - ``left`` ('closed' | 'open'): Whether the lower bound is inclusive. - ``right`` ('closed' | 'open'): Whether the upper bound is inclusive. - ``label`` (str, optional): Output variable name. Defaults to - ``{var}_{start}_to_{stop}``. + ``{var}_{start}_to_{stop}``, with an absent bound written as + ``min``/``max`` and any decimal point written as ``p`` so the name + stays a valid identifier (a range starting at 1.5 gives + ``hs_1p5_to_max``). Returns: Dataset with one probability variable per range specification, each a @@ -98,8 +102,8 @@ def range_probability( lfunc = left_funcs[left] rfunc = right_funcs[right] - llabel = f"{start:g}" if spec.get("start") is not None else "min" - rlabel = f"{stop:g}" if spec.get("stop") is not None else "max" + llabel = format_number(start) if spec.get("start") is not None else "min" + rlabel = format_number(stop) if spec.get("stop") is not None else "max" varname = spec.get("label", f"{dvar}_{llabel}_to_{rlabel}") # The denominator depends only on the variable, so it is computed once diff --git a/gridstats/ops/rpv.py b/gridstats/ops/rpv.py index 241930f..3d9c325 100644 --- a/gridstats/ops/rpv.py +++ b/gridstats/ops/rpv.py @@ -89,7 +89,6 @@ def rpv( percentile: float = 95, distribution: str = "gumbel_r", duration: float = 24, - group: str | None = None, **kwargs, ) -> xr.Dataset: """Compute return period values using extreme value analysis. @@ -101,8 +100,12 @@ def rpv( percentile: Percentile threshold for peaks-over-threshold selection. distribution: scipy.stats continuous distribution name. duration: Minimum time (hours) between detected peaks. - group: Time component to group by before computing RPV - (e.g. 'month'). Not commonly used for RPV. + + Note: + Grouping is not supported. A per-group extreme-value fit would be made + on a fraction of the samples (a twelfth of them for ``group: month``), + which is not statistically defensible, so requesting ``group`` on this + stat raises rather than silently returning the ungrouped result. Returns: Dataset with a 'period' dimension containing the estimated values. @@ -116,8 +119,8 @@ def rpv( f"Distribution '{distribution}' not found in scipy.stats." ) - dt_ns = float(np.diff(data[dim].values[:2])[0]) - dt_hours = dt_ns / 3.6e12 + # Timestep in hours, via Timedelta so that any datetime64 resolution works. + dt_hours = pd.Timedelta(np.diff(data[dim].values[:2])[0]).total_seconds() / 3600 if dt_hours > duration: raise ValueError( f"Timestep ({dt_hours:.1f}h) is larger than duration ({duration}h). " diff --git a/gridstats/output.py b/gridstats/output.py index 1d75928..780ffa6 100644 --- a/gridstats/output.py +++ b/gridstats/output.py @@ -28,6 +28,41 @@ # Attribute helpers # --------------------------------------------------------------------------- +#: Adjective prefixed to the ``long_name`` of a variable carrying each group +#: suffix. Spelled out rather than derived, because no rule turns all three of +#: 'month', 'season' and 'year' into the word an English reader expects. +_GROUP_ADJECTIVES = {"month": "monthly", "season": "seasonal", "year": "annual"} + +#: Suffix ``Pipeline._apply`` appends for a directionally sectorised call. +_DIRECTIONAL_SUFFIX = "direc" + + +def _split_suffixes(varname: str) -> tuple[str, str | None]: + """Strip the suffixes the pipeline appends from an output variable name. + + ``Pipeline._apply`` builds the output name as + ``{var}_{stat}[_{group}][_direc]``, so both extra tokens have to come off + before the remainder can be read as ``{var}_{stat}``. Only a *trailing* + token counts: a parent variable whose own name contains a group word (say + ``seasonal_hs``) is left alone. + + Args: + varname: Output variable name. + + Returns: + Tuple of the name with any group and directional suffix removed, and + the group token that was removed (None when there was none). + """ + base, sep, last = varname.rpartition("_") + if sep and last == _DIRECTIONAL_SUFFIX: + varname = base + + base, sep, last = varname.rpartition("_") + if sep and last in _GROUP_ADJECTIVES: + return base, last + return varname, None + + def _stat_units(varname: str, stat_units: dict[str, str]) -> str | None: """Resolve the units a statistic imposes on an output variable. @@ -56,8 +91,15 @@ def set_variable_attributes(dsout: xr.Dataset, extra_metadata: dict[str, Any] = and 'stat_units' keys. Statistics listed under 'stat_units' produce a quantity whose units differ - from those of the parent variable (a probability, a percentage); their - declared units replace whatever units the variable would otherwise inherit. + from those of the parent variable (a probability, a percentage, a count); + their declared units replace whatever units the variable would otherwise + inherit. + + The group and directional suffixes the pipeline appends to an output name + (``hs_mean_month``, ``hs_mean_month_direc``) are stripped before the + statistic and parent variable are looked up, and a grouped variable has its + 'long_name' prefixed with the matching adjective ("monthly mean significant + wave height"). Args: dsout: Dataset to annotate. @@ -84,9 +126,12 @@ def set_variable_attributes(dsout: xr.Dataset, extra_metadata: dict[str, Any] = if (value := entry.get(key)) is not None } - # Data variables — variable name is expected to be "{base_var}_{stat}" + # Data variables — variable name is expected to be "{base_var}_{stat}", + # optionally followed by the group and directional suffixes the pipeline + # appends. Those are stripped first so the lookup still resolves. for v, dvar in dsout.data_vars.items(): - parts = v.split("_") + base, group = _split_suffixes(v) + parts = base.split("_") var_name = "_".join(parts[:-1]) stat_name = parts[-1] units = _stat_units(v, stat_units) @@ -97,10 +142,8 @@ def set_variable_attributes(dsout: xr.Dataset, extra_metadata: dict[str, Any] = lng = attrs.get("long_name", std.replace("_", " ")) attrs["standard_name"] = f"{std}_{stat}" attrs["long_name"] = f"{stat} {lng}" - for group_tag in ("month", "year", "season"): - if group_tag in v: - attrs["long_name"] = f"{group_tag}ly " + attrs["long_name"] - break + if group is not None: + attrs["long_name"] = f"{_GROUP_ADJECTIVES[group]} {attrs['long_name']}" dvar.attrs = attrs except KeyError: logger.debug("No metadata found for variable '%s' in attributes.yml", v) diff --git a/gridstats/pipeline.py b/gridstats/pipeline.py index 83201f0..dbadfd6 100644 --- a/gridstats/pipeline.py +++ b/gridstats/pipeline.py @@ -17,10 +17,22 @@ from gridstats.config import CallConfig, PipelineConfig from gridstats.loaders.xarray import _ds_summary from gridstats.output import finalise, upload, write -from gridstats.registry import get_loader, get_stat +from gridstats.registry import get_loader, get_stat, list_group_stats, supports_group logger = logging.getLogger(__name__) +#: Stats that bin the data into directional sectors themselves and name their +#: own outputs after the statistic that produced each one (``statdir`` returns +#: ``hs_mean``, ``hs_max``). ``_apply`` treats them specially in two ways: the +#: sector settings are forwarded to the stat rather than used to wrap it in +#: :func:`_apply_directional`, which would bin the data twice; and the default +#: suffix leaves out the ``_{func}`` token, which would otherwise append the +#: name of the delegating stat to names that already carry the name of the +#: delegated one (``hs_mean_statdir_direc``). The result is that a ``statdir`` +#: call writes the same variable names as the equivalent per-function calls +#: with ``nsector`` set: ``hs_mean_direc``, ``hs_max_direc``. +_SELF_DIRECTIONAL_STATS = frozenset({"statdir"}) + # --------------------------------------------------------------------------- # Dask cluster helpers @@ -220,7 +232,7 @@ def run(self) -> xr.Dataset: dsout = finalise( dsout, source_ds, - chunks=self.config.output.__dict__.get("chunks", {}), + chunks=self.config.output.chunks, metadata=self.config.metadata, global_attrs=self.config.output.global_attrs, mask_config=self.config.output.mask, @@ -266,6 +278,16 @@ def _apply(self, call: CallConfig) -> xr.Dataset: """Run a single stat call and return the renamed result dataset.""" fn = get_stat(call.func) + # Reject grouping up front, before any data is loaded or computed, so a + # stat that cannot honour `group` fails loudly instead of returning the + # ungrouped result under a name the pipeline suffixes with the group. + if call.group and not supports_group(call.func): + raise ValueError( + f"Stat '{call.func}' does not support grouping, but " + f"group='{call.group}' was requested. Stats that support " + f"grouping: {', '.join(list_group_stats())}." + ) + # Load with source-level chunks only (native zarr chunks), then select # the required variables before applying call-level rechunking. This # avoids creating rechunk tasks for variables that are not needed by @@ -304,6 +326,14 @@ def _apply(self, call: CallConfig) -> xr.Dataset: **call.extra_kwargs(), ) + # Stats that sectorise internally take the sector settings as their own + # arguments; the pipeline-level sectoriser is skipped for them below. + self_directional = call.func in _SELF_DIRECTIONAL_STATS + if self_directional: + fn_kwargs["dir_var"] = call.dir_var + if call.nsector is not None: + fn_kwargs["nsector"] = call.nsector + # --- Dask cluster --- cluster_ctx = ( _local_cluster( @@ -319,7 +349,7 @@ def _apply(self, call: CallConfig) -> xr.Dataset: with cluster_ctx, xr.set_options(use_flox=call.use_flox): # --- Directional sectorisation --- - if call.nsector is not None: + if call.nsector is not None and not self_directional: result = _apply_directional( lambda d, **kw: fn(d, **kw), data, @@ -340,10 +370,13 @@ def _apply(self, call: CallConfig) -> xr.Dataset: # --- Rename output variables with suffix --- suffix = call.suffix if suffix is None: - suffix = f"_{call.func}" + # A self-directional stat has already named each output after the + # statistic that produced it, so only the group and sector tokens + # are still missing. + suffix = "" if self_directional else f"_{call.func}" if call.group: suffix += f"_{call.group}" - if call.nsector is not None: + if call.nsector is not None or self_directional: suffix += "_direc" return result.rename({v: f"{v}{suffix}" for v in result.data_vars}) diff --git a/gridstats/registry.py b/gridstats/registry.py index 4e67473..9d6d1c1 100644 --- a/gridstats/registry.py +++ b/gridstats/registry.py @@ -15,6 +15,7 @@ from __future__ import annotations import functools +import inspect import logging from typing import Callable @@ -143,6 +144,37 @@ def list_stats() -> list[str]: return sorted(_STATS) +def supports_group(name: str) -> bool: + """Return whether the named stat can honour a ``group`` argument. + + Support is determined by inspecting the function signature rather than by a + hand-maintained list, so it cannot drift out of sync with the code: a stat + supports grouping when it declares an explicit ``group`` parameter. A stat + that merely absorbs ``group`` into ``**kwargs`` does not, and neither does + one that names its variadic keyword parameter ``group``. + + ``register_stat`` wraps functions with ``functools.wraps``, so + ``inspect.signature`` follows ``__wrapped__`` back to the original function + and sees its real parameters rather than the wrapper's ``*args, **kwargs``. + + Args: + name: Registered stat name. + + Returns: + True if the stat declares an explicit ``group`` parameter. + + Raises: + KeyError: If no stat is registered under that name. + """ + param = inspect.signature(get_stat(name)).parameters.get("group") + return param is not None and param.kind is not inspect.Parameter.VAR_KEYWORD + + +def list_group_stats() -> list[str]: + """Return sorted names of the registered stats that support grouping.""" + return sorted(name for name in _STATS if supports_group(name)) + + def list_loaders() -> list[str]: """Return sorted names of all registered loaders.""" return sorted(_LOADERS) diff --git a/tests/test_config.py b/tests/test_config.py index b70d1ef..49479a5 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -325,6 +325,63 @@ def test_mask_via_pipeline_yaml(self, tmp_path): assert config.output.mask.isel == {"time": 0} +class TestOutputConfig: + def test_defaults(self): + c = OutputConfig(outfile="out.zarr") + assert c.chunks == {} + assert c.updir is None + assert c.append is False + assert c.consolidate is False + assert c.global_attrs == {} + + def test_chunks(self): + c = OutputConfig(outfile="out.zarr", chunks={"latitude": 100, "longitude": 50}) + assert c.chunks == {"latitude": 100, "longitude": 50} + + def test_chunks_from_yaml(self, tmp_path): + yaml_content = textwrap.dedent("""\ + source: + type: xarray + urlpath: data.zarr + output: + outfile: out.zarr + chunks: + latitude: 100 + longitude: 50 + calls: + - func: mean + """) + config_file = tmp_path / "config.yml" + config_file.write_text(yaml_content) + config = PipelineConfig.from_yaml(config_file) + assert config.output.chunks == {"latitude": 100, "longitude": 50} + + def test_chunks_must_be_integers(self): + with pytest.raises(ValidationError): + OutputConfig(outfile="out.zarr", chunks={"latitude": "big"}) + + def test_unknown_field_raises(self): + """A misspelled option must fail rather than be silently dropped.""" + with pytest.raises(ValidationError, match="chunk"): + OutputConfig(outfile="out.zarr", chunk={"latitude": 100}) + + def test_unknown_field_from_yaml_raises(self, tmp_path): + yaml_content = textwrap.dedent("""\ + source: + type: xarray + urlpath: data.zarr + output: + outfile: out.zarr + consolidated: true + calls: + - func: mean + """) + config_file = tmp_path / "config.yml" + config_file.write_text(yaml_content) + with pytest.raises(ValidationError, match="consolidated"): + PipelineConfig.from_yaml(config_file) + + class TestDerivedVarConfig: def test_shorthand_string(self): from gridstats.config import DerivedVarConfig diff --git a/tests/test_ops.py b/tests/test_ops.py index 97a22a7..624d401 100644 --- a/tests/test_ops.py +++ b/tests/test_ops.py @@ -95,6 +95,90 @@ def test_pcount_with_nans(self, ds): out = pcount(ds_nan[["hs"]], dim="time") assert float(out["hs"].max()) < 100.0 + def test_pcount_ungrouped_is_unchanged(self, ds): + from gridstats.ops.aggregations import pcount + ds_nan = ds[["hs"]].copy(deep=True) + ds_nan["hs"].values[: 24 * 10] = np.nan + out = pcount(ds_nan, dim="time") + assert "month" not in out.dims + expected = 100 * (ds.sizes["time"] - 24 * 10) / ds.sizes["time"] + assert float(out["hs"].min()) == pytest.approx(expected) + + def test_pcount_group_adds_dimension(self, ds): + from gridstats.ops.aggregations import pcount + out = pcount(ds[["hs"]], dim="time", group="month") + assert "month" in out.dims + assert "time" not in out.dims + assert list(out.month.values) == list(range(1, 13)) + + def test_pcount_group_complete_series_is_100_everywhere(self, ds): + """A gap-free variable must read 100 % in every group. + + The denominator is the group's own number of timesteps, so this fails if + the whole-record count is used instead. + """ + from gridstats.ops.aggregations import pcount + out = pcount(ds[["hs"]], dim="time", group="month") + np.testing.assert_allclose(out["hs"].values, 100.0) + + def test_pcount_group_denominator_is_per_group(self, ds): + """Gaps in one month must not be smeared across the others. + + January is blanked entirely and February half-blanked, which a + whole-record denominator cannot reproduce: it would report the same + ~84 % everywhere rather than 0 / 50 / 100. + """ + from gridstats.ops.aggregations import pcount + ds_nan = ds[["hs"]].copy(deep=True) + month = ds_nan.time.dt.month.values + + jan = np.where(month == 1)[0] + ds_nan["hs"].values[jan] = np.nan + + feb = np.where(month == 2)[0] + ds_nan["hs"].values[feb[: feb.size // 2]] = np.nan + + out = pcount(ds_nan, dim="time", group="month") + assert float(out["hs"].sel(month=1).max()) == pytest.approx(0.0) + assert float(out["hs"].sel(month=2).min()) == pytest.approx(50.0) + for m in range(3, 13): + assert float(out["hs"].sel(month=m).min()) == pytest.approx(100.0) + + # The whole-record denominator would have given every month this value. + wrong = float(pcount(ds_nan, dim="time")["hs"].min()) + assert wrong == pytest.approx( + 100 + * (ds_nan.sizes["time"] - jan.size - feb.size // 2) + / ds_nan.sizes["time"] + ) + for right in (0.0, 50.0, 100.0): + assert wrong != pytest.approx(right) + + def test_pcount_group_matches_ungrouped_on_each_subset(self, ds): + from gridstats.ops.aggregations import pcount + ds_nan = ds[["hs"]].copy(deep=True) + rng = np.random.default_rng(7) + mask = rng.random(ds_nan.sizes["time"]) < 0.3 + ds_nan["hs"].values[mask] = np.nan + + out = pcount(ds_nan, dim="time", group="month") + for m in (1, 6, 12): + subset = ds_nan.sel(time=ds_nan.time.dt.month == m) + expected = pcount(subset, dim="time") + xr.testing.assert_allclose( + out["hs"].sel(month=m).drop_vars("month"), expected["hs"] + ) + + def test_pcount_group_keeps_percentage_attributes(self, ds): + from gridstats.ops.aggregations import pcount + ds = ds.copy() + ds["hs"].attrs = {"long_name": "significant wave height", "units": "m"} + out = pcount(ds[["hs"]], dim="time", group="month") + assert out["hs"].attrs["units"] == "%" + assert out["hs"].attrs["long_name"] == ( + "percentage of valid significant wave height samples" + ) + def test_pcount_does_not_inherit_parent_units(self, ds): from gridstats.ops.aggregations import pcount ds = ds.copy() @@ -320,6 +404,22 @@ def test_rpv_values_increasing(self, ds): v100 = float(out["hs"].sel(period=100).mean()) assert v100 >= v10 + def test_rpv_microsecond_time_resolution(self, ds): + """The timestep must be derived without assuming nanosecond datetime64.""" + from gridstats.ops.rpv import rpv + ds_us = ds[["hs"]].assign_coords( + time=ds.time.values.astype("datetime64[us]") + ) + assert ds_us.time.dtype == np.dtype("datetime64[us]") + out = rpv(ds_us, dim="time", return_periods=[10, 100]) + expected = rpv(ds[["hs"]], dim="time", return_periods=[10, 100]) + xr.testing.assert_allclose(out["hs"], expected["hs"]) + + def test_rpv_timestep_larger_than_duration_raises(self, ds): + from gridstats.ops.rpv import rpv + with pytest.raises(ValueError, match="larger than duration"): + rpv(ds[["hs"]], dim="time", duration=0.5) + # --------------------------------------------------------------------------- # Distribution @@ -382,6 +482,78 @@ def test_statdir_missing_dir_var_raises(self, ds): with pytest.raises(ValueError, match="not found"): statdir(ds, funcs=["mean"], dir_var="nonexistent", dim="time") + def test_statdir_single_func_names_and_values(self, ds): + """One function: every variable is suffixed with it, values unchanged.""" + from gridstats.ops.directional import statdir + out = statdir(ds, funcs=["mean"], dir_var="dpm", nsector=4, dim="time") + assert set(out.data_vars) == {"hs_mean", "tp_mean", "dpm_mean"} + # Sector centred on 0 wraps across the 0/360 boundary. + mask = (ds.dpm >= 315.0) | (ds.dpm < 45.0) + expected = ds.hs.where(mask).mean(dim="time") + np.testing.assert_allclose( + out.hs_mean.sel(direction=0).values, expected.values, rtol=1e-6 + ) + + def test_statdir_multiple_funcs_do_not_collide(self, ds): + """Two functions used to raise MergeError on the shared base name.""" + from gridstats.ops.directional import statdir + out = statdir(ds, funcs=["mean", "max"], dir_var="dpm", nsector=4, dim="time") + assert set(out.data_vars) == { + "hs_mean", "tp_mean", "dpm_mean", "hs_max", "tp_max", "dpm_max", + } + assert out.direction.size == 4 + assert (out.hs_max >= out.hs_mean).all() + + def test_statdir_three_funcs(self, ds): + from gridstats.ops.directional import statdir + out = statdir( + ds[["hs", "dpm"]], funcs=["mean", "max", "min"], dir_var="dpm", + nsector=4, dim="time", + ) + assert {"hs_mean", "hs_max", "hs_min"} <= set(out.data_vars) + assert len(out.data_vars) == 6 + + def test_statdir_matches_the_delegated_stat_per_sector(self, ds): + """Delegating is equivalent to calling each stat on the masked data.""" + from gridstats.ops.aggregations import max as _max + from gridstats.ops.directional import statdir + out = statdir(ds[["hs", "dpm"]], funcs=["mean", "max"], dir_var="dpm", + nsector=4, dim="time") + mask = (ds.dpm >= 45.0) & (ds.dpm < 135.0) + expected = _max(ds[["hs", "dpm"]].where(mask), dim="time") + np.testing.assert_allclose( + out.hs_max.sel(direction=90).values, expected.hs.values, rtol=1e-6 + ) + + def test_statdir_group_is_passed_through(self, ds): + from gridstats.ops.directional import statdir + out = statdir( + ds[["hs", "dpm"]], funcs=["mean", "max"], dir_var="dpm", nsector=4, + dim="time", group="month", + ) + assert out.month.size == 12 + assert out.direction.size == 4 + assert set(out.hs_mean.dims) == {"month", "direction", "latitude", "longitude"} + + def test_statdir_group_on_unsupporting_func_raises(self, ds): + from gridstats.ops.directional import statdir + with pytest.raises(ValueError, match="do not support grouping"): + statdir( + ds[["hs", "dpm"]], funcs=["mean", "rpv"], dir_var="dpm", + dim="time", group="month", + ) + + def test_statdir_supports_group(self): + """statdir declares `group`, so the pipeline guard must let it through.""" + from gridstats.registry import list_group_stats, supports_group + assert supports_group("statdir") + assert "statdir" in list_group_stats() + + def test_statdir_empty_funcs_raises(self, ds): + from gridstats.ops.directional import statdir + with pytest.raises(ValueError, match="at least one function"): + statdir(ds, funcs=[], dir_var="dpm", dim="time") + # --------------------------------------------------------------------------- # Probability @@ -527,6 +699,64 @@ def test_range_probability_long_name_describes_bounds(self, ds): assert out["band"].attrs["long_name"] == "probability of hs >= 1 and < 2" assert out["all"].attrs["long_name"] == "probability of hs over the full range" + def test_default_label_fractional_bound_has_no_dot(self, ds): + """A dot makes the name unusable as an identifier or in a path.""" + from gridstats.ops.probability import range_probability + out = range_probability( + ds, + dim="time", + data_ranges=[ + {"var": "hs", "start": 1.5, "stop": None}, + {"var": "hs", "start": None, "stop": 2.5}, + {"var": "hs", "start": 0.25, "stop": 1.5}, + ], + ) + assert "hs_1p5_to_max" in out + assert "hs_min_to_2p5" in out + assert "hs_0p25_to_1p5" in out + assert not any("." in str(v) for v in out.data_vars) + + def test_default_label_integral_bound_unchanged(self, ds): + from gridstats.ops.probability import range_probability + out = range_probability( + ds, + dim="time", + data_ranges=[{"var": "hs", "start": 1.0, "stop": 2.0}], + ) + assert "hs_1_to_2" in out + + def test_explicit_label_passed_through_untouched(self, ds): + """An explicit label is the caller's business, dots and all.""" + from gridstats.ops.probability import range_probability + out = range_probability( + ds, + dim="time", + data_ranges=[{"var": "hs", "start": 1.5, "stop": None, "label": "hs_1.5+"}], + ) + assert "hs_1.5+" in out + + +# --------------------------------------------------------------------------- +# Output variable naming +# --------------------------------------------------------------------------- + +class TestNameFormatting: + """The numeric-bound formatter shared by the naming call sites.""" + + @pytest.mark.parametrize( + "value,expected", + [(1.5, "1p5"), (2.0, "2"), (2, "2"), (0.25, "0p25"), (-1.5, "-1p5")], + ) + def test_format_number(self, value, expected): + from gridstats.ops._naming import format_number + assert format_number(value) == expected + + def test_both_call_sites_share_the_formatter(self): + """`exceedance` and `probability` must not drift apart on naming.""" + from gridstats.ops import _naming, exceedance, probability + assert exceedance.format_number is _naming.format_number + assert probability.format_number is _naming.format_number + # --------------------------------------------------------------------------- # Derived variables diff --git a/tests/test_output.py b/tests/test_output.py index 531eae1..b20c715 100644 --- a/tests/test_output.py +++ b/tests/test_output.py @@ -76,6 +76,15 @@ def _ds(self, *varnames): ds[v].attrs = dict(self.HS_ATTRS) return ds + def test_count_is_dimensionless(self): + """A count of samples is a number, not a length.""" + out = set_variable_attributes(self._ds("hs_count")) + assert out["hs_count"].attrs["units"] == "1" + + def test_grouped_count_is_dimensionless(self): + out = set_variable_attributes(self._ds("hs_count_month")) + assert out["hs_count_month"].attrs["units"] == "1" + def test_pcount_is_a_percentage(self): out = set_variable_attributes(self._ds("hs_pcount")) assert out["hs_pcount"].attrs["units"] == "%" @@ -119,6 +128,129 @@ def test_extra_metadata_can_declare_stat_units(self): assert out["hs_myprob"].attrs["units"] == "1" +class TestGroupedVariableAttributes: + """Attributes of variables carrying the suffixes ``Pipeline._apply`` appends. + + A grouped call is written out as ``{var}_{stat}_{group}``, optionally with a + further ``_direc`` for a directionally sectorised call, so the statistic and + the parent variable can only be resolved once those suffixes come off. + """ + + def _ds(self, *varnames): + return xr.Dataset( + {v: (["x"], np.array([0.5], dtype="float32")) for v in varnames} + ) + + def test_monthly_variable_keeps_parent_units(self): + out = set_variable_attributes(self._ds("hs_mean_month")) + attrs = out["hs_mean_month"].attrs + assert attrs["units"] == "m" + assert attrs["standard_name"] == "sea_surface_wave_significant_height_mean" + assert attrs["long_name"].startswith("monthly ") + assert attrs["long_name"] == "monthly mean significant wave height" + + @pytest.mark.parametrize( + "varname,adjective", + [ + ("hs_mean_month", "monthly"), + ("hs_mean_season", "seasonal"), + ("hs_mean_year", "annual"), + ], + ) + def test_group_adjectives(self, varname, adjective): + out = set_variable_attributes(self._ds(varname)) + assert out[varname].attrs["long_name"].startswith(f"{adjective} ") + + def test_seasonly_is_never_produced(self): + """The adjective is a lookup, not 'month'/'season'/'year' plus 'ly'.""" + out = set_variable_attributes( + self._ds("hs_mean_season", "hs_mean_year", "hs_max_season") + ) + for v in out.data_vars: + assert "seasonly" not in out[v].attrs["long_name"] + assert "yearly" not in out[v].attrs["long_name"] + + def test_directional_suffix_resolves(self): + out = set_variable_attributes(self._ds("hs_mean_direc")) + attrs = out["hs_mean_direc"].attrs + assert attrs["units"] == "m" + assert attrs["long_name"] == "mean significant wave height" + + def test_group_and_directional_suffixes_together(self): + """`_apply` composes the group suffix before the directional one.""" + out = set_variable_attributes(self._ds("hs_mean_month_direc")) + attrs = out["hs_mean_month_direc"].attrs + assert attrs["units"] == "m" + assert attrs["long_name"] == "monthly mean significant wave height" + + def test_grouped_count_overrides_parent_units(self): + """Composes with the stat_units override rather than fighting it.""" + out = set_variable_attributes(self._ds("hs_count_month")) + attrs = out["hs_count_month"].attrs + assert attrs["units"] == "1" + assert attrs["long_name"] == "monthly count significant wave height" + + def test_grouped_probability_stays_dimensionless(self): + """The KeyError fallback path must still apply the stat_units override.""" + ds = self._ds("hs_gt_1p5_range_probability_month") + ds["hs_gt_1p5_range_probability_month"].attrs = { + "standard_name": "sea_surface_wave_significant_height", + "long_name": "significant wave height", + "units": "m", + } + out = set_variable_attributes(ds) + attrs = out["hs_gt_1p5_range_probability_month"].attrs + assert attrs["units"] == "1" + assert "standard_name" not in attrs + assert attrs.get("long_name") != "significant wave height" + + @pytest.mark.parametrize("varname", ["hs_mean", "hs_quantile", "hs_max"]) + def test_ungrouped_variables_unaffected(self, varname): + out = set_variable_attributes(self._ds(varname)) + attrs = out[varname].attrs + assert attrs["units"] == "m" + assert attrs["long_name"] == { + "hs_mean": "mean significant wave height", + "hs_quantile": "quantile significant wave height", + "hs_max": "maximum significant wave height", + }[varname] + + def test_parent_name_containing_a_group_word_is_not_mangled(self): + """Only a *trailing* group token is a group; a substring is not.""" + extra = { + "data_vars": { + "season_length": { + "standard_name": "storm_season_length", + "long_name": "storm season length", + "units": "d", + } + } + } + out = set_variable_attributes(self._ds("season_length_mean"), extra) + attrs = out["season_length_mean"].attrs + assert attrs["units"] == "d" + assert attrs["long_name"] == "mean storm season length" + assert not attrs["long_name"].startswith("seasonal") + + def test_group_word_as_parent_suffix_still_grouped(self): + """`{parent ending in a group word}_{stat}_{group}` resolves both.""" + extra = { + "data_vars": { + "wave_year": { + "standard_name": "wave_year", + "long_name": "wave year", + "units": "1", + } + } + } + out = set_variable_attributes(self._ds("wave_year_mean_year"), extra) + assert out["wave_year_mean_year"].attrs["long_name"] == "annual mean wave year" + + def test_unknown_grouped_variable_does_not_raise(self): + out = set_variable_attributes(self._ds("foo_bar_month")) + assert "foo_bar_month" in out + + class TestSetGlobalAttributes: def test_sets_institution(self, source_ds, dsout): out = set_global_attributes(source_ds, dsout) diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py index d52ccb9..54bb1ca 100644 --- a/tests/test_pipeline.py +++ b/tests/test_pipeline.py @@ -278,3 +278,350 @@ def test_apply_tiled(self, tmp_path, netcdf_source): assert "hs_mean" in dsout # Should cover the same grid as without tiling assert dsout["hs_mean"].latitude.size == 3 + + +class TestGroupSupportGuard: + """`group` on a stat that cannot honour it must raise, not be ignored.""" + + @staticmethod + def _config(tmp_path, netcdf_source, **call_kwargs): + from gridstats.config import ( + CallConfig, OutputConfig, PipelineConfig, XarraySourceConfig, + ) + + return PipelineConfig( + source=XarraySourceConfig( + type="xarray", urlpath=netcdf_source, engine="netcdf4" + ), + output=OutputConfig(outfile=str(tmp_path / "out.nc")), + calls=[CallConfig(**call_kwargs)], + ) + + @pytest.mark.parametrize("func", ["rpv", "hmo", "distribution3_timestep"]) + def test_group_on_unsupporting_stat_raises(self, tmp_path, netcdf_source, func): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, netcdf_source, func=func, dim="time", data_vars=["hs"], + group="month", + ) + with pytest.raises(ValueError, match="does not support grouping") as excinfo: + Pipeline(config).run() + + message = str(excinfo.value) + assert f"'{func}'" in message + assert "group='month'" in message + # The message points at the stats that do support it. + assert "mean" in message + assert "pcount" in message + assert func not in message.split("Stats that support grouping:")[1] + + def test_guard_fires_before_any_data_is_loaded(self, tmp_path): + """The check must not depend on the source being readable.""" + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, str(tmp_path / "does_not_exist.nc"), func="rpv", dim="time", + data_vars=["hs"], group="month", + ) + with pytest.raises(ValueError, match="does not support grouping"): + Pipeline(config).run() + + def test_guard_fires_for_tiled_calls(self, tmp_path, netcdf_source): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, netcdf_source, func="rpv", dim="time", data_vars=["hs"], + group="month", tiles={"latitude": 2}, + ) + with pytest.raises(ValueError, match="does not support grouping"): + Pipeline(config).run() + + def test_guard_fires_for_directional_calls(self, tmp_path, netcdf_source): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, netcdf_source, func="rpv", dim="time", data_vars=["hs"], + group="month", nsector=4, dir_var="dpm", + ) + with pytest.raises(ValueError, match="does not support grouping"): + Pipeline(config).run() + + def test_unsupporting_stat_without_group_is_allowed(self, tmp_path, netcdf_source): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, netcdf_source, func="rpv", dim="time", data_vars=["hs"], + return_periods=[10, 100], + ) + dsout = Pipeline(config).run() + assert "hs_rpv" in dsout + + @pytest.mark.parametrize("func", ["mean", "pcount"]) + def test_group_on_supporting_stat_does_not_raise( + self, tmp_path, netcdf_source, func + ): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, netcdf_source, func=func, dim="time", data_vars=["hs"], + group="month", + ) + dsout = Pipeline(config).run() + assert f"hs_{func}_month" in dsout + assert "month" in dsout[f"hs_{func}_month"].dims + + def test_group_on_range_probability_does_not_raise(self, tmp_path, netcdf_source): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, netcdf_source, func="range_probability", dim="time", + data_vars=["hs"], group="month", + data_ranges=[{"var": "hs", "start": 1.0, "label": "hs_above_1"}], + ) + dsout = Pipeline(config).run() + assert "hs_above_1_range_probability_month" in dsout + assert "month" in dsout["hs_above_1_range_probability_month"].dims + + def test_grouped_pcount_denominator_is_per_group_end_to_end( + self, tmp_path, sample_ds + ): + """The pipeline result must be the grouped value, not the record-wide one.""" + from gridstats.pipeline import Pipeline + + ds = sample_ds.copy(deep=True) + jan = np.where(ds.time.dt.month.values == 1)[0] + ds["hs"].values[jan] = np.nan + path = str(tmp_path / "gappy.nc") + ds.to_netcdf(path) + + config = self._config( + tmp_path, path, func="pcount", dim="time", data_vars=["hs"], group="month", + ) + dsout = Pipeline(config).run() + pc = dsout["hs_pcount_month"] + assert float(pc.sel(month=1).max()) == pytest.approx(0.0) + assert float(pc.sel(month=7).min()) == pytest.approx(100.0) + + +# --------------------------------------------------------------------------- +# Directional sector stats through the pipeline +# --------------------------------------------------------------------------- + +@pytest.fixture +def directional_ds(sample_ds): + """Sample dataset with a direction variable, needed by directional calls.""" + rng = np.random.default_rng(3) + ds = sample_ds.copy(deep=True) + ds["dpm"] = (ds.hs.dims, (rng.random(ds.hs.shape) * 360).astype("float32")) + return ds + + +@pytest.fixture +def directional_source(tmp_path, directional_ds): + path = tmp_path / "source_dir.nc" + directional_ds.to_netcdf(path) + return str(path) + + +class TestStatdirPipeline: + """`statdir` output names must survive `_apply`'s renaming intact.""" + + @staticmethod + def _config(tmp_path, source, outfile="out.nc", **call_kwargs): + from gridstats.config import ( + CallConfig, OutputConfig, PipelineConfig, XarraySourceConfig, + ) + + return PipelineConfig( + source=XarraySourceConfig( + type="xarray", urlpath=source, engine="netcdf4" + ), + output=OutputConfig(outfile=str(tmp_path / outfile)), + calls=[CallConfig(**call_kwargs)], + ) + + def test_multi_func_statdir_names(self, tmp_path, directional_source): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, directional_source, func="statdir", dim="time", + data_vars=["hs", "dpm"], dir_var="dpm", funcs=["mean", "max"], + ) + dsout = Pipeline(config).run() + assert set(dsout.data_vars) == { + "hs_mean_direc", "hs_max_direc", "dpm_mean_direc", "dpm_max_direc", + } + # The delegating stat's own name is not appended on top of the + # delegated one: no 'hs_mean_statdir_direc'. + assert not any("statdir" in v for v in dsout.data_vars) + assert dsout.direction.size == 4 + + def test_single_func_statdir_names(self, tmp_path, directional_source): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, directional_source, func="statdir", dim="time", + data_vars=["hs", "dpm"], dir_var="dpm", funcs=["mean"], + ) + dsout = Pipeline(config).run() + assert set(dsout.data_vars) == {"hs_mean_direc", "dpm_mean_direc"} + + def test_nsector_is_forwarded_not_applied_twice(self, tmp_path, directional_source): + """`nsector` on a statdir call must bin the data once, at the requested width.""" + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, directional_source, func="statdir", dim="time", + data_vars=["hs", "dpm"], dir_var="dpm", nsector=8, funcs=["mean", "max"], + ) + dsout = Pipeline(config).run() + assert dsout.direction.size == 8 + np.testing.assert_allclose( + dsout.direction.values, np.arange(0, 360, 45), atol=1e-6 + ) + + def test_statdir_matches_equivalent_nsector_calls(self, tmp_path, directional_source): + """statdir is a one-pass equivalent of the same stats called with nsector.""" + from gridstats.config import CallConfig + from gridstats.pipeline import Pipeline + + via_statdir = Pipeline( + self._config( + tmp_path, directional_source, outfile="statdir.nc", func="statdir", + dim="time", data_vars=["hs", "dpm"], dir_var="dpm", nsector=8, + funcs=["mean", "max"], + ) + ).run() + + config = self._config( + tmp_path, directional_source, outfile="nsector.nc", func="mean", + dim="time", data_vars=["hs", "dpm"], dir_var="dpm", nsector=8, + ) + config.calls.append( + CallConfig( + func="max", dim="time", data_vars=["hs", "dpm"], dir_var="dpm", + nsector=8, + ) + ) + via_nsector = Pipeline(config).run() + + assert set(via_statdir.data_vars) == set(via_nsector.data_vars) + for var in via_statdir.data_vars: + xr.testing.assert_allclose(via_statdir[var], via_nsector[var]) + + def test_grouped_statdir(self, tmp_path, directional_source): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, directional_source, func="statdir", dim="time", + data_vars=["hs", "dpm"], dir_var="dpm", group="month", + funcs=["mean", "max"], + ) + dsout = Pipeline(config).run() + assert "hs_mean_month_direc" in dsout + assert "hs_max_month_direc" in dsout + assert set(dsout["hs_mean_month_direc"].dims) == { + "month", "direction", "latitude", "longitude" + } + + def test_statdir_output_carries_attributes(self, tmp_path, directional_source): + """The name must still resolve through attributes.yml after renaming.""" + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, directional_source, func="statdir", dim="time", + data_vars=["hs", "dpm"], dir_var="dpm", funcs=["mean", "max"], + ) + dsout = Pipeline(config).run() + attrs = dsout["hs_mean_direc"].attrs + assert attrs["units"] == "m" + assert attrs["long_name"] == "mean significant wave height" + + def test_explicit_suffix_still_wins(self, tmp_path, directional_source): + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, directional_source, func="statdir", dim="time", + data_vars=["hs", "dpm"], dir_var="dpm", funcs=["mean"], suffix="", + ) + dsout = Pipeline(config).run() + assert "hs_mean" in dsout + + +# --------------------------------------------------------------------------- +# Output chunking +# --------------------------------------------------------------------------- + +class TestOutputChunks: + """`output.chunks` must reach `finalise` and change the written store.""" + + @staticmethod + def _config(tmp_path, source, chunks, outfile="out.zarr"): + from gridstats.config import ( + CallConfig, OutputConfig, PipelineConfig, XarraySourceConfig, + ) + + return PipelineConfig( + source=XarraySourceConfig( + type="xarray", urlpath=source, engine="netcdf4" + ), + output=OutputConfig(outfile=str(tmp_path / outfile), chunks=chunks), + calls=[CallConfig(func="mean", dim="time", data_vars=["hs"])], + ) + + def test_chunks_reach_the_written_store(self, tmp_path, netcdf_source): + from gridstats.pipeline import Pipeline + + config = self._config(tmp_path, netcdf_source, {"latitude": 1, "longitude": 2}) + dsout = Pipeline(config).run() + assert dsout["hs_mean"].chunks == ((1, 1, 1), (2, 1)) + + written = xr.open_zarr(str(tmp_path / "out.zarr")) + assert written["hs_mean"].encoding["chunks"] == (1, 2) + + def test_without_chunks_the_output_is_a_single_chunk(self, tmp_path, netcdf_source): + """Baseline for the test above: the grid is 3x3 unless chunks say otherwise.""" + from gridstats.pipeline import Pipeline + + config = self._config(tmp_path, netcdf_source, {}, outfile="unchunked.zarr") + Pipeline(config).run() + + written = xr.open_zarr(str(tmp_path / "unchunked.zarr")) + assert written["hs_mean"].encoding["chunks"] == (3, 3) + + def test_chunks_from_yaml_reach_the_written_store(self, tmp_path, netcdf_source): + from gridstats.pipeline import Pipeline + + yaml_content = textwrap.dedent(f"""\ + source: + type: xarray + urlpath: {netcdf_source} + engine: netcdf4 + output: + outfile: {tmp_path}/from_yaml.zarr + chunks: + latitude: 2 + longitude: 1 + calls: + - func: mean + dim: time + data_vars: [hs] + """) + cfg_file = tmp_path / "config.yml" + cfg_file.write_text(yaml_content) + Pipeline.from_yaml(cfg_file).run() + + written = xr.open_zarr(str(tmp_path / "from_yaml.zarr")) + assert written["hs_mean"].encoding["chunks"] == (2, 1) + + def test_chunks_for_absent_dims_are_ignored(self, tmp_path, netcdf_source): + """`finalise` only applies chunks for dims present in the output.""" + from gridstats.pipeline import Pipeline + + config = self._config( + tmp_path, netcdf_source, {"time": 10, "latitude": 1}, outfile="partial.zarr" + ) + dsout = Pipeline(config).run() + assert "time" not in dsout.dims + assert dsout["hs_mean"].chunks == ((1, 1, 1), (3,)) diff --git a/tests/test_registry.py b/tests/test_registry.py index 37d9104..631feb8 100644 --- a/tests/test_registry.py +++ b/tests/test_registry.py @@ -7,10 +7,12 @@ _STATS, get_loader, get_stat, + list_group_stats, list_loaders, list_stats, register_loader, register_stat, + supports_group, ) @@ -75,6 +77,67 @@ def a(data, **kwargs): assert names == sorted(names) +class TestSupportsGroup: + """Group support is read off the signature, so it cannot drift from the code.""" + + def test_explicit_group_parameter_supports_group(self): + @register_stat("with_group") + def fn(data, *, dim="time", group=None, **kwargs): + return data + + assert supports_group("with_group") is True + + def test_group_absorbed_by_kwargs_does_not_support_group(self): + @register_stat("without_group") + def fn(data, *, dim="time", **kwargs): + return data + + assert supports_group("without_group") is False + + def test_var_keyword_named_group_does_not_support_group(self): + @register_stat("var_keyword_group") + def fn(data, *, dim="time", **group): + return data + + assert supports_group("var_keyword_group") is False + + def test_signature_is_read_through_the_registration_wrapper(self): + """``functools.wraps`` must expose the wrapped function's parameters.""" + import inspect + + @register_stat("wrapped_probe") + def fn(data, *, dim="time", group=None, **kwargs): + return data + + params = inspect.signature(get_stat("wrapped_probe")).parameters + assert "group" in params + assert list(params) != ["data", "args", "kwargs"] + + def test_builtin_supporting_and_non_supporting_stats(self): + import gridstats.ops # noqa: F401 (trigger registration) + + assert supports_group("mean") is True + assert supports_group("range_probability") is True + assert supports_group("pcount") is True + assert supports_group("rpv") is False + assert supports_group("hmo") is False + assert supports_group("distribution3_timestep") is False + + def test_unknown_stat_raises_key_error(self): + with pytest.raises(KeyError, match="not registered"): + supports_group("nonexistent_xyz") + + def test_list_group_stats_sorted_and_filtered(self): + import gridstats.ops # noqa: F401 (trigger registration) + + names = list_group_stats() + assert names == sorted(names) + assert "mean" in names + assert "pcount" in names + assert "rpv" not in names + assert set(names) <= set(list_stats()) + + class TestRegisterLoader: def test_register_and_get(self): @register_loader("test_loader")