From 80b8c3d6f0243d0d01c1445e543ed9bfd51d6fb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Wed, 1 Jul 2026 12:00:17 +0200 Subject: [PATCH 01/13] First ensemble functionality --- sunflow/forecast.py | 70 +++++++++++++++++++++++++++++++++++++-------- sunflow/main.py | 11 ++----- 2 files changed, 61 insertions(+), 20 deletions(-) diff --git a/sunflow/forecast.py b/sunflow/forecast.py index 3abe325..d55a7ce 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -4,6 +4,7 @@ import numpy as np import xarray as xr from Models.ProbabilisticAdvection import ProbabilisticAdvection +from loguru import logger from .geospatial import get_coordinates @@ -93,10 +94,6 @@ def simple_advection_forecast( # Run probabilistic advection using the correct method name forecast = pa.maps_forecast(n_steps, ratio_data, motion_field) - # Remove ensemble dimension if present (squeeze to get shape: [time, lat, lon]) - if forecast.ndim == 4: # [ensemble, time, lat, lon] - forecast = forecast[0] # Take first (and only) ensemble member - return forecast @@ -115,7 +112,7 @@ def multiply_clearsky( Args: ratio_forecast: Forecast array of shape (n_steps, lat, lon) - containing SDS/SDS_CS ratios. + or (ensemble, n_steps, lat, lon) containing SDS/SDS_CS ratios. clearsky_data: xarray Dataset with a 'time' dimension containing the clearsky variable for each forecast step. previous_day_time_steps: List of datetimes (one per forecast step) @@ -124,26 +121,75 @@ def multiply_clearsky( NetCDF variable name in the datasets. Returns: - Solar irradiance forecast array of shape (n_steps, lat, lon) - in W m⁻². + Solar irradiance forecast array with the same shape as + ratio_forecast, in W m⁻². Raises: RuntimeError: If clearsky data is missing for any forecast timestep. """ - solar_forecast = np.zeros_like(ratio_forecast) + clearsky_steps: list[np.ndarray] = [] - for i, time_step in enumerate(previous_day_time_steps): + for time_step in previous_day_time_steps: try: sds_cs = clearsky_data.sel(time=time_step.replace(tzinfo=None))[ nc_variable_names["sds_cs"] ].values - - # Multiply ratio by clearsky - solar_forecast[i] = ratio_forecast[i] * sds_cs + clearsky_steps.append(sds_cs) except KeyError: raise RuntimeError( f"No clearsky data for {time_step.strftime('%Y-%m-%dT%H:%M:%SZ')}, " "cannot compute solar forecast for this step." ) + clearsky_stack = np.stack(clearsky_steps, axis=0) + + if ratio_forecast.ndim == 3: + if ratio_forecast.shape[0] != clearsky_stack.shape[0]: + raise ValueError( + "ratio_forecast time dimension does not match clearsky timesteps " + f"({ratio_forecast.shape[0]} != {clearsky_stack.shape[0]})." + ) + return ratio_forecast * clearsky_stack + + if ratio_forecast.ndim == 4: + if ratio_forecast.shape[1] != clearsky_stack.shape[0]: + raise ValueError( + "ratio_forecast time dimension does not match clearsky timesteps " + f"({ratio_forecast.shape[1]} != {clearsky_stack.shape[0]})." + ) + return ratio_forecast * clearsky_stack[np.newaxis, :, :, :] + + raise ValueError( + "ratio_forecast must have shape (time, lat, lon) or " + "(ensemble, time, lat, lon)." + ) + + +def prepend_t0(clearsky_data: xr.Dataset, ratio_data: np.ndarray, solar_forecast: np.ndarray, config: dict, clearsky_t0_time: datetime) -> np.ndarray: + # Prepend timestep 0: current observation (ratio_data[-1]) × clearsky at t=0 + sds_cs_t0 = clearsky_data.sel(time=clearsky_t0_time.replace(tzinfo=None))[ + config["nc_variable_names"]["sds_cs"] + ].values + solar_t0 = ratio_data[-1] * sds_cs_t0 + + if solar_forecast.ndim == 3: + solar_forecast = np.concatenate( + [solar_t0[np.newaxis, :, :], solar_forecast], + axis=0, + ) + elif solar_forecast.ndim == 4: + # Broadcast the same t=0 clearsky-based analysis field to all ensembles. + solar_t0_ens = np.broadcast_to( + solar_t0, + (solar_forecast.shape[0],) + solar_t0.shape, + ) + solar_forecast = np.concatenate( + [solar_t0_ens[:, np.newaxis, :, :], solar_forecast], + axis=1, + ) + else: + raise ValueError( + "solar_forecast must have shape (time, lat, lon) or " + "(ensemble, time, lat, lon)." + ) return solar_forecast diff --git a/sunflow/main.py b/sunflow/main.py index 6034bab..3c57aad 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -22,7 +22,7 @@ save_forecast, ) from .downloaders import download_past_data -from .forecast import multiply_clearsky, preprocess_data, simple_advection_forecast +from .forecast import multiply_clearsky, preprocess_data, simple_advection_forecast, prepend_t0 from .geospatial import check_solar_elevation, get_bbox from .time_handler import generate_time_steps, round_time from .validation import ( @@ -331,12 +331,7 @@ def run_nowcast( config["nc_variable_names"], ) - # Prepend timestep 0: current observation (ratio_data[-1]) × clearsky at t=0 - sds_cs_t0 = clearsky_data.sel(time=clearsky_t0_time.replace(tzinfo=None))[ - config["nc_variable_names"]["sds_cs"] - ].values - solar_t0 = ratio_data[-1] * sds_cs_t0 - solar_forecast = np.concatenate([solar_t0[np.newaxis, :, :], solar_forecast], axis=0) + solar_forecast = prepend_t0(clearsky_data, ratio_data, solar_forecast, config, clearsky_t0_time) # Save forecast (now contains actual solar irradiance, not ratios) filename = save_forecast( @@ -413,7 +408,7 @@ def cli() -> None: validate_run_mode(run_mode, dataset_name) validate_config(config, dataset_name) - validate_nowcast_config(nowcast_config) + # validate_nowcast_config(nowcast_config) verify_environment_variables(run_mode, dataset_name) # Determine the time steps to run From 5406030ee2068fb9355cb876a9ecb8b0b0528949 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Wed, 1 Jul 2026 12:07:47 +0200 Subject: [PATCH 02/13] Remove forcing of only a single ensemble member --- sunflow/main.py | 1 - sunflow/validation.py | 21 --------------------- 2 files changed, 22 deletions(-) diff --git a/sunflow/main.py b/sunflow/main.py index 3c57aad..40356ca 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -408,7 +408,6 @@ def cli() -> None: validate_run_mode(run_mode, dataset_name) validate_config(config, dataset_name) - # validate_nowcast_config(nowcast_config) verify_environment_variables(run_mode, dataset_name) # Determine the time steps to run diff --git a/sunflow/validation.py b/sunflow/validation.py index 21a1cae..ff1680a 100644 --- a/sunflow/validation.py +++ b/sunflow/validation.py @@ -46,27 +46,6 @@ def validate_config(config: dict[str, Any], dataset_name: str) -> None: sys.exit(1) -def validate_nowcast_config(nowcast_config: NowcastConfig) -> None: - """Validate that the options selected for the nowcast config are valid. - - Checks the nowcast config created from imported environment variables. - Exits immediately for invalid choices. - - Args: - nowcast_config: Instance of the NowcastConfig class - loaded from environment variables in config.py. - - Raises: - SystemExit: If any invalid choice is detected. - """ - if nowcast_config.ens_members != 1: - logger.error( - f"Invalid nowcast configuration: Currently, only ens_members=1 is supported. " - f"Current value: {nowcast_config.ens_members}. Exiting.\n" - ) - sys.exit(1) - - def validate_run_mode(run_mode: str, dataset_name: str) -> None: """Validate that the run mode is compatible with the dataset. From 09728fed43d77e670dc99b50cd8ef1a0efb80d3d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Wed, 1 Jul 2026 13:28:51 +0200 Subject: [PATCH 03/13] Configurable noise parameters --- sunflow/config.py | 18 ++++++++- sunflow/data_io.py | 4 -- sunflow/forecast.py | 90 ++++++++++++++++++++++++++++----------------- sunflow/main.py | 23 ++++++++++-- 4 files changed, 92 insertions(+), 43 deletions(-) diff --git a/sunflow/config.py b/sunflow/config.py index 6cb2e7f..bd13802 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -45,6 +45,8 @@ class NowcastConfig: nowcast_directory: str ens_members: int + alpha: float + beta: float past_steps: int future_steps: int input_data_availability_delay_minutes: int @@ -61,6 +63,8 @@ def from_env(cls) -> Self: - NOWCAST_DIRECTORY (default: .) - ENS_MEMBERS (default: 1) + - ALPHA (default: 0.0 for ENS_MEMBERS=1, 9.29 for ENS_MEMBERS>1) + - BETA (default: 0.0 for ENS_MEMBERS=1, 0.17 for ENS_MEMBERS>1) - PAST_STEPS (default: 4) - FUTURE_STEPS (default: 24) - INPUT_DATA_AVAILABILITY_DELAY_MINUTES (default: 24) @@ -69,9 +73,21 @@ def from_env(cls) -> Self: - SATELLITE_DATA_DIRECTORY (default: .) - MAX_CLEARSKY_FALLBACK_DAYS (default: 3) """ + + ens_members = int(os.getenv("ENS_MEMBERS", "1")) + # Reference for default noise values: + # A. Carpentieri, D. Folini, D. Nerini, S. Pulkkinen, M. Wild, A. Meyer, + # "Intraday probabilistic forecasts of surface solar radiation with cloud scale-dependent autoregressive advection," + # Applied Energy, Volume 351, 2023 + default_alpha = 0.0 if ens_members == 1 else 9.29 + default_beta = 0.0 if ens_members == 1 else 0.17 + + return cls( nowcast_directory=os.getenv("NOWCAST_DIRECTORY", "."), - ens_members=int(os.getenv("ENS_MEMBERS", "1")), + ens_members=ens_members, + alpha=float(os.getenv("ALPHA", str(default_alpha))), + beta=float(os.getenv("BETA", str(default_beta))), past_steps=int(os.getenv("PAST_STEPS", "4")), future_steps=int(os.getenv("FUTURE_STEPS", "24")), input_data_availability_delay_minutes=int( diff --git a/sunflow/data_io.py b/sunflow/data_io.py index 70973d3..2fcebd0 100644 --- a/sunflow/data_io.py +++ b/sunflow/data_io.py @@ -473,10 +473,6 @@ def save_forecast( ens_members = nowcast_config.ens_members filename = f"SolarNowcast_{time_step.strftime('%Y%m%d%H%M')}.nc" - # Add ensemble dimension if needed (forecast should be [ensemble, time, lat, lon]) - if forecast.ndim == 3: - forecast = forecast[np.newaxis, :, :, :] # Now [1, time, lat, lon] - # Build time coordinate (CF-convention: minutes since forecast reference time) time_step_naive = time_step.replace(tzinfo=None) diff --git a/sunflow/forecast.py b/sunflow/forecast.py index d55a7ce..5d151d8 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -4,7 +4,6 @@ import numpy as np import xarray as xr from Models.ProbabilisticAdvection import ProbabilisticAdvection -from loguru import logger from .geospatial import get_coordinates @@ -62,16 +61,19 @@ def preprocess_data( ) -def simple_advection_forecast( - ratio_data: np.ndarray, motion_field: np.ndarray, n_steps: int, ens_members: int +def probabilistic_advection_forecast( + ratio_data: np.ndarray, + motion_field: np.ndarray, + n_steps: int, + ens_members: int, + alpha: float, + beta: float, ) -> np.ndarray: - """Run a deterministic advection forecast on solar irradiance ratios. + """Run a probabilistic advection forecast on solar irradiance ratios. - Uses ProbabilisticAdvection with noise parameters alpha=0 and beta=0, - which disables Gaussian noise on the motion field norm and - von Mises noise on the direction, yielding a purely deterministic - advection result. The ensemble dimension added by the model is removed - before returning. + Uses ProbabilisticAdvection with configurable noise parameters: + alpha controls Gaussian noise on motion field norm and beta controls + von Mises noise on motion field direction. Args: ratio_data: Input array of shape (time, lat, lon) containing @@ -79,15 +81,18 @@ def simple_advection_forecast( motion_field: Optical flow field of shape (2, lat, lon) as produced by dense_lucaskanade. n_steps: Number of forecast timesteps to produce. + ens_members: Number of ensemble members. + alpha: Gaussian noise strength on motion field norm. + beta: von Mises noise strength on motion field angle. Returns: Forecast array of shape (n_steps, lat, lon). """ - # Initialize ProbabilisticAdvection with NO noise (alpha=0, beta=0) + # Initialize ProbabilisticAdvection with configured noise settings. pa = ProbabilisticAdvection( - alpha=0.0, # No Gaussian noise on motion field norm - beta=0.0, # No von Mises noise on motion field angle + alpha=alpha, + beta=beta, return_motion_field=False, ens_members=ens_members, ) @@ -149,7 +154,8 @@ def multiply_clearsky( "ratio_forecast time dimension does not match clearsky timesteps " f"({ratio_forecast.shape[0]} != {clearsky_stack.shape[0]})." ) - return ratio_forecast * clearsky_stack + solar_forecast = ratio_forecast * clearsky_stack + return solar_forecast[np.newaxis, :, :, :] # Add ensemble dimension for consistency if ratio_forecast.ndim == 4: if ratio_forecast.shape[1] != clearsky_stack.shape[0]: @@ -165,31 +171,47 @@ def multiply_clearsky( ) -def prepend_t0(clearsky_data: xr.Dataset, ratio_data: np.ndarray, solar_forecast: np.ndarray, config: dict, clearsky_t0_time: datetime) -> np.ndarray: +def prepend_t0( + clearsky_data: xr.Dataset, + ratio_data: np.ndarray, + solar_forecast: np.ndarray, + config: dict, + clearsky_t0_time: datetime, +) -> np.ndarray: + """Prepend analysis timestep (t=0) to an ensemble solar forecast. + + Computes the t=0 solar field as the latest observed ratio + (ratio_data[-1]) multiplied by clearsky irradiance at clearsky_t0_time, + then prepends that field to all ensemble members in solar_forecast. + + Args: + clearsky_data: Dataset containing clearsky irradiance values on + a time axis. + ratio_data: Ratio history array with shape (time, lat, lon). + solar_forecast: Forecast array with shape + (ensemble, forecast_time, lat, lon). + config: Runtime configuration dict containing + config["nc_variable_names"]["sds_cs"]. + clearsky_t0_time: Timestamp for the analysis clearsky field, + typically one day before the nowcast time. + + Returns: + Array of shape (ensemble, forecast_time + 1, lat, lon) + with the analysis field inserted at index 0 along the time axis. + """ # Prepend timestep 0: current observation (ratio_data[-1]) × clearsky at t=0 sds_cs_t0 = clearsky_data.sel(time=clearsky_t0_time.replace(tzinfo=None))[ config["nc_variable_names"]["sds_cs"] ].values solar_t0 = ratio_data[-1] * sds_cs_t0 - if solar_forecast.ndim == 3: - solar_forecast = np.concatenate( - [solar_t0[np.newaxis, :, :], solar_forecast], - axis=0, - ) - elif solar_forecast.ndim == 4: - # Broadcast the same t=0 clearsky-based analysis field to all ensembles. - solar_t0_ens = np.broadcast_to( - solar_t0, - (solar_forecast.shape[0],) + solar_t0.shape, - ) - solar_forecast = np.concatenate( - [solar_t0_ens[:, np.newaxis, :, :], solar_forecast], - axis=1, - ) - else: - raise ValueError( - "solar_forecast must have shape (time, lat, lon) or " - "(ensemble, time, lat, lon)." - ) + # Broadcast the same t=0 clearsky-based analysis field to all ensembles. + solar_t0_ens = np.broadcast_to( + solar_t0, + (solar_forecast.shape[0],) + solar_t0.shape, + ) + solar_forecast = np.concatenate( + [solar_t0_ens[:, np.newaxis, :, :], solar_forecast], + axis=1, + ) return solar_forecast diff --git a/sunflow/main.py b/sunflow/main.py index 40356ca..cbf87ae 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -22,7 +22,7 @@ save_forecast, ) from .downloaders import download_past_data -from .forecast import multiply_clearsky, preprocess_data, simple_advection_forecast, prepend_t0 +from .forecast import multiply_clearsky, preprocess_data, probabilistic_advection_forecast, prepend_t0 from .geospatial import check_solar_elevation, get_bbox from .time_handler import generate_time_steps, round_time from .validation import ( @@ -31,7 +31,6 @@ validate_clearsky_shapes, validate_config, validate_data_shape, - validate_nowcast_config, validate_run_mode, verify_environment_variables, ) @@ -274,12 +273,14 @@ def run_nowcast( # Compute motion field motion_field = dense_lucaskanade(ratio_data) - # Simple forecast (ratio forecast) - ratio_forecast = simple_advection_forecast( + # Probabilistic advection forecast (ratio forecast) + ratio_forecast = probabilistic_advection_forecast( ratio_data, motion_field, nowcast_config.future_steps, ens_members=nowcast_config.ens_members, + alpha=nowcast_config.alpha, + beta=nowcast_config.beta, ) # Generate previous day time steps for clearsky lookup @@ -405,6 +406,20 @@ def cli() -> None: logger.info(f"Running in {run_mode} mode") logger.info(f"Using {dataset_name} dataset") logger.info(f"Using {bbox_choice} bbox: {bbox}") + logger.info( + f"Using probabilistic advection noise parameters alpha={nowcast_config.alpha}, " + f"beta={nowcast_config.beta}" + ) + + if nowcast_config.ens_members == 1 and ( + nowcast_config.alpha != 0.0 or nowcast_config.beta != 0.0 + ): + logger.warning( + "Running with a single ensemble member, but non-zero probabilistic advection noise " \ + "parameters alpha and/or beta. This is generally not recommended as it " \ + "simply adds noise to the nowcast without providing any ensemble spread. " \ + "Consider setting alpha=0.0 and beta=0.0 for a single-member run." + ) validate_run_mode(run_mode, dataset_name) validate_config(config, dataset_name) From 64e17e1b30d6cff6340b50232851bf7601fb2515 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Wed, 1 Jul 2026 13:39:49 +0200 Subject: [PATCH 04/13] Linting --- sunflow/config.py | 6 +++--- sunflow/forecast.py | 4 +++- sunflow/main.py | 18 ++++++++++++------ sunflow/validation.py | 2 -- 4 files changed, 18 insertions(+), 12 deletions(-) diff --git a/sunflow/config.py b/sunflow/config.py index bd13802..bfa8a89 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -75,14 +75,14 @@ def from_env(cls) -> Self: """ ens_members = int(os.getenv("ENS_MEMBERS", "1")) - # Reference for default noise values: + # Reference for default noise values: # A. Carpentieri, D. Folini, D. Nerini, S. Pulkkinen, M. Wild, A. Meyer, - # "Intraday probabilistic forecasts of surface solar radiation with cloud scale-dependent autoregressive advection," + # "Intraday probabilistic forecasts of surface solar radiation with cloud + # scale-dependent autoregressive advection," # Applied Energy, Volume 351, 2023 default_alpha = 0.0 if ens_members == 1 else 9.29 default_beta = 0.0 if ens_members == 1 else 0.17 - return cls( nowcast_directory=os.getenv("NOWCAST_DIRECTORY", "."), ens_members=ens_members, diff --git a/sunflow/forecast.py b/sunflow/forecast.py index 5d151d8..439b394 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -155,7 +155,9 @@ def multiply_clearsky( f"({ratio_forecast.shape[0]} != {clearsky_stack.shape[0]})." ) solar_forecast = ratio_forecast * clearsky_stack - return solar_forecast[np.newaxis, :, :, :] # Add ensemble dimension for consistency + return solar_forecast[ + np.newaxis, :, :, : + ] # Add ensemble dimension for consistency if ratio_forecast.ndim == 4: if ratio_forecast.shape[1] != clearsky_stack.shape[0]: diff --git a/sunflow/main.py b/sunflow/main.py index cbf87ae..a1fbf4d 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -7,7 +7,6 @@ from enum import Enum import isodate -import numpy as np import yaml from loguru import logger from pysteps.motion.lucaskanade import dense_lucaskanade @@ -22,7 +21,12 @@ save_forecast, ) from .downloaders import download_past_data -from .forecast import multiply_clearsky, preprocess_data, probabilistic_advection_forecast, prepend_t0 +from .forecast import ( + multiply_clearsky, + prepend_t0, + preprocess_data, + probabilistic_advection_forecast, +) from .geospatial import check_solar_elevation, get_bbox from .time_handler import generate_time_steps, round_time from .validation import ( @@ -332,7 +336,9 @@ def run_nowcast( config["nc_variable_names"], ) - solar_forecast = prepend_t0(clearsky_data, ratio_data, solar_forecast, config, clearsky_t0_time) + solar_forecast = prepend_t0( + clearsky_data, ratio_data, solar_forecast, config, clearsky_t0_time + ) # Save forecast (now contains actual solar irradiance, not ratios) filename = save_forecast( @@ -415,9 +421,9 @@ def cli() -> None: nowcast_config.alpha != 0.0 or nowcast_config.beta != 0.0 ): logger.warning( - "Running with a single ensemble member, but non-zero probabilistic advection noise " \ - "parameters alpha and/or beta. This is generally not recommended as it " \ - "simply adds noise to the nowcast without providing any ensemble spread. " \ + "Running with a single ensemble member, but non-zero probabilistic advection" + "noise parameters alpha and/or beta. This is generally not recommended as it" + "simply adds noise to the nowcast without providing any ensemble spread. " "Consider setting alpha=0.0 and beta=0.0 for a single-member run." ) diff --git a/sunflow/validation.py b/sunflow/validation.py index ff1680a..09c9e2c 100644 --- a/sunflow/validation.py +++ b/sunflow/validation.py @@ -9,8 +9,6 @@ import xarray as xr from loguru import logger -from .config import NowcastConfig - class MissingClearskyDataError(RuntimeError): pass From d6460fd165f9509cb8742ac0fa942207c38a8de1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Wed, 1 Jul 2026 15:14:34 +0200 Subject: [PATCH 05/13] Add support for ensemble median and have that be default --- sunflow/data_io.py | 18 +++++++++++++----- sunflow/main.py | 43 ++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 55 insertions(+), 6 deletions(-) diff --git a/sunflow/data_io.py b/sunflow/data_io.py index 2fcebd0..4c65aaa 100644 --- a/sunflow/data_io.py +++ b/sunflow/data_io.py @@ -443,6 +443,7 @@ def save_forecast( dataset_name: str, nowcast_config: NowcastConfig, model_version: str, + output_mode: str, run_mode: str = "files", s3_config: S3Config | None = None, ) -> str: @@ -453,16 +454,18 @@ def save_forecast( numeric values (float64, minutes since the forecast reference time). Args: - forecast: Forecast array, shape [time, lat, lon] or - [ensemble, time, lat, lon]. + forecast: Forecast array, shape [ensemble, time, lat, lon]. time_step: Forecast reference time (start of the forecast window). n_steps: Number of forecast time steps to write. latitudes: 1-D array of latitude values (degrees). longitudes: 1-D array of longitude values (degrees). dataset_name: Name of the source dataset (options: KNMI, DWD). - nowcast_config: NowcastConfig object supplying output directory, - ensemble size, and input data frequency. + nowcast_config: NowcastConfig object supplying output directory + and input data frequency. model_version: Model version string written as a global attribute. + output_mode: Output aggregation mode label written to global + NetCDF attrs (expected: 'deterministic', 'median', + or 'full_ensemble'). run_mode: One of 'files' (local) or 's3'. Defaults to 'files'. s3_config: S3Config object; required when run_mode is 's3'. @@ -470,9 +473,13 @@ def save_forecast( Filename (basename only) of the written NetCDF file. """ input_data_frequency_minutes = nowcast_config.input_data_frequency_minutes - ens_members = nowcast_config.ens_members filename = f"SolarNowcast_{time_step.strftime('%Y%m%d%H%M')}.nc" + if forecast.ndim != 4: + raise ValueError("forecast must have shape (ensemble, time, lat, lon).") + + ens_members = forecast.shape[0] + # Build time coordinate (CF-convention: minutes since forecast reference time) time_step_naive = time_step.replace(tzinfo=None) @@ -521,6 +528,7 @@ def save_forecast( f"Simple Probabilistic Advection solar forecast " f"using {dataset_name} data" ), + "output_mode": output_mode, "history": ( f"Created " f"{datetime.now(timezone.utc).strftime('%Y-%m-%d %H:%M:%S')} UTC" diff --git a/sunflow/main.py b/sunflow/main.py index a1fbf4d..1043ca0 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -7,6 +7,7 @@ from enum import Enum import isodate +import numpy as np import yaml from loguru import logger from pysteps.motion.lucaskanade import dense_lucaskanade @@ -136,6 +137,14 @@ def parse_datetime_with_timezone(datetime_str: str) -> datetime: help="End of time span in ISO8601 format (inclusive). Use with --start-time.", default=None, ) + parser.add_argument( + "--full_ensemble", + action="store_true", + help=( + "Save full ensemble output. By default, the pixel-wise median across " + "ensemble members is saved." + ), + ) args = parser.parse_args() @@ -175,6 +184,7 @@ def run_nowcast( bbox_choice: str, nowcast_config: NowcastConfig, s3_config: S3Config, + full_ensemble: bool = False, custom_time: bool = True, ) -> RunResult: """Run a single nowcast for the given (already-rounded) time step. @@ -188,6 +198,8 @@ def run_nowcast( bbox_choice: Bounding box identifier. nowcast_config: NowcastConfig object. s3_config: S3Config object. + full_ensemble: If True, save all ensemble members. If False, + save pixel-wise median over ensemble members. custom_time: If True, skip the retry wait loop on missing data. Returns: @@ -340,9 +352,29 @@ def run_nowcast( clearsky_data, ratio_data, solar_forecast, config, clearsky_t0_time ) + if full_ensemble: + output_forecast = solar_forecast + output_mode = "full_ensemble" + logger.info("Saving full ensemble forecast") + else: + if solar_forecast.shape[0] == 1: + output_forecast = solar_forecast + output_mode = "deterministic" + logger.info( + "Saving deterministic forecast " + "(single ensemble member, kept as singleton ensemble dimension)" + ) + else: + output_forecast = np.median(solar_forecast, axis=0, keepdims=True) + output_mode = "median" + logger.info( + "Saving pixel-wise median forecast across ensemble members " + "(with singleton ensemble dimension)" + ) + # Save forecast (now contains actual solar irradiance, not ratios) filename = save_forecast( - solar_forecast, + output_forecast, time_step, nowcast_config.future_steps + 1, # +1 for the t=0 analysis step latitudes, @@ -350,6 +382,7 @@ def run_nowcast( dataset_name, nowcast_config, model_version, + output_mode, run_mode, s3_config, ) @@ -412,6 +445,13 @@ def cli() -> None: logger.info(f"Running in {run_mode} mode") logger.info(f"Using {dataset_name} dataset") logger.info(f"Using {bbox_choice} bbox: {bbox}") + logger.info(f"Number of ensemble members: {nowcast_config.ens_members}") + if args.full_ensemble: + logger.info("Output mode: full ensemble") + elif nowcast_config.ens_members == 1: + logger.info("Output mode: deterministic") + else: + logger.info("Output mode: median") logger.info( f"Using probabilistic advection noise parameters alpha={nowcast_config.alpha}, " f"beta={nowcast_config.beta}" @@ -479,6 +519,7 @@ def cli() -> None: bbox_choice, nowcast_config, s3_config, + full_ensemble=args.full_ensemble, custom_time=custom_time, ) results.append(result) From 3ea58373169911f4164ca92af7e8bc39af02f71d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Wed, 1 Jul 2026 15:34:03 +0200 Subject: [PATCH 06/13] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a51720..3c3231a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Added support for running an ensemble. Default output is the ensemble median, but --full_ensemble can be specified [!17](https://github.com/dmidk/sunflow/pull/17), @KristianHMoller - Added a check for the number of ensemble members, as the code currently supports only one [!13](https://github.com/dmidk/sunflow/pull/13), @KristianHMoller - Subsetting to bounding box is now also done in the `s3` and `files` code paths [!11](https://github.com/dmidk/sunflow/pull/11), @JoachimKoenigslieb - Subsetting to bounding box now correctly handles both ascending and descending lat/lon coordinates [!11](https://github.com/dmidk/sunflow/pull/11), @JoachimKoenigslieb From 63513d5b2d161641da6cd2490d5c5777af87da6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Sat, 11 Jul 2026 21:08:58 +0200 Subject: [PATCH 07/13] Make ensemble_members a CLI argument --- README.md | 3 ++- sunflow/config.py | 5 ++--- sunflow/main.py | 10 ++++++++-- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index da9dc8f..f936084 100644 --- a/README.md +++ b/README.md @@ -101,7 +101,6 @@ podman run -it --rm --entrypoint="" sunflow bash |----------|---------|-------------| | `NOWCAST_DIRECTORY` | `.` | Directory for forecast output files | | `SATELLITE_DATA_DIRECTORY` | `.` | Directory for input satellite data archive | -| `ENS_MEMBERS` | `1` | Number of ensemble members | | `PAST_STEPS` | `4` | Number of past time steps for motion field | | `FUTURE_STEPS` | `24` | Number of forecast time steps | | `INPUT_DATA_AVAILABILITY_DELAY_MINUTES` | `24` | Data availability delay (minutes) | @@ -140,6 +139,8 @@ podman run -it --rm --entrypoint="" sunflow bash - `--start-time` - Start of a time range in ISO8601 format (use with `--end-time`) - `--end-time` - End of a time range in ISO8601 format, inclusive (use with `--start-time`) - `--run_mode` - Specify run mode: `download` (fetch from API), `files` (local files), or `s3` (object storage) +- `--ensemble_members` - Number of ensemble members (Default 1) +- `--full_ensemble` - Specify that the full ensemble is the desired output rather than ensemble statistics ## Data Sources diff --git a/sunflow/config.py b/sunflow/config.py index bfa8a89..e718d22 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -56,13 +56,12 @@ class NowcastConfig: max_clearsky_fallback_days: int @classmethod - def from_env(cls) -> Self: + def from_env(cls, ensemble_members: int = 1) -> Self: """Load nowcast configuration from environment variables with defaults. Reads the following environment variables: - NOWCAST_DIRECTORY (default: .) - - ENS_MEMBERS (default: 1) - ALPHA (default: 0.0 for ENS_MEMBERS=1, 9.29 for ENS_MEMBERS>1) - BETA (default: 0.0 for ENS_MEMBERS=1, 0.17 for ENS_MEMBERS>1) - PAST_STEPS (default: 4) @@ -74,7 +73,7 @@ def from_env(cls) -> Self: - MAX_CLEARSKY_FALLBACK_DAYS (default: 3) """ - ens_members = int(os.getenv("ENS_MEMBERS", "1")) + ens_members = ensemble_members # Reference for default noise values: # A. Carpentieri, D. Folini, D. Nerini, S. Pulkkinen, M. Wild, A. Meyer, # "Intraday probabilistic forecasts of surface solar radiation with cloud diff --git a/sunflow/main.py b/sunflow/main.py index 1043ca0..da60dd3 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -137,6 +137,12 @@ def parse_datetime_with_timezone(datetime_str: str) -> datetime: help="End of time span in ISO8601 format (inclusive). Use with --start-time.", default=None, ) + parser.add_argument( + "--ensemble-members", + type=int, + default=1, + help="Number of ensemble members (default: 1)", + ) parser.add_argument( "--full_ensemble", action="store_true", @@ -428,9 +434,9 @@ def cli() -> None: ) # Load configuration - nowcast_config = NowcastConfig.from_env() - s3_config = S3Config.from_env() args = parse_arguments() + nowcast_config = NowcastConfig.from_env(ensemble_members=args.ensemble_members) + s3_config = S3Config.from_env() run_mode = args.run_mode dataset_name = args.dataset From 7ae579e74754fc18c0660266a6ac5f3471abbb5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Sat, 11 Jul 2026 21:58:53 +0200 Subject: [PATCH 08/13] Adding ensemble statistics to default output with more than one ensemble member --- CHANGELOG.md | 4 ++- README.md | 5 +++ sunflow/config.py | 48 ++++++++++++++++++++++++++++ sunflow/data_io.py | 76 ++++++++++++++++++++++++++++++++++----------- sunflow/forecast.py | 24 ++++++++++++++ sunflow/main.py | 26 ++++++++++------ 6 files changed, 154 insertions(+), 29 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c3231a..08f115e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,13 +9,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- Added support for running an ensemble. Default output is the ensemble median, but --full_ensemble can be specified [!17](https://github.com/dmidk/sunflow/pull/17), @KristianHMoller +- Added support for running an ensemble. Default output is ensemble statistics, but --full_ensemble can be specified [!17](https://github.com/dmidk/sunflow/pull/17), @KristianHMoller - Added a check for the number of ensemble members, as the code currently supports only one [!13](https://github.com/dmidk/sunflow/pull/13), @KristianHMoller - Subsetting to bounding box is now also done in the `s3` and `files` code paths [!11](https://github.com/dmidk/sunflow/pull/11), @JoachimKoenigslieb - Subsetting to bounding box now correctly handles both ascending and descending lat/lon coordinates [!11](https://github.com/dmidk/sunflow/pull/11), @JoachimKoenigslieb ### Changed +- Change number of ensemble members from being read from environment variable to being a CLI argument. [!17](https://github.com/dmidk/sunflow/pull/17), @KristianHMoller +- Ensemble output statistics are now configurable via the `ENSEMBLE_STATISTICS` environment variable, and selected statistics are written as separate NetCDF variables. [!17](https://github.com/dmidk/sunflow/pull/17), @KristianHMoller - Pass down number of ensembles from configurations to `ProbabilisticAdvection`. This gives a roughly 3x speedup in the no ensembles case [!12](https://github.com/dmidk/sunflow/pull/12), @JoachimKoenigslieb ## [v1.1.0] diff --git a/README.md b/README.md index f936084..da1164b 100644 --- a/README.md +++ b/README.md @@ -107,6 +107,7 @@ podman run -it --rm --entrypoint="" sunflow bash | `INPUT_DATA_FREQUENCY_MINUTES` | `15` | Data frequency (minutes) | | `MAX_WAITING_TIME_MINUTES` | `27` | Maximum wait time for data (minutes) | | `MAX_CLEARSKY_FALLBACK_DAYS` | `3` | Days back to search for fallback clear-sky data | +| `ENSEMBLE_STATISTICS` | `median,mean,p10,p25,p75,p90` | Comma-separated list of statistics for ensemble output. Allowed: `median`, `mean`, `p10`, `p25`, `p75`, `p90` (aliases `10th_percentile`, `25th_percentile`, `75th_percentile`, `90th_percentile` are accepted). | #### Data Source Configuration @@ -142,6 +143,10 @@ podman run -it --rm --entrypoint="" sunflow bash - `--ensemble_members` - Number of ensemble members (Default 1) - `--full_ensemble` - Specify that the full ensemble is the desired output rather than ensemble statistics +For ensemble runs (`--ensemble_members > 1`), default output is the configured +ensemble statistics from `ENSEMBLE_STATISTICS`. Use `--full_ensemble` to output +all members instead. + ## Data Sources - **KNMI**: MSG-CPP products — requires a free API key from the [KNMI Data Platform](https://dataplatform.knmi.nl/dataset/access/msg-cpp-products-1-0) diff --git a/sunflow/config.py b/sunflow/config.py index e718d22..b8c7377 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -3,6 +3,26 @@ from dataclasses import dataclass from typing import Self +DEFAULT_ENSEMBLE_STATISTICS = "median,mean,p10,p25,p75,p90" +_ALLOWED_STATISTICS = { + "median", + "mean", + "p10", + "p25", + "p75", + "p90", + "10th_percentile", + "25th_percentile", + "75th_percentile", + "90th_percentile", +} +_STATISTIC_ALIASES = { + "10th_percentile": "p10", + "25th_percentile": "p25", + "75th_percentile": "p75", + "90th_percentile": "p90", +} + # Predefined Bounding Box (BBOX) options BBOX_OPTIONS: dict[str, str | None] = { "DENMARK": "4,50,18,62", @@ -54,6 +74,7 @@ class NowcastConfig: max_waiting_time_minutes: int satellite_data_directory: str max_clearsky_fallback_days: int + ensemble_statistics: list[str] @classmethod def from_env(cls, ensemble_members: int = 1) -> Self: @@ -71,6 +92,7 @@ def from_env(cls, ensemble_members: int = 1) -> Self: - MAX_WAITING_TIME_MINUTES (default: 27) - SATELLITE_DATA_DIRECTORY (default: .) - MAX_CLEARSKY_FALLBACK_DAYS (default: 3) + - ENSEMBLE_STATISTICS (default: median,mean,p10,p25,p75,p90) """ ens_members = ensemble_members @@ -81,6 +103,9 @@ def from_env(cls, ensemble_members: int = 1) -> Self: # Applied Energy, Volume 351, 2023 default_alpha = 0.0 if ens_members == 1 else 9.29 default_beta = 0.0 if ens_members == 1 else 0.17 + statistics = _parse_ensemble_statistics( + os.getenv("ENSEMBLE_STATISTICS", DEFAULT_ENSEMBLE_STATISTICS) + ) return cls( nowcast_directory=os.getenv("NOWCAST_DIRECTORY", "."), @@ -98,4 +123,27 @@ def from_env(cls, ensemble_members: int = 1) -> Self: max_waiting_time_minutes=int(os.getenv("MAX_WAITING_TIME_MINUTES", "27")), satellite_data_directory=os.getenv("SATELLITE_DATA_DIRECTORY", "."), max_clearsky_fallback_days=int(os.getenv("MAX_CLEARSKY_FALLBACK_DAYS", "3")), + ensemble_statistics=statistics, ) + + +def _parse_ensemble_statistics(raw_statistics: str) -> list[str]: + """Parse and validate requested ensemble statistics from environment.""" + statistics = [ + token.strip().lower() + for token in raw_statistics.split(",") + if token.strip() + ] + if not statistics: + raise ValueError("ENSEMBLE_STATISTICS must contain at least one statistic") + + invalid = [stat for stat in statistics if stat not in _ALLOWED_STATISTICS] + if invalid: + raise ValueError( + "Invalid ENSEMBLE_STATISTICS value(s): " + f"{', '.join(invalid)}. Allowed values are: " + "median, mean, p10, p25, p75, p90, " + "10th_percentile, 25th_percentile, 75th_percentile, 90th_percentile" + ) + + return [_STATISTIC_ALIASES.get(stat, stat) for stat in statistics] diff --git a/sunflow/data_io.py b/sunflow/data_io.py index 4c65aaa..b67d537 100644 --- a/sunflow/data_io.py +++ b/sunflow/data_io.py @@ -435,7 +435,7 @@ def fetch_clearsky_with_fallback( def save_forecast( - forecast: np.ndarray, + forecast: np.ndarray | dict[str, np.ndarray], time_step: datetime, n_steps: int, latitudes: np.ndarray, @@ -447,14 +447,15 @@ def save_forecast( run_mode: str = "files", s3_config: S3Config | None = None, ) -> str: - """Save forecast array to a CF-compliant NetCDF4 file. + """Save forecast data to a CF-compliant NetCDF4 file. Writes the probabilistic advection forecast to either a local file or S3, depending on `run_mode`. The time coordinate is stored as CF-convention numeric values (float64, minutes since the forecast reference time). Args: - forecast: Forecast array, shape [ensemble, time, lat, lon]. + forecast: Forecast array with shape [ensemble, time, lat, lon], + or mapping of statistic name to arrays with the same shape. time_step: Forecast reference time (start of the forecast window). n_steps: Number of forecast time steps to write. latitudes: 1-D array of latitude values (degrees). @@ -464,7 +465,7 @@ def save_forecast( and input data frequency. model_version: Model version string written as a global attribute. output_mode: Output aggregation mode label written to global - NetCDF attrs (expected: 'deterministic', 'median', + NetCDF attrs (expected: 'deterministic', 'ensemble_statistics', or 'full_ensemble'). run_mode: One of 'files' (local) or 's3'. Defaults to 'files'. s3_config: S3Config object; required when run_mode is 's3'. @@ -475,23 +476,49 @@ def save_forecast( input_data_frequency_minutes = nowcast_config.input_data_frequency_minutes filename = f"SolarNowcast_{time_step.strftime('%Y%m%d%H%M')}.nc" - if forecast.ndim != 4: - raise ValueError("forecast must have shape (ensemble, time, lat, lon).") + data_vars: dict[str, tuple[list[str], np.ndarray, dict[str, str]]] + statistics_attr = "" - ens_members = forecast.shape[0] + if isinstance(forecast, dict): + if not forecast: + raise ValueError("forecast statistics mapping cannot be empty") - # Build time coordinate (CF-convention: minutes since forecast reference time) + data_vars = {} + first_shape: tuple[int, ...] | None = None + for statistic, values in forecast.items(): + if values.ndim != 4: + raise ValueError( + "Each statistic array must have shape " + "(ensemble, time, lat, lon)." + ) - time_step_naive = time_step.replace(tzinfo=None) - _time_units = f"minutes since {time_step_naive.strftime('%Y-%m-%d %H:%M:%S')}" - time_datetimes = [ - time_step_naive + timedelta(minutes=input_data_frequency_minutes * i) - for i in range(0, n_steps) - ] + if first_shape is None: + first_shape = values.shape + elif values.shape != first_shape: + raise ValueError("All statistic arrays must share the same shape") - ds = xr.Dataset( - { - "probabilistic_advection": ( + variable_name = f"GHI_probabilistic_advection_{statistic}" + data_vars[variable_name] = ( + ["ensemble", "time", "lat", "lon"], + values, + { + "description": ( + f"Probabilistic advection solar forecast ({statistic})" + ), + "long_name": "Surface downwelling solar radiation", + "units": "W m-2", + }, + ) + + ens_members = first_shape[0] if first_shape is not None else 0 + statistics_attr = ",".join(forecast.keys()) + else: + if forecast.ndim != 4: + raise ValueError("forecast must have shape (ensemble, time, lat, lon).") + + ens_members = forecast.shape[0] + data_vars = { + "GHI_probabilistic_advection": ( ["ensemble", "time", "lat", "lon"], forecast, { @@ -500,7 +527,19 @@ def save_forecast( "units": "W m-2", }, ), - }, + } + + # Build time coordinate (CF-convention: minutes since forecast reference time) + + time_step_naive = time_step.replace(tzinfo=None) + _time_units = f"minutes since {time_step_naive.strftime('%Y-%m-%d %H:%M:%S')}" + time_datetimes = [ + time_step_naive + timedelta(minutes=input_data_frequency_minutes * i) + for i in range(0, n_steps) + ] + + ds = xr.Dataset( + data_vars, coords={ "time": ( ["time"], @@ -534,6 +573,7 @@ def save_forecast( f"{datetime.now(timezone.utc).strftime('%Y-%m-%d %H:%M:%S')} UTC" ), "model_version": model_version, + **({"statistics": statistics_attr} if statistics_attr else {}), }, ) diff --git a/sunflow/forecast.py b/sunflow/forecast.py index 439b394..1fb6d93 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -217,3 +217,27 @@ def prepend_t0( axis=1, ) return solar_forecast + +def compute_ensemble_statistics( + forecast: np.ndarray, statistics: list[str] +) -> dict[str, np.ndarray]: + """Compute requested statistics over ensemble members (axis 0).""" + computed: dict[str, np.ndarray] = {} + for statistic in statistics: + match statistic: + case "median": + computed["median"] = np.median(forecast, axis=0, keepdims=True) + case "mean": + computed["mean"] = np.mean(forecast, axis=0, keepdims=True) + case "p10": + computed["p10"] = np.percentile(forecast, 10, axis=0, keepdims=True) + case "p25": + computed["p25"] = np.percentile(forecast, 25, axis=0, keepdims=True) + case "p75": + computed["p75"] = np.percentile(forecast, 75, axis=0, keepdims=True) + case "p90": + computed["p90"] = np.percentile(forecast, 90, axis=0, keepdims=True) + case _: # Defensive check; config parsing validates these values. + raise ValueError(f"Unsupported ensemble statistic: {statistic}") + + return computed diff --git a/sunflow/main.py b/sunflow/main.py index da60dd3..4e25117 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -27,6 +27,7 @@ prepend_t0, preprocess_data, probabilistic_advection_forecast, + compute_ensemble_statistics, ) from .geospatial import check_solar_elevation, get_bbox from .time_handler import generate_time_steps, round_time @@ -57,7 +58,6 @@ class RunResult: # Model version model_version = __version__ - def parse_arguments() -> argparse.Namespace: """Parse and validate command line arguments. @@ -138,7 +138,7 @@ def parse_datetime_with_timezone(datetime_str: str) -> datetime: default=None, ) parser.add_argument( - "--ensemble-members", + "--ensemble_members", type=int, default=1, help="Number of ensemble members (default: 1)", @@ -147,8 +147,8 @@ def parse_datetime_with_timezone(datetime_str: str) -> datetime: "--full_ensemble", action="store_true", help=( - "Save full ensemble output. By default, the pixel-wise median across " - "ensemble members is saved." + "Save full ensemble output. By default, ensemble statistics are " + "saved for ensemble runs." ), ) @@ -205,7 +205,7 @@ def run_nowcast( nowcast_config: NowcastConfig object. s3_config: S3Config object. full_ensemble: If True, save all ensemble members. If False, - save pixel-wise median over ensemble members. + save configured ensemble statistics over ensemble members. custom_time: If True, skip the retry wait loop on missing data. Returns: @@ -371,11 +371,14 @@ def run_nowcast( "(single ensemble member, kept as singleton ensemble dimension)" ) else: - output_forecast = np.median(solar_forecast, axis=0, keepdims=True) - output_mode = "median" + output_forecast = compute_ensemble_statistics( + solar_forecast, nowcast_config.ensemble_statistics + ) + output_mode = "ensemble_statistics" logger.info( - "Saving pixel-wise median forecast across ensemble members " - "(with singleton ensemble dimension)" + "Saving ensemble statistics across members: " + f"{', '.join(nowcast_config.ensemble_statistics)} " + "(each with singleton ensemble dimension)" ) # Save forecast (now contains actual solar irradiance, not ratios) @@ -457,7 +460,10 @@ def cli() -> None: elif nowcast_config.ens_members == 1: logger.info("Output mode: deterministic") else: - logger.info("Output mode: median") + logger.info( + "Output mode: ensemble statistics " + f"({', '.join(nowcast_config.ensemble_statistics)})" + ) logger.info( f"Using probabilistic advection noise parameters alpha={nowcast_config.alpha}, " f"beta={nowcast_config.beta}" From f406c0d7024429d583f9ccbd198ae0f4c841225e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Sat, 11 Jul 2026 22:00:06 +0200 Subject: [PATCH 09/13] Linting --- sunflow/config.py | 4 +--- sunflow/data_io.py | 3 +-- sunflow/forecast.py | 1 + sunflow/main.py | 4 ++-- 4 files changed, 5 insertions(+), 7 deletions(-) diff --git a/sunflow/config.py b/sunflow/config.py index b8c7377..c1f82be 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -130,9 +130,7 @@ def from_env(cls, ensemble_members: int = 1) -> Self: def _parse_ensemble_statistics(raw_statistics: str) -> list[str]: """Parse and validate requested ensemble statistics from environment.""" statistics = [ - token.strip().lower() - for token in raw_statistics.split(",") - if token.strip() + token.strip().lower() for token in raw_statistics.split(",") if token.strip() ] if not statistics: raise ValueError("ENSEMBLE_STATISTICS must contain at least one statistic") diff --git a/sunflow/data_io.py b/sunflow/data_io.py index b67d537..7fbeebd 100644 --- a/sunflow/data_io.py +++ b/sunflow/data_io.py @@ -488,8 +488,7 @@ def save_forecast( for statistic, values in forecast.items(): if values.ndim != 4: raise ValueError( - "Each statistic array must have shape " - "(ensemble, time, lat, lon)." + "Each statistic array must have shape " "(ensemble, time, lat, lon)." ) if first_shape is None: diff --git a/sunflow/forecast.py b/sunflow/forecast.py index 1fb6d93..139e306 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -218,6 +218,7 @@ def prepend_t0( ) return solar_forecast + def compute_ensemble_statistics( forecast: np.ndarray, statistics: list[str] ) -> dict[str, np.ndarray]: diff --git a/sunflow/main.py b/sunflow/main.py index 4e25117..7a30633 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -7,7 +7,6 @@ from enum import Enum import isodate -import numpy as np import yaml from loguru import logger from pysteps.motion.lucaskanade import dense_lucaskanade @@ -23,11 +22,11 @@ ) from .downloaders import download_past_data from .forecast import ( + compute_ensemble_statistics, multiply_clearsky, prepend_t0, preprocess_data, probabilistic_advection_forecast, - compute_ensemble_statistics, ) from .geospatial import check_solar_elevation, get_bbox from .time_handler import generate_time_steps, round_time @@ -58,6 +57,7 @@ class RunResult: # Model version model_version = __version__ + def parse_arguments() -> argparse.Namespace: """Parse and validate command line arguments. From 30b2b1ac5c7beb68be95bb3adadfd64fab52dc13 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Mon, 13 Jul 2026 21:15:34 +0200 Subject: [PATCH 10/13] Crop also in ensemble mode --- sunflow/forecast.py | 52 +++++++++++++++++++++++++++++++++++-------- sunflow/geospatial.py | 13 ++++++----- sunflow/main.py | 30 +++++++++++++++---------- 3 files changed, 68 insertions(+), 27 deletions(-) diff --git a/sunflow/forecast.py b/sunflow/forecast.py index 139e306..f3f8bc4 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -5,7 +5,7 @@ import xarray as xr from Models.ProbabilisticAdvection import ProbabilisticAdvection -from .geospatial import get_coordinates +from .geospatial import crop_forecast_to_domain, get_coordinates def preprocess_data( @@ -220,25 +220,59 @@ def prepend_t0( def compute_ensemble_statistics( - forecast: np.ndarray, statistics: list[str] + forecast: np.ndarray, + statistics: list[str], + latitudes: np.ndarray, + longitudes: np.ndarray, + domain_nowcast: str, ) -> dict[str, np.ndarray]: """Compute requested statistics over ensemble members (axis 0).""" computed: dict[str, np.ndarray] = {} for statistic in statistics: match statistic: case "median": - computed["median"] = np.median(forecast, axis=0, keepdims=True) + computed["median"], latitudes, longitudes = crop_forecast_to_domain( + np.median(forecast, axis=0, keepdims=True), + latitudes, + longitudes, + domain_nowcast, + ) case "mean": - computed["mean"] = np.mean(forecast, axis=0, keepdims=True) + computed["mean"], latitudes, longitudes = crop_forecast_to_domain( + np.mean(forecast, axis=0, keepdims=True), + latitudes, + longitudes, + domain_nowcast, + ) case "p10": - computed["p10"] = np.percentile(forecast, 10, axis=0, keepdims=True) + computed["p10"], latitudes, longitudes = crop_forecast_to_domain( + np.percentile(forecast, 10, axis=0, keepdims=True), + latitudes, + longitudes, + domain_nowcast, + ) case "p25": - computed["p25"] = np.percentile(forecast, 25, axis=0, keepdims=True) + computed["p25"], latitudes, longitudes = crop_forecast_to_domain( + np.percentile(forecast, 25, axis=0, keepdims=True), + latitudes, + longitudes, + domain_nowcast, + ) case "p75": - computed["p75"] = np.percentile(forecast, 75, axis=0, keepdims=True) + computed["p75"], latitudes, longitudes = crop_forecast_to_domain( + np.percentile(forecast, 75, axis=0, keepdims=True), + latitudes, + longitudes, + domain_nowcast, + ) case "p90": - computed["p90"] = np.percentile(forecast, 90, axis=0, keepdims=True) + computed["p90"], latitudes, longitudes = crop_forecast_to_domain( + np.percentile(forecast, 90, axis=0, keepdims=True), + latitudes, + longitudes, + domain_nowcast, + ) case _: # Defensive check; config parsing validates these values. raise ValueError(f"Unsupported ensemble statistic: {statistic}") - return computed + return computed, latitudes, longitudes diff --git a/sunflow/geospatial.py b/sunflow/geospatial.py index 1fc1560..2b91f9a 100644 --- a/sunflow/geospatial.py +++ b/sunflow/geospatial.py @@ -234,10 +234,10 @@ def crop_forecast_to_domain( longitudes: np.ndarray, domain_bbox: str, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Crop a [time, lat, lon] forecast and its coordinates to a domain bbox. + """Crop a forecast and its coordinates to a domain bbox. Args: - forecast: Forecast array with shape [time, lat, lon]. + forecast: Forecast array with shape [ensemble, time, lat, lon]. latitudes: 1-D latitude array. longitudes: 1-D longitude array. domain_bbox: Requested bbox string lon_min,lat_min,lon_max,lat_max. @@ -246,12 +246,13 @@ def crop_forecast_to_domain( Tuple (cropped_forecast, cropped_latitudes, cropped_longitudes). Raises: - RuntimeError: If forecast dimensionality is not [time, lat, lon] or the + RuntimeError: If forecast dimensionality is not [ensemble, time, lat, lon] or the requested domain has no overlap with the provided coordinates. """ - if forecast.ndim != 3: + if forecast.ndim != 4: raise RuntimeError( - f"Expected forecast shape [time, lat, lon], got {forecast.shape}." + "Expected forecast shape [ensemble, time, lat, lon], " + f"got {forecast.shape}." ) lon_min, lat_min, lon_max, lat_max = parse_bbox(domain_bbox) @@ -278,7 +279,7 @@ def crop_forecast_to_domain( f"Requested domain_nowcast={domain_bbox} does not overlap forecast grid." ) - cropped_forecast = forecast[:, lat_idx, :][:, :, lon_idx] + cropped_forecast = forecast[:, :, lat_idx, :][:, :, :, lon_idx] return cropped_forecast, latitudes[lat_idx], longitudes[lon_idx] diff --git a/sunflow/main.py b/sunflow/main.py index 61d337b..d5539d6 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -28,7 +28,6 @@ preprocess_data, probabilistic_advection_forecast, ) - from .geospatial import ( check_solar_elevation, crop_forecast_to_domain, @@ -402,20 +401,34 @@ def run_nowcast( ) if full_ensemble: - output_forecast = solar_forecast + output_forecast, latitudes, longitudes = crop_forecast_to_domain( + solar_forecast, + latitudes, + longitudes, + domain_nowcast, + ) output_mode = "full_ensemble" logger.info("Saving full ensemble forecast") else: if solar_forecast.shape[0] == 1: - output_forecast = solar_forecast + output_forecast, latitudes, longitudes = crop_forecast_to_domain( + solar_forecast, + latitudes, + longitudes, + domain_nowcast, + ) output_mode = "deterministic" logger.info( "Saving deterministic forecast " "(single ensemble member, kept as singleton ensemble dimension)" ) else: - output_forecast = compute_ensemble_statistics( - solar_forecast, nowcast_config.ensemble_statistics + output_forecast, latitudes, longitudes = compute_ensemble_statistics( + solar_forecast, + nowcast_config.ensemble_statistics, + latitudes, + longitudes, + domain_nowcast, ) output_mode = "ensemble_statistics" logger.info( @@ -424,13 +437,6 @@ def run_nowcast( "(each with singleton ensemble dimension)" ) - solar_forecast, latitudes, longitudes = crop_forecast_to_domain( - solar_forecast, - latitudes, - longitudes, - domain_nowcast, - ) - # Save forecast (now contains actual solar irradiance, not ratios) filename = save_forecast( output_forecast, From 0349137153f9154edc5d9779c009203cbe883ad5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Mon, 13 Jul 2026 21:21:38 +0200 Subject: [PATCH 11/13] Change PA ensemble noise to match paper --- sunflow/config.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/sunflow/config.py b/sunflow/config.py index d4daeac..b8faa5c 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -89,8 +89,8 @@ def from_env(cls, ensemble_members: int = 1) -> Self: Reads the following environment variables: - NOWCAST_DIRECTORY (default: .) - - ALPHA (default: 0.0 for ENS_MEMBERS=1, 9.29 for ENS_MEMBERS>1) - - BETA (default: 0.0 for ENS_MEMBERS=1, 0.17 for ENS_MEMBERS>1) + - ALPHA (default: 0.0 for ENS_MEMBERS=1, 9.23 for ENS_MEMBERS>1) + - BETA (default: 0.0 for ENS_MEMBERS=1, 0.15 for ENS_MEMBERS>1) - PAST_STEPS (default: 4) - FUTURE_STEPS (default: 24) - INPUT_DATA_AVAILABILITY_DELAY_MINUTES (default: 24) @@ -107,8 +107,8 @@ def from_env(cls, ensemble_members: int = 1) -> Self: # "Intraday probabilistic forecasts of surface solar radiation with cloud # scale-dependent autoregressive advection," # Applied Energy, Volume 351, 2023 - default_alpha = 0.0 if ens_members == 1 else 9.29 - default_beta = 0.0 if ens_members == 1 else 0.17 + default_alpha = 0.0 if ens_members == 1 else 9.23 + default_beta = 0.0 if ens_members == 1 else 0.15 statistics = _parse_ensemble_statistics( os.getenv("ENSEMBLE_STATISTICS", DEFAULT_ENSEMBLE_STATISTICS) ) From ba0a5a685f0844b3db318b7fe6ee7e33f6793c97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Mon, 13 Jul 2026 21:24:59 +0200 Subject: [PATCH 12/13] Linting --- sunflow/validation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sunflow/validation.py b/sunflow/validation.py index a378932..7e7ad62 100644 --- a/sunflow/validation.py +++ b/sunflow/validation.py @@ -10,6 +10,9 @@ import xarray as xr from loguru import logger +from .geospatial import parse_bbox + + class MissingClearskyDataError(RuntimeError): pass From e2af071f757b34ea735ec967d244d0bb2a518735 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= <150110122+KristianHMoller@users.noreply.github.com> Date: Sun, 19 Jul 2026 22:29:01 +0200 Subject: [PATCH 13/13] Update to using the version of SolarSTEPS with arctan2 (#19) --- CHANGELOG.md | 1 + pyproject.toml | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 90b1952..90b7918 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Change number of ensemble members from being read from environment variable to being a CLI argument. [!17](https://github.com/dmidk/sunflow/pull/17), @KristianHMoller - Ensemble output statistics are now configurable via the `ENSEMBLE_STATISTICS` environment variable, and selected statistics are written as separate NetCDF variables. [!17](https://github.com/dmidk/sunflow/pull/17), @KristianHMoller +- Update SolarSTEPS dependency to use the version with arctan2 [!19](https://github.com/dmidk/sunflow/pull/19), @KristianHMoller - Modified subset_to_bbox function to slice based on edges rather than centers [!10](https://github.com/dmidk/sunflow/pull/10), @KristianHMoller - Pass down number of ensembles from configurations to `ProbabilisticAdvection`. This gives a roughly 3x speedup in the no ensembles case [!12](https://github.com/dmidk/sunflow/pull/12), @JoachimKoenigslieb diff --git a/pyproject.toml b/pyproject.toml index 57f5e34..fe97637 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,7 +18,7 @@ dependencies = [ "pvlib>=0.10.0", "pysteps>=1.7.0", "opencv-python-headless>=4.5.0", - "SolarSTEPS @ git+https://github.com/dmidk/SolarSTEPS@18012d7b11c56895e14abc39c3800f68a9a44ecb", + "SolarSTEPS @ git+https://github.com/dmidk/SolarSTEPS@3def2b36974bc2fb92f2ff14fc40dc8da2a9c547", "loguru>=0.7.3", "isodate>=0.7.2", "fsspec>=2023.1.0",