wave energy converters and new wind and wave data modules

.gitignore

@@ -27,3 +27,4 @@ paper
 # Ignore IDE project files
 .idea/
 .vscode
+.vs

atlite/convert.py

@@ -23,6 +23,7 @@
 from dask.diagnostics import ProgressBar
 from numpy import pi
 from scipy.sparse import csr_matrix
+from tqdm import tqdm
 
 from atlite import csp as cspm
 from atlite import hydro as hydrom
@@ -36,6 +37,7 @@
 from atlite.resource import (
     get_cspinstallationconfig,
     get_solarpanelconfig,
+    get_waveenergyconverter,
     get_windturbineconfig,
     windturbine_smooth,
 )
@@ -653,7 +655,100 @@ def wind(
     )
 
 
-# irradiation
+# wave
+def convert_wave(ds, converter, time_chunk_size: int = 100) -> xr.DataArray:
+    r"""
+    Convert wave height (Hs) and wave peak period (Tp) data into normalized power output
+    using the device-specific Wave Energy Converter (WEC) power matrix.
+
+    This function matches each combination of significant wave height and peak period
+    in the dataset to a corresponding power output from the WEC power matrix.
+    The resulting power output is normalized by the maximum possible output (capacity)
+    to obtain the specific generation profile.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        Input dataset (cutout) containing two variables:
+        wave_height: significant wave height (m)
+        wave_period: peak wave period (s)
+    converter : dict
+        Dictionary defining the WEC characteristics, including:
+        Power_Matrix: a power matrix dictionary stored in `resources\wecgenerator`
+    time_chunk_size : int
+        Size of time chunks for processing large datasets, to limit memory spikes. Default is 100.
+
+    Returns
+    -------
+    xarray.DataArray
+        DataArray of specific power generation values (normalized power output).
+
+    Notes
+    -----
+    A progress message is printed every one million cases to track computation.
+    """
+    power_matrix = (
+        pd.DataFrame.from_dict(converter["Power_Matrix"])
+        .stack()
+        .rename_axis(index=["wave_height", "wave_period"])
+        .where(lambda x: x > 0)
+        .dropna()
+        .to_xarray()
+    )
+
+    results = []
+    steps = np.arange(0, len(ds.time), step=100)
+
+    for step in tqdm(
+        steps, desc="Processing wave data chunks", total=len(steps), unit="time chunk"
+    ):
+        ds_ = ds.isel(time=slice(step, step + time_chunk_size))
+        cf = power_matrix.interp(
+            {"wave_height": ds_.wave_height, "wave_period": ds_.wave_period},
+            method="nearest",
+        )
+        results.append(cf)
+
+    da = xr.concat(results, dim="time")
+    da.attrs["units"] = "kWh/kWp"
+    da = da.rename("specific generation")
+    da = da.fillna(0)
+
+    return da
+
+
+def wave(cutout, converter, **params):
+    """
+    Compute wave energy generation time series for a given cutout and Wave Energy Converter (WEC) type.
+
+    Parameters
+    ----------
+    cutout : atlite.Cutout
+        Atlite cutout object containing wave-related data (e.g., `wave_height`, `wave_period`).
+    wec_type : str, pathlib.Path, or dict
+        WEC configuration describing the device's power characteristics.
+
+    Returns
+    -------
+    xarray.DataArray
+        Time series of normalized wave power generation for the entire cutout area, with units of "kWh/kWp".
+        The dimensions and resolution follow the input cutout and aggregation parameters.
+
+    References
+    ----------
+    [1] Lavidas G., Mezilis L., Alday M., Baki H., Tan J., Jain A., Engelfried T. and Raghavan V.,
+    Marine renewables in Energy Systems: Impacts of climate data, generators, energy policies,
+    opportunities, and untapped potential for 100% decarbonised systems. Energy, Volume 336, 2025,
+    138359, ISSN 0360-5442, https://doi.org/10.1016/j.energy.2025.138359.
+    """
+    if isinstance(converter, str | Path):
+        converter = get_waveenergyconverter(converter)
+
+    return cutout.convert_and_aggregate(
+        convert_func=convert_wave, converter=converter, **params
+    )
+
+
 def convert_irradiation(
     ds,
     orientation,

atlite/cutout.py

@@ -44,6 +44,7 @@
     soil_temperature,
     solar_thermal,
     temperature,
+    wave,
     wind,
 )
 from atlite.data import available_features, cutout_prepare
