diff --git a/atlite/convert.py b/atlite/convert.py index c165cfa3..585346fe 100644 --- a/atlite/convert.py +++ b/atlite/convert.py @@ -478,19 +478,28 @@ def solar_thermal( # wind -def convert_wind(ds, turbine): +def convert_wind( + ds: xr.Dataset, turbine, windspeed_bias_correction=None, from_height=None +): """ Convert wind speeds for turbine to wind energy generation. """ V, POW, hub_height, P = itemgetter("V", "POW", "hub_height", "P")(turbine) - wnd_hub = windm.extrapolate_wind_speed(ds, to_height=hub_height) + if windspeed_bias_correction is not None: + ds = ds.assign( + {f"wnd{from_height}m": ds[f"wnd{from_height}m"] * windspeed_bias_correction} + ) + + wnd_hub = windm.extrapolate_wind_speed( + ds, to_height=hub_height, from_height=from_height + ) - def _interpolate(da): + def apply_power_curve(da): return np.interp(da, V, POW / P) da = xr.apply_ufunc( - _interpolate, + apply_power_curve, wnd_hub, input_core_dims=[[]], output_core_dims=[[]], @@ -503,7 +512,15 @@ def _interpolate(da): return da -def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params): +def wind( + cutout, + turbine, + smooth=False, + add_cutout_windspeed=False, + real_long_run_average_windspeed=None, + real_long_run_average_height=100, + **params, +): """ Generate wind generation time-series. @@ -529,6 +546,10 @@ def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params): output at the highest wind speed in the power curve. If False, a warning will be raised if the power curve does not have a cut-out wind speed. The default is False. + real_long_run_average_speed : Path or rasterio dataset or None + Raster dataset with wind speeds to bias correct long run average wind speeds + real_long_run_average_height : int = 100 + Height in meters of provided real_long_average_speed dataset Note ---- @@ -549,8 +570,24 @@ def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params): if smooth: turbine = windturbine_smooth(turbine, params=smooth) + windspeed_bias_correction = None + from_height = None + + if real_long_run_average_windspeed is not None: + windspeed_bias_correction = windm.calculate_windspeed_bias_correction( + cutout.data, + real_long_run_average_windspeed, + real_long_run_average_height, + crs=cutout.crs, + ) + from_height = real_long_run_average_height + return cutout.convert_and_aggregate( - convert_func=convert_wind, turbine=turbine, **params + convert_func=convert_wind, + turbine=turbine, + windspeed_bias_correction=windspeed_bias_correction, + from_height=from_height, + **params, ) diff --git a/atlite/datasets/era5.py b/atlite/datasets/era5.py index 786e31b4..eb47aa01 100644 --- a/atlite/datasets/era5.py +++ b/atlite/datasets/era5.py @@ -8,6 +8,7 @@ https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation """ +import datetime import logging import os import warnings @@ -45,6 +46,7 @@ def nullcontext(): features = { "height": ["height"], "wind": ["wnd100m", "wnd_azimuth", "roughness"], + "wind_lra": ["wnd100m_lra"], "influx": [ "influx_toa", "influx_direct", @@ -58,6 +60,7 @@ def nullcontext(): } static_features = {"height"} +longrunaverage_features = {"wind_lra"} def _add_height(ds): @@ -101,6 +104,32 @@ def _rename_and_clean_coords(ds, add_lon_lat=True): return ds +def get_data_wind_lra(retrieval_params): + """ + Get long run average wind data for given retrieval parameters. + """ + ds = retrieve_data( + variable=[ + "100m_u_component_of_wind", + "100m_v_component_of_wind", + ], + **retrieval_params, + ) + ds = _rename_and_clean_coords(ds) + + ds["wnd100m_lra"] = ( + sqrt(ds["u100"] ** 2 + ds["v100"] ** 2) + .mean("date") + .assign_attrs( + units=ds["u100"].attrs["units"], + long_name="100 metre wind speed as long run average", + ) + ) + ds = ds.drop_vars(["u100", "v100"]) + + return ds + + def get_data_wind(retrieval_params): """ Get wind data for given retrieval parameters. @@ -255,7 +284,7 @@ def _area(coords): return [y1, x0, y0, x1] -def retrieval_times(coords, static=False, monthly_requests=False): +def retrieval_times(coords, static=False, monthly_requests=False, longrunaverage=False): """ Get list of retrieval cdsapi arguments for time dimension in coordinates. @@ -287,6 +316,12 @@ def retrieval_times(coords, static=False, monthly_requests=False): "day": str(time[0].day), "time": time[0].strftime("%H:00"), } + elif longrunaverage: + return { + "year": [str(y) for y in range(1980, datetime.date.today().year)], + "month": [f"{m:02}" for m in range(1, 12 + 1)], + "time": ["00:00"], + } # Prepare request for all months and years times = [] @@ -323,7 +358,7 @@ def noisy_unlink(path): logger.error(f"Unable to delete file {path}, as it is still in use.") -def retrieve_data(product, chunks=None, tmpdir=None, lock=None, **updates): +def retrieve_data(dataset, chunks=None, tmpdir=None, lock=None, **updates): """ Download data like ERA5 from the Climate Data Store (CDS). @@ -340,7 +375,7 @@ def retrieve_data(product, chunks=None, tmpdir=None, lock=None, **updates): client = cdsapi.Client( info_callback=logger.debug, debug=logging.DEBUG >= logging.root.level ) - result = client.retrieve(product, request) + result = client.retrieve(dataset, request) if lock is None: lock = nullcontext() @@ -418,7 +453,8 @@ def get_data( sanitize = creation_parameters.get("sanitize", True) retrieval_params = { - "product": "reanalysis-era5-single-levels", + "dataset": "reanalysis-era5-single-levels", + "product_type": "reanalysis", "area": _area(coords), "chunks": cutout.chunks, "grid": [cutout.dx, cutout.dy], @@ -431,7 +467,7 @@ def get_data( logger.info(f"Requesting data for feature {feature}...") - def retrieve_once(time): + def retrieve_once(time, longrunaverage=False): ds = func({**retrieval_params, **time}) if sanitize and sanitize_func is not None: ds = sanitize_func(ds) @@ -439,6 +475,10 @@ def retrieve_once(time): if feature in static_features: return retrieve_once(retrieval_times(coords, static=True)).squeeze() + elif feature in longrunaverage_features: + retrieval_params["dataset"] += "-monthly-means" + retrieval_params["product_type"] = "monthly_averaged_reanalysis" + return retrieve_once(retrieval_times(coords, longrunaverage=True)) time_chunks = retrieval_times(coords, monthly_requests=monthly_requests) if concurrent_requests: diff --git a/atlite/wind.py b/atlite/wind.py index 07f17d51..6aa2d7a8 100644 --- a/atlite/wind.py +++ b/atlite/wind.py @@ -7,8 +7,13 @@ import logging import re +from pathlib import Path import numpy as np +import rasterio as rio +import xarray as xr + +from .gis import _as_transform logger = logging.getLogger(__name__) @@ -51,14 +56,11 @@ def extrapolate_wind_speed(ds, to_height, from_height=None): Retrieved 2019-02-15. """ - # Fast lane - to_name = f"wnd{int(to_height):0d}m" - if to_name in ds: - return ds[to_name] - if from_height is None: # Determine closest height to to_name - heights = np.asarray([int(s[3:-1]) for s in ds if re.match(r"wnd\d+m", s)]) + heights = np.asarray( + [int(m.group(1)) for s in ds if (m := re.match(r"wnd(\d+)m", s))] + ) if len(heights) == 0: raise AssertionError("Wind speed is not in dataset") @@ -67,17 +69,46 @@ def extrapolate_wind_speed(ds, to_height, from_height=None): from_name = f"wnd{int(from_height):0d}m" + # Fast lane + if from_height == to_height: + return ds[from_name] + # Wind speed extrapolation wnd_spd = ds[from_name] * ( np.log(to_height / ds["roughness"]) / np.log(from_height / ds["roughness"]) ) - wnd_spd.attrs.update( + return wnd_spd.assign_attrs( { "long name": f"extrapolated {to_height} m wind speed using logarithmic " f"method with roughness and {from_height} m wind speed", "units": "m s**-1", } - ) + ).rename(f"wnd{to_height}m") + + +def calculate_windspeed_bias_correction(ds, real_lra, lra_height, crs): + data_lra = ds[f"wnd{lra_height}m_lra"] + + if isinstance(real_lra, (str, Path)): + real_lra = rio.open(real_lra) + + if isinstance(real_lra, rio.DatasetReader): + real_lra = rio.band(real_lra, 1) + + if isinstance(real_lra, rio.Band): + real_lra, transform = rio.warp.reproject( + real_lra, + np.empty(data_lra.shape), + dst_crs=crs, + dst_transform=_as_transform( + x=data_lra.indexes["x"], y=data_lra.indexes["y"] + ), + resampling=rio.enums.Resampling.average, + ) + + real_lra = xr.DataArray( + real_lra, [data_lra.indexes["y"], data_lra.indexes["x"]] + ) - return wnd_spd.rename(to_name) + return real_lra / data_lra