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

@@ -36,6 +36,7 @@
 from atlite.resource import (
     get_cspinstallationconfig,
     get_solarpanelconfig,
+    get_wecgeneratorconfig,
     get_windturbineconfig,
     windturbine_smooth,
 )
@@ -653,6 +654,67 @@ def wind(
     )
 
 
+# #wave
+def convert_wave(ds, wec_type):
+    power_matrix = pd.DataFrame.from_dict(wec_type["Power_Matrix"])
+
+    max_pow = power_matrix.to_numpy().max()
+
+    Hs = np.ceil(ds["wave_height"] * 2) / 2
+    Tp = np.ceil(ds["wave_period"] * 2) / 2
+
+    Hs_list = Hs.to_numpy().flatten().tolist()
+    Tp_list = Tp.to_numpy().flatten().tolist()
+
+    # empty list for result
+    power_list = []
+    cases = len(Hs_list)
+    count = 0
+
+    # for loop to loop through Hs and Tp pairs and get the power output and capacity factor
+    for Hs_ind, Tp_ind in zip(Hs_list, Tp_list):
+        if count % 1000000 == 0:
+            print(f"Case {count} of {cases}: {count / cases * 100} %")
+
+        if np.isnan(Hs_ind) or np.isnan(Tp_ind):
+            power_list.append(0)
+        elif Hs_ind > 10 or Tp_ind > 18:
+            power_list.append(0)
+        else:
+            generated_power = power_matrix.loc[Hs_ind, Tp_ind]
+            power_list.append(generated_power / max_pow)
+        count += 1
+
+    # results list to numpy array
+    power_list_np = np.array(power_list)
+
+    power_list_np = power_list_np.reshape(Hs.shape)
+
+    da = xr.DataArray(
+        power_list_np, coords=Hs.coords, dims=Hs.dims, name="Power generated"
+    )
+    da.attrs["units"] = "kWh/kWp"
+    da = da.rename("specific generation")
+    da = da.fillna(0)
+
+    return da
+
+
+def wave(cutout, wec_type, **params):
+    """
+    Generate wave generation time series
+
+    evaluates the significant wave height (Hs) and wave peak period (Tp)
+    and assesses the power output with the chosen power matrix for each time step and grid cell
+    """
+    if isinstance(wec_type, (str, Path)):
+        wec_type = get_wecgeneratorconfig(wec_type)
+
+    return cutout.convert_and_aggregate(
+        convert_func=convert_wave, wec_type=wec_type, **params
+    )
+
+
 # irradiation
 def convert_irradiation(
     ds,

atlite/cutout.py

@@ -44,6 +44,7 @@
     soil_temperature,
     solar_thermal,
     temperature,
+    wave,
     wind,
 )
 from atlite.data import available_features, cutout_prepare
