Add ensemble functionality

[KristianHMoller]

The underlying probabilistic advection method has the option of running an ensemble with noise. However, until now, this option has not been implemented into sunflow.

This PR allows for running as an ensemble by specifying the CLI argument --ensemble_members to be different from one. Per default, the output will be selected ensemble statistics, but specifying --full_ensemble when calling sunflow will instead provide the full ensemble as output. The output_mode global attribute in the NetCDF file will specify the type of output: full ensemble, ensemble_statistic or deterministic.

If only a subset of the ensemble statistics are desired, this can be specified by exporting the environment variable ENSEMBLE_STATISTICS with the desired variables comma separated.

The default values of the noise parameters alpha and beta are 0 for deterministic mode and using other values will result in a warning.
For ensemble mode, currently, the default values are taken from Carpentieri et al.: https://www.sciencedirect.com/science/article/pii/S030626192301139X

[KristianHMoller]

Scientific question 1: Which values to use for alpha and beta?
Currently using optimal values from Carpentieri et al., but these may be country-specific and are presumably derived with the potential error in the implementation (dmidk/SolarSTEPS#2).
Also: @irenelivia , we used slightly different values in our previous code. Do you recall why?

[KristianHMoller]

Implementation question 1: Do we want warning or failure if running deterministically (only one ensemble member) with noise?

[KristianHMoller]

Implementation question 2: Do we want the number of ensemble members and noise parameters to be CLI arguments rather than environemnt variables? I.e. parameters called when running sunflow: sunflow --ens_members 10?

Current decision: CLI argument

[KristianHMoller]

Ideas for ensemble spread variables: Quantiles, 10/90 percentile
Implemented as default options: 10, 25, 75 and 90th percentiles

[KristianHMoller]

Mean vs. median - which is more useful?
Current decision: Default uses both

[KristianHMoller]

@irenelivia If you want, you can try this out to see if anything obvious stands out as ideas for improvements or things that should be changed.

Ideas for simple commands to test:
source .venv/bin/activate

sunflow --ensemble_members 10

sunflow --ensemble_members 10 --full_ensemble

export ENSEMBLE_STATISTICS=mean,median
sunflow --ensemble_members 10

[irenelivia]

This might be a little bug that my friend Claude helped me find!
With this line, we are checking/validating the clearsky field only for the "previous_day_time_steps" = the forecast time steps, but not for clearsky_t0_time. This might lead to an unexplained crash later on in forecast.py, if clearsky data for clearsky_t0_time is missing, instead of failing at the "validation step"

Here's a walk-through of what happens:

# Line 352-353: t=0 time is created
clearsky_t0_time = time_step - timedelta(days=1)
all_clearsky_time_steps = [clearsky_t0_time] + previous_day_time_steps

# Line 357-367: t=0 data is fetched
clearsky_data = fetch_clearsky_with_fallback(all_clearsky_time_steps, ...)

# Line 378: BUG - only validates forecast steps, not t=0!
validate_clearsky_completeness(clearsky_data, previous_day_time_steps)

# Line 400: prepend_t0() uses clearsky_t0_time
solar_forecast = prepend_t0(clearsky_data, ratio_data, ..., clearsky_t0_time)

# forecast.py:205: CRASH - KeyError if missing
sds_cs_t0 = clearsky_data.sel(time=clearsky_t0_time.replace(tzinfo=None))[...].values

[irenelivia]

Fix:

# OLD: Only validates forecast steps
validate_clearsky_completeness(clearsky_data, previous_day_time_steps)

# NEW: Validates all times including t=0
validate_clearsky_completeness(clearsky_data, all_clearsky_time_steps)

[irenelivia]

Check my comment on Line 378 of main.py
This will fail if clearsky_data for t0 is missing

[irenelivia]

Suggestion to make sure that the shape of the forecast matches the ensemble members chosen in the config:

if forecast.ndim == 4:
assert forecast.shape[0] == nowcast_config.ens_members, \
      f"Ensemble mismatch: array has {forecast.shape[0]} but config has {nowcast_config.ens_members}"

[irenelivia]

Suggestion to update this docstring to

Returns:
    Forecast array of shape (n_steps, lat, lon) for single-member runs,
    or (ens_members, n_steps, lat, lon) for multi-member ensemble runs.

[irenelivia]

Suggestion to put {statistic} first in the name, so it reads:

mean_GHI_probabilistic_advection
median_GHI_probabilistic_advection
and so on...

This makes it easier when making a fast visualization with ncview 😄 Otherwise they're all called GHI_probabilistic_advection in the viewer

[irenelivia]

Awesome work @KristianHMoller .

Three improvements suggested to strengthen this PR:

(1) Little validation bug in main.py:378 - the clearsky completeness check should include the t=0 timestamp (all_clearsky_time_steps instead of previous_day_time_steps), preventing crashes when t=0 data is missing.

(2) Clarify docstring in forecast.py:89-90 to document that return shape varies by ens_members (single-member: 3D, multi-member: 4D), eliminating confusion for future users.

(3) Add dimension check in data_io.py to validate that ensemble array dimensions match the configured ens_members count, catching programmer errors early if validation logic changes.

I tested:

1. sunflow --ensemble_members 10 --time 2026-07-17T09:15Z --full_ensemble

Works seamlessly with the ensemble dimension. Nice! Very interesting to skim through the different noise images.
(Note: this still has the "arctan artifacts")

2. sunflow --ensemble_members 10 --time 2026-07-17T09:15Z

Also works fine, just a small comment to change the output names of the statistical variables, so the statistic comes first in the name (easier to see on ncview for oldies like me 😄)
(Note: this still has the "arctan artifacts")

3. export ENSEMBLE_STATISTICS=mean,median
   sunflow --ensemble_members 10 --time 2026-07-17T09:15Z

Also works fine, outputting only mean and median statistics.
(Note: this still has the "arctan artifacts")

4. sunflow --ensemble_members 0 --time 2026-07-17T09:15Z

Works fine, outputs 1 ensemble member, no arctan artifacts

6. sunflow --ensemble_members 1 --time 2026-07-17T09:15Z

Works fine, outputs 1 ensemble member, no arctan artifacts

[KristianHMoller]

Mean and percentiles do not seem to cover full domain, while median does. Investigate why!

[KristianHMoller]

Significant speed-up should be possible by doing cropping to domain_nowcast directly following the forecast calculation