diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a51720..3c6c62e 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index da9dc8f..3e3143d 100644 --- a/README.md +++ b/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 diff --git a/config.yaml b/config.yaml index da1625d..f5d0c05 100644 --- a/config.yaml +++ b/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: diff --git a/sunflow/config.py b/sunflow/config.py index 6cb2e7f..ddddcf8 100644 --- a/sunflow/config.py +++ b/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") + ), ) diff --git a/sunflow/data_io.py b/sunflow/data_io.py index 70973d3..b2c8c43 100644 --- a/sunflow/data_io.py +++ b/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, ) diff --git a/sunflow/forecast.py b/sunflow/forecast.py index 3abe325..d164819 100644 --- a/sunflow/forecast.py +++ b/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], diff --git a/sunflow/geospatial.py b/sunflow/geospatial.py index 80712ff..5c1dc29 100644 --- a/sunflow/geospatial.py +++ b/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 diff --git a/sunflow/main.py b/sunflow/main.py index 6034bab..4a88357 100644 --- a/sunflow/main.py +++ b/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)