@@ -673,6 +674,8 @@ def layout_from_capacity_list(self, data, col="Capacity"):
 
     wind = wind
 
+    wave = wave
+
     irradiation = irradiation
 
     pv = pv

atlite/datasets/__init__.py

@@ -6,6 +6,11 @@
 atlite datasets.
 """
 
-from atlite.datasets import era5, gebco, sarah
+from atlite.datasets import era5, gebco, mrel_wave, sarah
 
-modules = {"era5": era5, "sarah": sarah, "gebco": gebco}
+modules = {
+    "era5": era5,
+    "sarah": sarah,
+    "mrel_wave": mrel_wave,
+    "gebco": gebco,
+}

atlite/datasets/era5.py

@@ -57,6 +57,7 @@ def nullcontext():
     ],
     "temperature": ["temperature", "soil temperature", "dewpoint temperature"],
     "runoff": ["runoff"],
+    "wave": ["wave_height", "wave_period"],
 }
 
 static_features = {"height"}
@@ -244,6 +245,53 @@ def sanitize_runoff(ds):
     return ds
 
 
+def get_data_wave_height(retrieval_params):
+    """
+    Get wave height data for given retrieval parameters.
+    """
+    ds = retrieve_data(
+        variable=[
+            "significant_height_of_combined_wind_waves_and_swell",
+        ],
+        **retrieval_params,
+    )
+    ds = _rename_and_clean_coords(ds)
+    ds = ds.rename({"swh": "wave_height"})
+
+    return ds
+
+
+def sanitize_wave_height(ds):
+    """
+    Sanitize retrieved wave height data.
+    """
+    ds["wave_height"] = ds["wave_height"].clip(min=0.0)
+    return ds
+
+
+def get_data_wave_period(retrieval_params):
+    """
+    Get wave period data for given retrieval parameters.
+    """
+    ds = retrieve_data(
+        variable=["peak_wave_period"],
+        **retrieval_params,
+    )
+
+    ds = _rename_and_clean_coords(ds)
+    ds = ds.rename({"pp1d": "wave_period"})
+
+    return ds
+
+
+def sanitize_wave_period(ds):
+    """
+    Sanitize retrieved wave period data.
+    """
+    ds["wave_period"] = ds["wave_period"].clip(min=0.0)
+    return ds
+
+
 def get_data_height(retrieval_params):
     """
     Get height data for given retrieval parameters.

atlite/datasets/mrel_wave.py

