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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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 0abe53602cd3ae3e65733b918ee555d097f5b790 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Thu, 2 Jul 2026 11:45:00 +0200 Subject: [PATCH 07/10] Add support for solarSTEPS --- sunflow/config.py | 7 ++++ sunflow/forecast.py | 85 +++++++++++++++++++++++++++++++++++++++++++++ sunflow/main.py | 9 +++++ 3 files changed, 101 insertions(+) diff --git a/sunflow/config.py b/sunflow/config.py index bfa8a89..0b50430 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -10,6 +10,13 @@ "CUSTOM": None, } +# Centerpoint of lowest object resolved in number of grid points. +# Used for determining the optimal number of cascades for SolarSTEPS. +# Reverse engineered from the optimal values identified by Carpentieri et al. +# "Intraday probabilistic forecasts of surface solar radiation with +# cloud scale-dependent autoregressive advection" +TARGET_SMALLEST_RESOLUTION = 1.3515 + @dataclass class S3Config: diff --git a/sunflow/forecast.py b/sunflow/forecast.py index 439b394..cd1971f 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -3,8 +3,12 @@ import numpy as np import xarray as xr +from loguru import logger from Models.ProbabilisticAdvection import ProbabilisticAdvection +from Models.SolarSTEPS import SolarSTEPS +from pysteps.cascade.bandpass_filters import _gaussweights_1d +from .config import TARGET_SMALLEST_RESOLUTION from .geospatial import get_coordinates @@ -102,6 +106,87 @@ def probabilistic_advection_forecast( return forecast +def solarsteps_forecast( + ratio_data: np.ndarray, + motion_field: np.ndarray, + n_steps: int, + ens_members: int, +) -> np.ndarray: + """Run a SolarSTEPS ensemble forecast on solar irradiance ratios. + + Configures SolarSTEPS with fixed stochastic-noise and normalization + options, determines the number of cascade levels from the input grid + size, and generates an ensemble forecast for the requested lead times. + + Args: + ratio_data: Input array of shape (time, lat, lon) containing + SDS/SDS_CS ratios for the past timesteps. + 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. + + Returns: + Forecast array of shape (ensemble, n_steps, lat, lon). + """ + + n_cascade_levels = determine_cascade_levels(ratio_data) + solarsteps = SolarSTEPS( + ar_order=1, + n_cascade_levels=n_cascade_levels, + probmatching=True, + norm=True, + local=False, + noise_kwargs={ + "noise_win_size": 90, + "noise_std_win_size": 15, + "noise_method": "local-SSFT", + }, + norm_kwargs={"extra_normalization": True}, + verbose=False, + ) + forecast = solarsteps.ensemble_forecast( + ratio_data, motion_field, n_steps, seeds=np.arange(ens_members) + ) + + return forecast + + +def determine_cascade_levels(ratio_data: np.ndarray) -> int: + """Determine the number of cascade levels for SolarSTEPS based on grid size. + + Increases the number of cascade levels until the coarsest Gaussian + cascade resolution reaches or undershoots TARGET_SMALLEST_RESOLUTION + for the largest horizontal grid dimension. + + Args: + ratio_data: Input array of shape (time, lat, lon) containing + SDS/SDS_CS ratios for the past timesteps. + + Returns: + Number of cascade levels. + + Raises: + ValueError: If TARGET_SMALLEST_RESOLUTION is less than or equal to 1. + """ + + if TARGET_SMALLEST_RESOLUTION <= 1: + raise ValueError("TARGET_SMALLEST_RESOLUTION must be greater than 1.") + + largest_dimension = max(ratio_data.shape[1], ratio_data.shape[2]) + n_cascade_levels = 1 + while True: + _, resolution_centers = _gaussweights_1d( + largest_dimension, n_cascade_levels, gauss_scale=0.5 + ) + if resolution_centers[0] <= TARGET_SMALLEST_RESOLUTION: + logger.info( + f"Optimal number of cascades identified to be: {n_cascade_levels}" + ) + return n_cascade_levels + n_cascade_levels += 1 + + def multiply_clearsky( ratio_forecast: np.ndarray, clearsky_data: xr.Dataset, diff --git a/sunflow/main.py b/sunflow/main.py index 1043ca0..5a9201b 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -27,6 +27,7 @@ prepend_t0, preprocess_data, probabilistic_advection_forecast, + solarsteps_forecast, ) from .geospatial import check_solar_elevation, get_bbox from .time_handler import generate_time_steps, round_time @@ -299,6 +300,14 @@ def run_nowcast( beta=nowcast_config.beta, ) + # SolarSTEPS forecast (ratio forecast) + ratio_forecast = solarsteps_forecast( + ratio_data, + motion_field, + nowcast_config.future_steps, + ens_members=nowcast_config.ens_members, + ) + # Generate previous day time steps for clearsky lookup previous_day_time_steps = generate_time_steps( time_step, From c44503404e1b7a4eeda2189dee895c9c36508c52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Thu, 2 Jul 2026 12:53:59 +0200 Subject: [PATCH 08/10] Adding CLI argument for forecast model --- sunflow/data_io.py | 15 ++++++++--- sunflow/main.py | 63 ++++++++++++++++++++++++++++------------------ 2 files changed, 50 insertions(+), 28 deletions(-) diff --git a/sunflow/data_io.py b/sunflow/data_io.py index 4c65aaa..19b71cb 100644 --- a/sunflow/data_io.py +++ b/sunflow/data_io.py @@ -444,6 +444,7 @@ def save_forecast( nowcast_config: NowcastConfig, model_version: str, output_mode: str, + forecast_model: str, run_mode: str = "files", s3_config: S3Config | None = None, ) -> str: @@ -489,13 +490,20 @@ def save_forecast( for i in range(0, n_steps) ] + if forecast_model == "probabilistic_advection": + ghi_variable_name = "GHI_probabilistic_advection" + elif forecast_model == "solarsteps": + ghi_variable_name = "GHI_solarsteps" + else: + raise ValueError(f"Invalid forecast_model '{forecast_model}'. ") + ds = xr.Dataset( { - "probabilistic_advection": ( + ghi_variable_name: ( ["ensemble", "time", "lat", "lon"], forecast, { - "description": "Probabilistic advection solar forecast", + "description": f"{forecast_model} solar forecast", "long_name": "Surface downwelling solar radiation", "units": "W m-2", }, @@ -525,8 +533,7 @@ def save_forecast( }, attrs={ "description": ( - f"Simple Probabilistic Advection solar forecast " - f"using {dataset_name} data" + f"{forecast_model} solar forecast " f"using {dataset_name} data" ), "output_mode": output_mode, "history": ( diff --git a/sunflow/main.py b/sunflow/main.py index 5a9201b..1338324 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -141,10 +141,14 @@ def parse_datetime_with_timezone(datetime_str: str) -> datetime: parser.add_argument( "--full_ensemble", action="store_true", - help=( - "Save full ensemble output. By default, the pixel-wise median across " - "ensemble members is saved." - ), + help="Save full ensemble output. By default, the pixel-wise median across " + "ensemble members is saved.", + ) + parser.add_argument( + "--forecast_model", + choices=["probabilistic_advection", "solarsteps"], + default="probabilistic_advection", + help="Choose forecast model (default: probabilistic_advection)", ) args = parser.parse_args() @@ -179,6 +183,7 @@ def parse_datetime_with_timezone(datetime_str: str) -> datetime: def run_nowcast( time_step: datetime, run_mode: str, + forecast_model: str, config: dict, bbox: str, dataset_name: str, @@ -193,6 +198,7 @@ def run_nowcast( Args: time_step: The time step to produce a forecast for. run_mode: One of 'download', 'files', or 's3'. + forecast_model: Forecast model to use ('probabilistic_advection' or 'solarsteps'). config: Dataset configuration dict. bbox: Bounding box string. dataset_name: Name of dataset. @@ -290,23 +296,27 @@ def run_nowcast( # Compute motion field motion_field = dense_lucaskanade(ratio_data) - # 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, - ) + if forecast_model == "probabilistic_advection": - # SolarSTEPS forecast (ratio forecast) - ratio_forecast = solarsteps_forecast( - ratio_data, - motion_field, - nowcast_config.future_steps, - ens_members=nowcast_config.ens_members, - ) + # 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, + ) + elif forecast_model == "solarsteps": + # SolarSTEPS forecast (ratio forecast) + ratio_forecast = solarsteps_forecast( + ratio_data, + motion_field, + nowcast_config.future_steps, + ens_members=nowcast_config.ens_members, + ) + else: + raise ValueError(f"Unknown forecast model: {forecast_model}") # Generate previous day time steps for clearsky lookup previous_day_time_steps = generate_time_steps( @@ -392,6 +402,7 @@ def run_nowcast( nowcast_config, model_version, output_mode, + forecast_model, run_mode, s3_config, ) @@ -445,6 +456,7 @@ def cli() -> None: dataset_name = args.dataset bbox_choice = args.bbox bbox = get_bbox(bbox_choice, args.custom_bbox) + forecast_model = args.forecast_model config = yaml.safe_load(open("config.yaml"))[dataset_name] @@ -454,6 +466,7 @@ 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"Forecast model: {forecast_model}") logger.info(f"Number of ensemble members: {nowcast_config.ens_members}") if args.full_ensemble: logger.info("Output mode: full ensemble") @@ -461,10 +474,11 @@ def cli() -> None: 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}" - ) + if forecast_model == "probabilistic_advection": + logger.info( + "Using probabilistic advection noise parameters " + f"alpha={nowcast_config.alpha}, beta={nowcast_config.beta}" + ) if nowcast_config.ens_members == 1 and ( nowcast_config.alpha != 0.0 or nowcast_config.beta != 0.0 @@ -522,6 +536,7 @@ def cli() -> None: result = run_nowcast( time_step, run_mode, + forecast_model, config, bbox, dataset_name, From 132e3ee720edf43858088d5e15fb82d5f044627b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Thu, 2 Jul 2026 14:36:34 +0200 Subject: [PATCH 09/10] Make SolarSTEPS noise environment variables --- sunflow/config.py | 8 ++++++++ sunflow/forecast.py | 7 +++++-- sunflow/main.py | 8 ++++++++ 3 files changed, 21 insertions(+), 2 deletions(-) diff --git a/sunflow/config.py b/sunflow/config.py index 0b50430..0d0efa2 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -54,6 +54,8 @@ class NowcastConfig: ens_members: int alpha: float beta: float + noise_win_size: int + noise_std_win_size: int past_steps: int future_steps: int input_data_availability_delay_minutes: int @@ -89,12 +91,18 @@ def from_env(cls) -> 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 + default_noise_win_size = 0 if ens_members == 1 else 90 + default_noise_std_win_size = 0 if ens_members == 1 else 15 return cls( nowcast_directory=os.getenv("NOWCAST_DIRECTORY", "."), ens_members=ens_members, alpha=float(os.getenv("ALPHA", str(default_alpha))), beta=float(os.getenv("BETA", str(default_beta))), + noise_win_size=int(os.getenv("NOISE_WIN_SIZE", str(default_noise_win_size))), + noise_std_win_size=int( + os.getenv("NOISE_STD_WIN_SIZE", str(default_noise_std_win_size)) + ), 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/forecast.py b/sunflow/forecast.py index cd1971f..5e84f08 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -111,6 +111,8 @@ def solarsteps_forecast( motion_field: np.ndarray, n_steps: int, ens_members: int, + noise_win_size: float, + noise_std_win_size: float, ) -> np.ndarray: """Run a SolarSTEPS ensemble forecast on solar irradiance ratios. @@ -126,6 +128,7 @@ def solarsteps_forecast( n_steps: Number of forecast timesteps to produce. ens_members: Number of ensemble members. + Returns: Forecast array of shape (ensemble, n_steps, lat, lon). """ @@ -138,8 +141,8 @@ def solarsteps_forecast( norm=True, local=False, noise_kwargs={ - "noise_win_size": 90, - "noise_std_win_size": 15, + "noise_win_size": noise_win_size, + "noise_std_win_size": noise_std_win_size, "noise_method": "local-SSFT", }, norm_kwargs={"extra_normalization": True}, diff --git a/sunflow/main.py b/sunflow/main.py index 1338324..e5532a6 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -314,6 +314,8 @@ def run_nowcast( motion_field, nowcast_config.future_steps, ens_members=nowcast_config.ens_members, + noise_win_size=nowcast_config.noise_win_size, + noise_std_win_size=nowcast_config.noise_std_win_size, ) else: raise ValueError(f"Unknown forecast model: {forecast_model}") @@ -479,6 +481,12 @@ def cli() -> None: "Using probabilistic advection noise parameters " f"alpha={nowcast_config.alpha}, beta={nowcast_config.beta}" ) + elif forecast_model == "solarsteps": + logger.info( + "Using SolarSTEPS noise parameters " + f"noise_win_size={nowcast_config.noise_win_size}, " + f"noise_std_win_size={nowcast_config.noise_std_win_size}" + ) if nowcast_config.ens_members == 1 and ( nowcast_config.alpha != 0.0 or nowcast_config.beta != 0.0 From ee3785fa06ea5cfa0f13d348eb0fdc873fcc2d84 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20H=2E=20M=C3=B8ller?= Date: Fri, 3 Jul 2026 11:34:54 +0200 Subject: [PATCH 10/10] Improve descriptions for ensemble --- sunflow/config.py | 16 +++++++++------- sunflow/data_io.py | 11 +++++++++++ sunflow/forecast.py | 6 +++++- sunflow/main.py | 26 +++++++++++++++++--------- 4 files changed, 42 insertions(+), 17 deletions(-) diff --git a/sunflow/config.py b/sunflow/config.py index 0d0efa2..a36dcfc 100644 --- a/sunflow/config.py +++ b/sunflow/config.py @@ -54,6 +54,7 @@ class NowcastConfig: ens_members: int alpha: float beta: float + noise_method: str | None noise_win_size: int noise_std_win_size: int past_steps: int @@ -84,25 +85,26 @@ def from_env(cls) -> Self: """ ens_members = int(os.getenv("ENS_MEMBERS", "1")) - # Reference for default noise values: + # Probabilistic advection noise parameters + # (default values are based on the original implementation by Carpentieri et al.) + # Reference for default probabilistic advection 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 - default_noise_win_size = 0 if ens_members == 1 else 90 - default_noise_std_win_size = 0 if ens_members == 1 else 15 + + # Default SolarSTEPS noise parameters are based on the original implementation return cls( nowcast_directory=os.getenv("NOWCAST_DIRECTORY", "."), ens_members=ens_members, alpha=float(os.getenv("ALPHA", str(default_alpha))), beta=float(os.getenv("BETA", str(default_beta))), - noise_win_size=int(os.getenv("NOISE_WIN_SIZE", str(default_noise_win_size))), - noise_std_win_size=int( - os.getenv("NOISE_STD_WIN_SIZE", str(default_noise_std_win_size)) - ), + noise_method=os.getenv("NOISE_METHOD", "local-SSFT"), + noise_win_size=int(os.getenv("NOISE_WIN_SIZE", "90")), + noise_std_win_size=int(os.getenv("NOISE_STD_WIN_SIZE", "15")), 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 19b71cb..e6e334a 100644 --- a/sunflow/data_io.py +++ b/sunflow/data_io.py @@ -480,6 +480,14 @@ def save_forecast( raise ValueError("forecast must have shape (ensemble, time, lat, lon).") ens_members = forecast.shape[0] + if forecast_model == "probabilistic_advection": + noise_values = f"alpha: {nowcast_config.alpha}, beta: {nowcast_config.beta}" + elif forecast_model == "solarsteps": + noise_values = ( + f"noise_method: {nowcast_config.noise_method}, " + f"noise_win_size: {nowcast_config.noise_win_size}, " + f"noise_std_win_size: {nowcast_config.noise_std_win_size}" + ) # Build time coordinate (CF-convention: minutes since forecast reference time) @@ -541,6 +549,9 @@ def save_forecast( f"{datetime.now(timezone.utc).strftime('%Y-%m-%d %H:%M:%S')} UTC" ), "model_version": model_version, + "forecast_model": forecast_model, + "input_data_source": dataset_name, + "noise_values": noise_values, }, ) diff --git a/sunflow/forecast.py b/sunflow/forecast.py index 5e84f08..09480a3 100644 --- a/sunflow/forecast.py +++ b/sunflow/forecast.py @@ -111,6 +111,7 @@ def solarsteps_forecast( motion_field: np.ndarray, n_steps: int, ens_members: int, + noise_method: str | None, noise_win_size: float, noise_std_win_size: float, ) -> np.ndarray: @@ -127,6 +128,9 @@ def solarsteps_forecast( produced by dense_lucaskanade. n_steps: Number of forecast timesteps to produce. ens_members: Number of ensemble members. + noise_method: Noise method to use (e.g., "local-SSFT"). + noise_win_size: Window size for noise calculation. + noise_std_win_size: Standard deviation window size for noise calculation. Returns: @@ -143,7 +147,7 @@ def solarsteps_forecast( noise_kwargs={ "noise_win_size": noise_win_size, "noise_std_win_size": noise_std_win_size, - "noise_method": "local-SSFT", + "noise_method": noise_method, }, norm_kwargs={"extra_normalization": True}, verbose=False, diff --git a/sunflow/main.py b/sunflow/main.py index e5532a6..ba0ae47 100644 --- a/sunflow/main.py +++ b/sunflow/main.py @@ -314,6 +314,7 @@ def run_nowcast( motion_field, nowcast_config.future_steps, ens_members=nowcast_config.ens_members, + noise_method=nowcast_config.noise_method, noise_win_size=nowcast_config.noise_win_size, noise_std_win_size=nowcast_config.noise_std_win_size, ) @@ -484,19 +485,26 @@ def cli() -> None: elif forecast_model == "solarsteps": logger.info( "Using SolarSTEPS noise parameters " + f"noise_method={nowcast_config.noise_method}, " f"noise_win_size={nowcast_config.noise_win_size}, " f"noise_std_win_size={nowcast_config.noise_std_win_size}" ) - 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." - ) + if nowcast_config.ens_members == 1: + if 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." + ) + elif forecast_model == "solarsteps": + logger.warning( + "Running SolarSTEPS with a single ensemble member. " + "This is generally not recommended as it simply adds noise to the " + "nowcast without providing any ensemble spread." + ) validate_run_mode(run_mode, dataset_name) validate_config(config, dataset_name)