Add SolarSTEPS forecast functionality

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

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:
@@ -45,6 +52,11 @@ class NowcastConfig:
 
     nowcast_directory: str
     ens_members: int
+    alpha: float
+    beta: float
+    noise_method: str | None
+    noise_win_size: int
+    noise_std_win_size: int
     past_steps: int
     future_steps: int
     input_data_availability_delay_minutes: int
@@ -61,6 +73,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 +83,28 @@ def from_env(cls) -> Self:
         - SATELLITE_DATA_DIRECTORY (default: .)
         - MAX_CLEARSKY_FALLBACK_DAYS (default: 3)
         """
+
+        ens_members = int(os.getenv("ENS_MEMBERS", "1"))
+        # 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 SolarSTEPS noise parameters are based on the original implementation
+
         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))),
+            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(

sunflow/data_io.py

@@ -443,6 +443,8 @@ def save_forecast(
     dataset_name: str,
     nowcast_config: NowcastConfig,
     model_version: str,
+    output_mode: str,
+    forecast_model: str,
     run_mode: str = "files",
     s3_config: S3Config | None = None,
 ) -> str:
@@ -453,29 +455,39 @@ 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'.
 
     Returns:
         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"
 
-    # 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]
+    if forecast.ndim != 4:
+        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)
 
@@ -486,13 +498,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",
                 },
@@ -522,14 +541,17 @@ 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": (
                 f"Created "
                 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,
         },
     )