@@ -0,0 +1,116 @@
+# SPDX-FileCopyrightText: Contributors to atlite <https://github.com/pypsa/atlite>
+#
+# SPDX-License-Identifier: MIT
+"""
+Module for curating the already downloaded wave data of MREL (ECHOWAVE).
+
+For further reference see:
+[1] Matías A., George L., The ECHOWAVE Hindcast: A 30-years high resolution database
+for wave energy applications in North Atlantic European waters, Renewable Energy,
+Volume 236, 2024, 121391,ISSN 0960-1481, https://doi.org/10.1016/j.renene.2024.121391
+"""
+
+import logging
+
+import numpy as np
+import xarray as xr
+from rasterio.warp import Resampling
+
+from atlite.gis import regrid
+
+logger = logging.getLogger(__name__)
+
+crs = 4326
+dx = 0.03
+dy = 0.03
+
+features = {"hs": "wave_height", "fp": "wave_period"}
+
+
+def _rename_and_clean_coords(ds, cutout):
+    """
+    Rename 'longitude' and 'latitude' columns to 'x' and 'y', fix roundings and grid dimensions.
+    """
+    coords = cutout.coords
+
+    if "longitude" in ds and "latitude" in ds:
+        ds = ds.rename({"longitude": "x", "latitude": "y"})
+    # round coords since cds coords are float32 which would lead to mismatches
+    ds = ds.assign_coords(
+        x=np.round(ds.x.astype(float), 5), y=np.round(ds.y.astype(float), 5)
+    )
+    if (cutout.dx != dx) or (cutout.dy != dy):
+        ds = regrid(ds, coords["x"], coords["y"], resampling=Resampling.average)
+
+    return ds
+
+
+def sanitize_wave_height(ds):
+    """
+    Sanitize retrieved wave height data.
+    """
+    ds["wave_height"] = ds["wave_height"].clip(min=0.0)
+    return ds
+
+
+def sanitize_wave_period(ds):
+    """
+    Sanitize retrieved wave height data.
+    """
+    ds["wave_period"] = ds["wave_period"].clip(min=0.0)
+    return ds
+
+
+def _bounds(coords, pad: float = 0) -> dict[str, slice]:
+    """
+    Convert coordinate bounds to slice and pad if requested.
+    """
+    x0, x1 = coords["x"].min().item() - pad, coords["x"].max().item() + pad
+    y0, y1 = coords["y"].min().item() - pad, coords["y"].max().item() + pad
+
+    return {"x": slice(x0, x1), "y": slice(y0, y1)}
+
+
+def get_data(cutout, feature, tmpdir, **creation_parameters):
+    """
+    Load stored MREL (ECHOWAVE) data and reformat to matching the given cutout.
+
+    This function loads and resamples the stored MREL data for a given
+    `atlite.Cutout`.
+
+    Parameters
+    ----------
+    cutout : atlite.Cutout
+    feature : str
+        Name of the feature data to retrieve. Must be in
+        `atlite.datasets.mrel_wave.features`
+    **creation_parameters :
+        Mandatory arguments are:
+            * 'data_path', str. Directory of the stored MREL data.
+
+    Returns
+    -------
+    xarray.Dataset
+        Dataset of dask arrays of the retrieved variables.
+    """
+
+    if "data_path" not in creation_parameters:
+        logger.error('Argument "data_path" not defined')
+        raise ValueError('Argument "data_path" not defined')
+    path = creation_parameters["data_path"]
+
+    ds = xr.open_dataset(path)
+    ds = _rename_and_clean_coords(ds, cutout)
+    bounds = _bounds(cutout.coords, pad=creation_parameters.get("pad", 0))
+    ds = ds.sel(**bounds)
+
+    # invert the wave peak frequency to obrain wave peak period
+    ds["tp"] = 1 / ds["fp"]
+
+    ds = ds[list(features.keys())].rename(features)
+    for feature in features.values():
+        sanitize_func = globals().get(f"sanitize_{feature}")
+        if sanitize_func is not None:
+            ds = sanitize_func(ds)
+
+    return ds

atlite/resource.py

@@ -31,6 +31,7 @@
 WINDTURBINE_DIRECTORY = RESOURCE_DIRECTORY / "windturbine"
 SOLARPANEL_DIRECTORY = RESOURCE_DIRECTORY / "solarpanel"
 CSPINSTALLATION_DIRECTORY = RESOURCE_DIRECTORY / "cspinstallation"
+WAVEENERGYCONVERTER_DIRECTORY = RESOURCE_DIRECTORY / "waveenergyconverter"
 
 if TYPE_CHECKING:
     from typing import TypedDict
@@ -109,6 +110,26 @@ def get_windturbineconfig(
     return _validate_turbine_config_dict(conf, add_cutout_windspeed)
 
 
+def get_waveenergyconverter(converter):
+    """
+    Load the wave energy converter power matrix
+    the configuration can either be one from local storage then 'wec_type' is
+    considered part of the file base name '<wec_type>.yaml'
+    """
+    assert isinstance(converter, (str | Path))
+
+    if isinstance(converter, str):
+        converter_path = waveenergyconverter[converter.replace(".yaml", "")]
+
+    elif isinstance(converter, Path):
+        converter_path = converter
+
+    with open(converter_path) as f:
+        conf = yaml.safe_load(f)
+
+    return conf
+
+
 def get_solarpanelconfig(panel):
     """
     Load the 'panel'.yaml file from local disk and provide a solar panel dict.
@@ -512,6 +533,9 @@ def get_oedb_windturbineconfig(
 # Global caches
 _oedb_turbines = None
 windturbines = arrowdict({p.stem: p for p in WINDTURBINE_DIRECTORY.glob("*.yaml")})
+waveenergyconverter = arrowdict(
+    {p.stem: p for p in WAVEENERGYCONVERTER_DIRECTORY.glob("*.yaml")}
+)
 solarpanels = arrowdict({p.stem: p for p in SOLARPANEL_DIRECTORY.glob("*.yaml")})
 cspinstallations = arrowdict(
     {p.stem: p for p in CSPINSTALLATION_DIRECTORY.glob("*.yaml")}