rafa-guedes · Aug 12, 2026 · Aug 12, 2026 · Aug 12, 2026 · Aug 12, 2026
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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")