wave energy converters and new wind and wave data modules

.gitignore

@@ -29,4 +29,4 @@ paper
 .idea/
 .vscode
 AGENTS.md
-uv.lock
+uv.lock

atlite/convert.py

@@ -21,6 +21,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
@@ -34,6 +35,7 @@
 from atlite.resource import (
     get_cspinstallationconfig,
     get_solarpanelconfig,
+    get_waveenergyconverter,
     get_windturbineconfig,
     windturbine_smooth,
 )
@@ -1098,7 +1100,103 @@ def wind(
     )
 
 
-# irradiation
+def convert_wave(
+    ds: xr.Dataset,
+    converter: dict,
+    time_chunk_size: int = 100,
+) -> xr.DataArray:
+    """
+    Convert wave height and peak period data into normalized WEC power output.
+
+    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")
+
+    return da.fillna(0)
+
+
+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: xr.Dataset,
     orientation: Callable,

atlite/cutout.py

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

atlite/datasets/__init__.py

@@ -4,12 +4,13 @@
 
 """atlite datasets."""
 
-from atlite.datasets import era5, era5_edh, gebco, glofas, sarah
+from atlite.datasets import era5, era5_edh, gebco, glofas, mrel_wave, sarah
 
 modules = {
     "era5": era5,
     "era5-edh": era5_edh,
     "gebco": gebco,
     "glofas": glofas,
     "sarah": sarah,
+    "mrel_wave": mrel_wave,
 }

atlite/datasets/era5.py

@@ -70,6 +70,7 @@ def nullcontext():  # noqa: D103
     ],
     "temperature": ["temperature", "soil temperature", "dewpoint temperature"],
     "runoff": ["runoff"],
+    "wave": ["wave_height", "wave_period"],
 }
 
 static_features = {"height"}
@@ -384,6 +385,93 @@ def sanitize_runoff(ds: xr.Dataset) -> xr.Dataset:
     return ds
 
 
+def get_data_wave_height(retrieval_params):
+    """
+    Retrieve and compute wave height variable from ERA5.
+
+    Parameters
+    ----------
+    retrieval_params
+        CDS API retrieval parameters including area, time, and format.
+
+    Returns
+    -------
+    xr.Dataset
+        Dataset with variable: wave_height.
+    """
+    ds = retrieve_data(
+        variable=[
+            "significant_height_of_combined_wind_waves_and_swell",
+        ],
+        **retrieval_params,
+    )
+    ds = _rename_and_clean_coords(ds)
+
+    return ds.rename({"swh": "wave_height"})
+
+
+def sanitize_wave_height(ds):
+    """
+    Clip wave height values to 0.0.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        Dataset containing the 'wave_height' variable.
+
+    Returns
+    -------
+    xr.Dataset
+        Dataset with negative 'wave_height' values clipped to zero.
+    """
+    ds["wave_height"] = ds["wave_height"].clip(min=0.0)
+
+    return ds
+
+
+def get_data_wave_period(retrieval_params):
+    """
+    Retrieve and compute wave period variable from ERA5.
+
+    Parameters
+    ----------
+    retrieval_params
+        CDS API retrieval parameters including area, time, and format.
+
+    Returns
+    -------
+    xr.Dataset
+        Dataset with variable: wave_period.
+    """
+    ds = retrieve_data(
+        variable=["peak_wave_period"],
+        **retrieval_params,
+    )
+
+    ds = _rename_and_clean_coords(ds)
+
+    return ds.rename({"pp1d": "wave_period"})
+
+
+def sanitize_wave_period(ds):
+    """
+    Clip wave period values to 0.0.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        Dataset containing the 'wave_period' variable.
+
+    Returns
+    -------
+    xr.Dataset
+        Dataset with negative 'wave_period' values clipped to zero.
+    """
+    ds["wave_period"] = ds["wave_period"].clip(min=0.0)
+
+    return ds
+
+
 def get_data_height(retrieval_params: dict[str, Any]) -> xr.Dataset:
     """
     Retrieve geopotential and convert to terrain height.