@@ -661,6 +662,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,12 @@
 atlite datasets.
 """
 
-from atlite.datasets import era5, gebco, sarah
+from atlite.datasets import cerra, era5, gebco, mrel_wave, sarah
 
-modules = {"era5": era5, "sarah": sarah, "gebco": gebco}
+modules = {
+    "era5": era5,
+    "sarah": sarah,
+    "mrel_wave": mrel_wave,
+    "cerra": cerra,
+    "gebco": gebco,
+}

atlite/datasets/cerra.py

@@ -0,0 +1,68 @@
+"""
+In order to create a CERRA cutout, the data must be manually downloaded from the Climate Data Store.
+The variable used is "10m wind speed" and there is not a direction component in it.
+This 10m wind speed was transformed into a 100m wind speed in order to follow the rest of atlite's processes.
+"""
+
+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.05
+dy = 0.05
+
+features = {"wind": ["wnd100m", "roughness"]}
+
+
+def as_slice(bounds, pad=True):
+    """
+    Convert coordinate bounds to slice and pad by 0.01.
+    """
+    if not isinstance(bounds, slice):
+        bounds = bounds + (-0.01, 0.01)
+        bounds = slice(*bounds)
+    return bounds
+
+
+def get_data(cutout, feature, tmpdir, **creation_parameters):
+    """
+    Retrieve data from a local CERRA dataset and process it.
+    """
+    coords = cutout.coords
+
+    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 = ds.sel(x=as_slice(cutout.extent[:2]), y=as_slice(cutout.extent[2:]))
+    ds = ds.assign_coords(x=ds.x.astype(float).round(4), y=ds.y.astype(float).round(4))
+
+    if (cutout.dx != dx) or (cutout.dy != dy):
+        ds = regrid(ds, coords["x"], coords["y"], resampling=Resampling.average)
+
+    if "sr" in ds:
+        ds = ds.rename({"sr": "roughness"})
+
+    logger.info("Calculating 100 metre wind speed")
+    if "si10" in ds and "roughness" in ds:
+        ds["wnd100m"] = (
+            ds["si10"] * (np.log(100 / ds["roughness"]) / np.log(10 / ds["roughness"]))
+        ).assign_attrs(units="m s**-1", long_name="100 metre wind speed")
+        ds = ds.drop_vars("si10")
+
+    ds = ds.assign_coords(x=ds.coords["x"], y=ds.coords["y"])
+
+    logger.info("Resampling to 1H.")
+    ds = ds.resample(time="1h").interpolate("linear")
+
+    return ds

atlite/datasets/era5.py

@@ -44,7 +44,7 @@ def nullcontext():
 
 features = {
     "height": ["height"],
-    "wind": ["wnd100m", "wnd_shear_exp", "wnd_azimuth", "roughness"],
+    "wind": ["wnd100m", "wnd_azimuth", "roughness"],
     "influx": [
         "influx_toa",
         "influx_direct",
@@ -55,6 +55,8 @@ def nullcontext():
     ],
     "temperature": ["temperature", "soil temperature", "dewpoint temperature"],
     "runoff": ["runoff"],
+    "wave_height": ["wave_height"],
+    "wave_period": ["wave_period"],
 }
 
 static_features = {"height"}
@@ -244,6 +246,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,108 @@
+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.0625
+dy = 0.04
+
+features = {"wave_height": ["wave_height"], "wave_period": ["wave_period"]}
+
+
+def _rename_and_clean_coords(ds):
+    """
+    Rename 'longitude' and 'latitude' columns to 'x' and 'y' and fix roundings.
+
+    Optionally (add_lon_lat, default:True) preserves latitude and
+    longitude columns as 'lat' and 'lon'.
+    """
+    ds = ds.rename({"longitude": "x", "latitude": "y"})
+
+    ds = ds.assign_coords(
+        x=np.round(ds.x.astype(float), 5), y=np.round(ds.y.astype(float), 5)
+    )
+
+    return ds
+
+
+def get_data_wave_height(ds):
+    ds = ds.rename({"hs": "wave_height"})
+    ds["wave_height"] = ds["wave_height"].clip(min=0.0)
+
+    return ds
+
+
+def get_data_wave_period(ds):
+    ds = ds.rename({"tp": "wave_period"})
+    # ds["wave_period"] = (1 / ds["wave_period"])
+    ds["wave_period"] = ds["wave_period"].clip(min=0.0)
+
+    return ds
+
+
+def as_slice(bounds, pad=True):
+    """
+    Convert coordinate bounds to slice and pad by 0.01.
+    """
+    if not isinstance(bounds, slice):
+        bounds = bounds + (-0.01, 0.01)
+        bounds = slice(*bounds)
+    return bounds
+
+
+def get_data(cutout, feature, tmpdir, **creation_parameters):
+    coords = cutout.coords
+
+    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)
+
+    if "longitude" in ds and "latitude" in ds:
+        ds = ds.rename({"longitude": "x", "latitude": "y"})
+
+    ds = ds.sel(x=as_slice(cutout.extent[:2]), y=as_slice(cutout.extent[2:]))
+    ds = ds.assign_coords(x=ds.x.astype(float).round(4), y=ds.y.astype(float).round(4))
+
+    if (cutout.dx != dx) or (cutout.dy != dy):
+        ds = regrid(ds, coords["x"], coords["y"], resampling=Resampling.average)
+
+    # coords = cutout.coords
+
+    # if "data_path" not in creation_parameters:
+    #     logger.error('Argument "data_path" not defined')
+    #     return None
+
+    # path = creation_parameters["data_path"]
+
+    # logger.info(f"Opening dataset from {path}")
+    # ds = xr.open_dataset(path, chunks=cutout.chunks)
+    # ds = _rename_and_clean_coords(ds)
+
+    variables = ds.data_vars
+
+    for var in variables:
+        if var not in ["hs", "tp"]:
+            ds = ds.drop_vars(var)
+
+    # ds = ds.sel(x=as_slice(cutout.extent[:2]), y=as_slice(cutout.extent[2:]))
+
+    # if (cutout.dx != dx) or (cutout.dy != dy):
+    #     ds = regrid(ds, coords["x"], coords["y"], resampling=Resampling.average)
+
+    logger.info("Obtaining wave data.")
+
+    ds = get_data_wave_height(ds)
+    ds = get_data_wave_period(ds)
+
+    # ds = ds.assign_coords(x=ds.coords["x"], y=ds.coords["y"])
+
+    return ds