Add smaller fixes to make running from files more robust

CHANGELOG.md

@@ -6,6 +6,13 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+- Made minutes available for formatting file names [!15](https://github.com/dmidk/sunflow/pull/15), @JoachimKoenigslieb
+- Added `{seconds}` as a filename formatting variable [!15](https://github.com/dmidk/sunflow/pull/15), @JoachimKoenigslieb
+- Use `.expand_dims` instead of `.assign_coords` to make sure we have both time dimension and time coordinates when loading from files [!15](https://github.com/dmidk/sunflow/pull/15), @JoachimKoenigslieb
+- Added a `clearsky` config object which can load clear-sky data from files via `clearsky.path` or generate it with `pvlib` via the `simplified_solis` method [!15](https://github.com/dmidk/sunflow/pull/15), @JoachimKoenigslieb
+- `check_solar_elevation` now does not assume location is in Copenhagen by default [!15](https://github.com/dmidk/sunflow/pull/15), @JoachimKoenigslieb
+- Added `MIN_SOLAR_ELEVATION_DEGREES` (defaulting to 6 degrees) to configure the minimum maximum-corner solar elevation required to run [!15](https://github.com/dmidk/sunflow/pull/15), @JoachimKoenigslieb
+
 
 ### Added
 

README.md

@@ -108,6 +108,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 |
+| `MIN_SOLAR_ELEVATION_DEGREES` | `6` | Minimum maximum-corner solar elevation required to run |
 
 #### Data Source Configuration
 

config.yaml

@@ -22,6 +22,7 @@ KNMI:
   #   {month}: Two-digit month (e.g., 03)
   #   {day}: Two-digit day of month (e.g., 09)
   #   {hour}: Two-digit hour (e.g., 12)
+  #   {seconds}: Two-digit seconds (e.g., 00)
   # Path separators are supported, e.g.: "{year}/{month}/{day}/{dataset_name}_{timestamp}.nc"
 
 DWD:

sunflow/config.py

@@ -52,6 +52,7 @@ class NowcastConfig:
     max_waiting_time_minutes: int
     satellite_data_directory: str
     max_clearsky_fallback_days: int
+    min_solar_elevation_degrees: float
 
     @classmethod
     def from_env(cls) -> Self:
@@ -68,6 +69,7 @@ def from_env(cls) -> Self:
         - MAX_WAITING_TIME_MINUTES (default: 27)
         - SATELLITE_DATA_DIRECTORY (default: .)
         - MAX_CLEARSKY_FALLBACK_DAYS (default: 3)
+        - MIN_SOLAR_ELEVATION_DEGREES (default: 6)
         """
         return cls(
             nowcast_directory=os.getenv("NOWCAST_DIRECTORY", "."),
@@ -83,4 +85,7 @@ def from_env(cls) -> 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")),
+            min_solar_elevation_degrees=float(
+                os.getenv("MIN_SOLAR_ELEVATION_DEGREES", "6")
+            ),
         )

sunflow/data_io.py

@@ -133,6 +133,8 @@ def generate_input_filename(
         {month}: Two-digit month (e.g. 03)
         {day}: Two-digit day (e.g. 09)
         {hour}: Two-digit hour (e.g. 12)
+        {minute}: Two-digit minute (e.g. 30)
+        {seconds}: Two-digit seconds (e.g. 00)
 
     Path separators in the result are supported, so a format like
     ``{year}/{month}/{day}/{dataset_name}_{timestamp}.nc`` resolves to a
@@ -154,6 +156,8 @@ def generate_input_filename(
         month=time_step.strftime("%m"),
         day=time_step.strftime("%d"),
         hour=time_step.strftime("%H"),
+        minute=time_step.strftime("%M"),
+        seconds=time_step.strftime("%S"),
     )
     return filename
 
@@ -226,7 +230,15 @@ def load_data_from_files(
 
         try:
             ds = xr.open_dataset(filepath)
-            collected.append(ds.assign_coords(time=[time_step.replace(tzinfo=None)]))
+            time_step_naive = time_step.replace(tzinfo=None)
+
+            # Ensure time is in both dimension and in coords.
+            if "time" in ds.dims:
+                ds = ds.assign_coords(time=[time_step_naive])
+            else:
+                ds = ds.expand_dims(time=[time_step_naive])
+
+            collected.append(ds)
             logger.info(f"Loaded {data_type} from {filepath}")
         except Exception as e:
             logger.error(f"Failed to load {data_type} {filepath}: {e}")
@@ -337,6 +349,7 @@ def fetch_clearsky_with_fallback(
     time_steps: list[datetime],
     run_mode: str,
     max_fallback_days: int,
+    filename_format: str,
     config: dict[str, Any],
     bbox: str,
     dataset_name: str,
@@ -355,6 +368,7 @@ def fetch_clearsky_with_fallback(
         time_steps: Requested clearsky times (typically forecast times - 1 day).
         run_mode: One of 'download', 'files', or 's3'.
         max_fallback_days: Maximum number of days back to search per time step.
+        filename_format: Template string for clear-sky files.
         config: Dataset configuration dict.
         bbox: Bounding box string.
         dataset_name: Name of dataset (options: KNMI, DWD).
@@ -396,7 +410,7 @@ def fetch_clearsky_with_fallback(
                             bbox_choice,
                             nowcast_config.satellite_data_directory,
                             "clearsky data",
-                            config["filename_format"],
+                            filename_format,
                             bbox=bbox,
                         )
                     case "s3":
@@ -406,7 +420,7 @@ def fetch_clearsky_with_fallback(
                             bbox_choice,
                             s3_config,
                             "clearsky data",
-                            config["filename_format"],
+                            filename_format,
                             bbox=bbox,
                         )
 

sunflow/forecast.py

@@ -2,12 +2,73 @@
 from datetime import datetime
 
 import numpy as np
+import pvlib
 import xarray as xr
 from Models.ProbabilisticAdvection import ProbabilisticAdvection
 
 from .geospatial import get_coordinates
 
 
+def make_pvlib_clearsky_dataset(
+    times: list[datetime],
+    latitudes: np.ndarray,
+    longitudes: np.ndarray,
+    variable_name: str,
+) -> xr.Dataset:
+    """Compute a pvlib simplified-Solis clear-sky GHI dataset on a lat/lon grid.
+
+    For each timestep in `times`, solar position is computed at every grid
+    point by repeating the timestamp across all (lat, lon) pairs. The
+    simplified Solis clear-sky model is then applied to the apparent solar
+    elevation to produce global horizontal irradiance (GHI) in W m⁻².
+    Longitudes are normalised to the range [-180, 180] before the solar
+    position calculation.
+
+    Args:
+        times: Ordered list of timezone-aware datetimes for which to
+            compute clear-sky irradiance.
+        latitudes: 1-D array of latitude values in degrees North.
+        longitudes: 1-D array of longitude values in degrees East
+            (values > 180 are wrapped to [-180, 180]).
+
+    Returns:
+        xr.Dataset with a single variable keyed by
+        PVLIB_CLEARSKY_VARIABLE_NAME of shape (time, latitude, longitude)
+        containing GHI in W m⁻². The time coordinate is stored without
+        timezone information.
+    """
+    lon_grid, lat_grid = np.meshgrid(longitudes, latitudes)
+    lon_grid = np.where(lon_grid > 180, lon_grid - 360, lon_grid)
+    flat_latitudes = lat_grid.ravel()
+    flat_longitudes = lon_grid.ravel()
+
+    fields = []
+    for time in times:
+        repeated_times = [time] * len(flat_latitudes)
+        solar_position = pvlib.solarposition.get_solarposition(
+            repeated_times,
+            flat_latitudes,
+            flat_longitudes,
+            method="nrel_numpy",
+        )
+        clearsky = pvlib.clearsky.simplified_solis(solar_position["apparent_elevation"])
+        fields.append(clearsky["ghi"].to_numpy().reshape(lat_grid.shape))
+
+    return xr.Dataset(
+        {
+            variable_name: (
+                ["time", "latitude", "longitude"],
+                np.stack(fields),
+            )
+        },
+        coords={
+            "time": [time.replace(tzinfo=None) for time in times],
+            "latitude": latitudes,
+            "longitude": longitudes,
+        },
+    )
+
+
 def preprocess_data(
     data: xr.Dataset,
     time_steps: list[datetime],

sunflow/geospatial.py

@@ -83,6 +83,9 @@ def get_coordinates(ds: xr.Dataset) -> tuple[np.ndarray, np.ndarray]:
     elif "lat" in ds.coords and "lon" in ds.coords:
         latitudes = ds.lat.values
         longitudes = ds.lon.values
+    elif "latitude" in ds.coords and "longitude" in ds.coords:
+        latitudes = ds.latitude.values
+        longitudes = ds.longitude.values
     else:
         raise RuntimeError("Could not find coordinates in dataset.")
 
@@ -97,23 +100,26 @@ def get_coordinates(ds: xr.Dataset) -> tuple[np.ndarray, np.ndarray]:
 
 def check_solar_elevation(
     time: datetime,
-    lat: float = 55.6761,
-    lon: float = 12.5683,
+    lat: float,
+    lon: float,
 ) -> float:
     """Compute the solar elevation angle at a given time and location.
 
-    Uses pvlib to calculate the solar position. Defaults to Copenhagen,
-    Denmark (55.68°N, 12.57°E).
+    Uses pvlib to calculate the solar position. Longitudes in 0..360 convention
+    are converted to the -180..180 convention expected by pvlib.
 
     Args:
         time: Datetime (timezone-aware) for which to compute the elevation.
-        lat: Latitude in decimal degrees. Default 55.6761 (Copenhagen).
-        lon: Longitude in decimal degrees. Default 12.5683 (Copenhagen).
+        lat: Latitude in decimal degrees.
+        lon: Longitude in decimal degrees.
 
     Returns:
         Solar elevation angle in degrees above the horizon. Negative values
         indicate the sun is below the horizon.
     """
+    if lon > 180:
+        lon -= 360
+
     location = pvlib.location.Location(lat, lon)
     solar_elevation = location.get_solarposition(time)["elevation"].values[0]
     return solar_elevation

sunflow/main.py

@@ -22,8 +22,13 @@
     save_forecast,
 )
 from .downloaders import download_past_data
-from .forecast import multiply_clearsky, preprocess_data, simple_advection_forecast
-from .geospatial import check_solar_elevation, get_bbox
+from .forecast import (
+    make_pvlib_clearsky_dataset,
+    multiply_clearsky,
+    preprocess_data,
+    simple_advection_forecast,
+)
+from .geospatial import check_solar_elevation, get_bbox, get_coordinates
 from .time_handler import generate_time_steps, round_time
 from .validation import (
     MissingClearskyDataError,
@@ -193,6 +198,11 @@ def run_nowcast(
     """
     time_step_str = time_step.strftime("%Y-%m-%dT%H:%M:%SZ")
     logger.info(f"--- Running nowcast for {time_step_str} ---")
+    nc_variable_names = config["nc_variable_names"].copy()
+    clearsky_config = config.get(
+        "clearsky",
+        {"method": "file", "path": config["filename_format"]},
+    )
 
     # Fetch current data (with retry loop in operational mode)
     fetch_current_data_with_retry(
@@ -209,9 +219,14 @@ def run_nowcast(
 
     # Check solar elevation
     try:
-        solar_elevation = check_solar_elevation(time_step)
-        logger.info(f"Solar elevation: {solar_elevation:.2f} degrees")
-        if solar_elevation < 1:
+        lon_min, lat_min, lon_max, lat_max = map(float, bbox.split(","))
+        solar_elevation = max(
+            check_solar_elevation(time_step, lat=lat, lon=lon)
+            for lat in (lat_min, lat_max)
+            for lon in (lon_min, lon_max)
+        )
+        logger.info(f"Maximum corner solar elevation: {solar_elevation:.2f} degrees")
+        if solar_elevation < nowcast_config.min_solar_elevation_degrees:
             reason = "sun too low"
             logger.warning(f"{reason.capitalize()}. Skipping.\n")
             return RunResult(
@@ -261,10 +276,36 @@ def run_nowcast(
     if n_loaded == 0:
         raise RuntimeError("No past data loaded. Cannot proceed.")
 
+    if clearsky_config["method"] == "pvlib":
+        latitudes, longitudes = get_coordinates(data)
+        data = data.merge(
+            make_pvlib_clearsky_dataset(
+                past_time_steps,
+                latitudes,
+                longitudes,
+                nc_variable_names["sds_cs"],
+            )
+        )
+    elif clearsky_config["path"] != config["filename_format"]:
+        logger.info("Loading separate clearsky data for past observations...")
+        clearsky_data = fetch_clearsky_with_fallback(
+            past_time_steps,
+            run_mode,
+            nowcast_config.max_clearsky_fallback_days,
+            clearsky_config["path"],
+            config,
+            bbox,
+            dataset_name,
+            bbox_choice,
+            nowcast_config,
+            s3_config,
+        )
+        data = data.merge(clearsky_data[[nc_variable_names["sds_cs"]]])
+
     # Preprocess
     logger.info("Preprocessing data...")
     ratio_data, latitudes, longitudes = preprocess_data(
-        data, past_time_steps, config["nc_variable_names"]
+        data, past_time_steps, nc_variable_names
     )
 
     # Validate data shape
@@ -294,19 +335,29 @@ def run_nowcast(
     clearsky_t0_time = time_step - timedelta(days=1)
     all_clearsky_time_steps = [clearsky_t0_time] + previous_day_time_steps
 
-    # Fetch clearsky data with fallback to earlier days if a file is missing
-    logger.info("Fetching clearsky data...")
-    clearsky_data = fetch_clearsky_with_fallback(
-        all_clearsky_time_steps,
-        run_mode,
-        nowcast_config.max_clearsky_fallback_days,
-        config,
-        bbox,
-        dataset_name,
-        bbox_choice,
-        nowcast_config,
-        s3_config,
-    )
+    if clearsky_config["method"] == "pvlib":
+        logger.info("Generating pvlib clearsky data...")
+        clearsky_data = make_pvlib_clearsky_dataset(
+            all_clearsky_time_steps,
+            latitudes,
+            longitudes,
+            nc_variable_names["sds_cs"],
+        )
+    else:
+        # Fetch clearsky data with fallback to earlier days if a file is missing
+        logger.info("Fetching clearsky data...")
+        clearsky_data = fetch_clearsky_with_fallback(
+            all_clearsky_time_steps,
+            run_mode,
+            nowcast_config.max_clearsky_fallback_days,
+            clearsky_config["path"],
+            config,
+            bbox,
+            dataset_name,
+            bbox_choice,
+            nowcast_config,
+            s3_config,
+        )
 
     # Validate clearsky data
     try:
@@ -319,7 +370,7 @@ def run_nowcast(
         clearsky_data,
         previous_day_time_steps,
         expected_spatial_shape,
-        config["nc_variable_names"],
+        nc_variable_names,
     )
 
     # Multiply by clearsky to get actual solar irradiance
@@ -328,12 +379,12 @@ def run_nowcast(
         ratio_forecast,
         clearsky_data,
         previous_day_time_steps,
-        config["nc_variable_names"],
+        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"]
+        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)