# (c) 2026 DTU Wind Energy
"""
Utility functions for working with wind data
"""
__all__ = [
"wind_speed",
"wind_direction",
"wind_speed_and_direction",
"wind_vectors",
"wind_direction_difference",
"wd_to_sector",
"vinterp_wind_direction",
"vinterp_wind_speed",
"rotor_equivalent_wind_speed",
"tswc_resample",
]
import numbers
import numpy as np
import pandas as pd
import xarray as xr
from scipy.integrate import trapezoid
from windkit.xarray_structures.sector import create_sector_coords
from windkit.spatial import _spatial_stack, _spatial_unstack
[docs]
def wind_speed(u, v):
"""
Calculate wind speed from wind vectors.
Parameters
----------
u, v : numpy.ndarray, xarray.DataArray
U and V wind vectors.
Returns
-------
ws : numpy.ndarray, xarray.DataArray
Wind speed.
"""
return np.sqrt(u * u + v * v)
[docs]
def wind_direction(u, v):
"""
Calculate wind directions from wind vectors.
Parameters
----------
u, v : np.ndarray, xr.DataArray
U and V wind vectors.
Returns
-------
wd : np.ndarray, xr.DataArray
Wind direction
"""
return 180.0 + np.arctan2(u, v) * 180.0 / np.pi
[docs]
def wind_speed_and_direction(u, v):
"""
Calculate wind speed and wind direction from wind vectors.
Parameters
----------
u, v : numpy.ndarray, xarray.DataArray
U and V wind vectors.
Returns
-------
speed : numpy.ndarray, xarray.DataArray
Wind speed.
direction : numpy.ndarray, xarray.DataArray
Wind direction.
"""
return wind_speed(u, v), wind_direction(u, v)
[docs]
def wind_vectors(ws, wd):
"""
Calculate wind vectors u,v from the speed and direction.
Parameters
----------
speed : numpy.ndarray, xarray.DataArray
Wind speed
direction : numpy.ndarray, xarray.DataArray
Wind direction
Returns
-------
u, v : numpy.ndarray, xarray.DataArray
Wind vectors u and v
"""
return (
-np.abs(ws) * np.sin(np.pi / 180.0 * wd),
-np.abs(ws) * np.cos(np.pi / 180.0 * wd),
)
[docs]
def wind_direction_difference(wd_obs, wd_mod):
"""
Calculate the circular (minimum) distance between
two directions (observed and modelled).
Parameters
----------
wd_obs : xarray.DataArray
observed direction arrays.
wd_mod: xarray.DataArray
modelled direction arrays.
Returns
-------
xarray.DataArray: circular (minimum) differences.
Examples
--------
>>> wd_obs = xr.DataArray([15.0, 345.0, 355.0], dims=('time',))
>>> wd_mod = xr.DataArray([345.0, 300.0, 5.0], dims=('time',))
>>> wind_direction_difference(wd_obs, wd_mod)
<xarray.DataArray (time: 3)>
array([-30., -45., 10.])
Dimensions without coordinates: time
"""
wd_diff = wd_mod - wd_obs
wd_diff = wd_diff.where(wd_diff < 180.0, wd_diff - 360.0)
wd_diff = wd_diff.where(wd_diff > -180.0, wd_diff + 360.0)
return wd_diff
[docs]
def wd_to_sector(wd, sectors=12, output_type="centers", quantiles=False):
"""
Convert wind directions to 0-based sector indices.
Parameters
----------
wd : xarray.DataArray, numpy.array
Wind directions. The function uses xarray.apply_ufunc, so the return value
will keep the shape of the input value.
sectors : int
Number of sectors. Defaults to 12.
output_type : str
If set to 'centers' the values in 'wd' are the sector centers. If set to
'indices', the values in 'wd' are the sector indices. Defaults to 'centers'.
quantiles : bool
Allows to use equal probability sectors (quantiles=True) instead of fixed
width sectors. Note that this is an experimental feature to be used only
together with the :py:mod:`windkit.ltc` module for now. Other :py:mod:`windkit` modules may
not be compatible with non fixed width sectors. Defaults to False.
Returns
-------
sector_centers : xarray.DataArray,np.array
wind speed sector centers.
sector_coords : xarray.DataArray
data array with sector coordinates incling center, ceiling and floor.
Examples
--------
>>> wd = xr.DataArray([355.0, 14.0, 25.0, 270.0,], dims=('time',))
>>> wd_to_sector(wd)
(<xarray.DataArray (time: 4)>
array([ 0., 0., 30., 270.])
Dimensions without coordinates: time,
<xarray.DataArray (sector: 12)>
array([ 0., 30., 60., 90., 120., 150., 180., 210., 240., 270., 300.,
330.])
Coordinates:
* sector (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
sector_ceil (sector) float64 15.0 45.0 75.0 105.0 ... 285.0 315.0 345.0
sector_floor (sector) float64 345.0 15.0 45.0 75.0 ... 255.0 285.0 315.0)
"""
def _wd_to_sector_constant(wd, n_sectors=12):
width = 360.0 / n_sectors
edges = np.linspace(0.0, 360.0, n_sectors + 1)
edges[0] = -0.1
edges[-1] = 360.1
sector = np.digitize(np.mod(wd + width / 2.0, 360.0), edges) - 1
sector = sector.astype(np.float64)
sector[sector >= n_sectors] = np.nan
return sector
def _wd_to_sector_quantiles(wd, n_sectors=12):
# TODO move this to xarray nor numpy so we can use apply_ufunc
sector_da = wd.copy()
sector_cat, edges = pd.qcut(wd.values.flatten(), n_sectors, retbins=True)
edges[0] = 0.0
edges[-1] = 360.0
sector_da.data = sector_cat.codes.reshape(wd.shape)
sector_coords_da = create_sector_coords(edges)
return sector_da, sector_coords_da
if output_type not in ["centers", "indices"]:
raise ValueError("unkown output type. Possible values are 'centers','indices'")
if not quantiles:
sector_indices = xr.apply_ufunc(
_wd_to_sector_constant, wd, kwargs={"n_sectors": sectors}
)
sector_coords = create_sector_coords(sectors)
sector_centers = sector_indices * 360.0 / sectors
else:
if (
type(wd) is not xr.DataArray
or ("point" not in wd.dims)
or (len(wd["point"]) > 1)
):
raise ValueError(
"For quantiles=True, only xarray.DataArray with point dimensions of length 1 are supported"
)
sector_indices, sector_coords = _wd_to_sector_quantiles(wd, sectors)
centers_values = sector_coords.isel(
sector=sector_indices.values.flatten()
).values
sector_centers = wd.copy()
sector_centers.values = centers_values.reshape(-1, 1)
if output_type == "indices":
return sector_indices, sector_coords
else:
return sector_centers, sector_coords
[docs]
def vinterp_wind_direction(wind_direction, height, **kwargs):
"""
Interpolate wind direction to a given height.
Parameters
----------
wind_direction : xarray.DataArray
Wind direction.
height : float
Height to interpolate wind direction to.
**kwargs : dict, optional
Additional keyword arguments passed to xarray.interp.
Returns
-------
wind_direction : xarray.DataArray
Interpolated wind direction.
"""
if not isinstance(wind_direction, xr.DataArray):
raise TypeError("wind_direction must be a xarray.DataArray")
if "height" not in wind_direction.dims:
raise ValueError("wind_direction must have a height dimension")
if not isinstance(height, (np.ScalarType, xr.DataArray)):
raise TypeError("height must be a scalar or xarray.DataArray")
wd_ref = wind_direction.isel(height=0)
wd_diff = wind_direction_difference(wind_direction, wd_ref)
wd_new = wd_ref - wd_diff.interp(height=height, **kwargs)
return np.mod(wd_new, 360.0)
[docs]
def vinterp_wind_speed(wind_speed, height, log_height=True, **kwargs):
"""
Vertically interpolate wind speed to a given height from other height levels.
Parameters
----------
wind_speed : xarray.DataArray
Wind speed. Must have a height dimension.
height : float, xarray.DataArray
Height to interpolate wind speed to.
log_height : bool, optional
If True, interpolate in log-height space. Defaults to True.
**kwargs : dict, optional
Additional keyword arguments passed to xarray.interp.
Returns
-------
wind_speed : xarray.DataArray
Interpolated wind speed.
"""
if not isinstance(wind_speed, xr.DataArray):
raise TypeError("wind_speed must be a xarray.DataArray")
if "height" not in wind_speed.dims:
raise ValueError("wind_speed must have a height dimension")
if not isinstance(height, (np.ScalarType, xr.DataArray)):
raise TypeError("height must be a scalar or xarray.DataArray")
wind_speed = wind_speed.copy()
if log_height:
wind_speed = wind_speed.assign_coords(height=np.log1p(wind_speed.height))
if isinstance(height, xr.DataArray):
height_ = height.copy()
height = np.log1p(height)
wind_speed = wind_speed.interp(height=height, **kwargs)
if log_height and isinstance(height_, xr.DataArray):
wind_speed = wind_speed.assign_coords(height=height_)
return wind_speed
[docs]
def rotor_equivalent_wind_speed(
wind_speed,
wind_direction,
hub_height,
rotor_diameter,
delta_z=1.0,
n_integrate=1001,
):
"""
Calculate the rotor equivalent wind speed (REWS) from given wind speed and directions
on height levels.
The procedure is as follows:
1. Find the area of each segment of the rotor spanned area.
2. Calculate the wind speed at the center of each segment by linearly interpolating
the wind speed to the height of the segment center in log-height.
3. Calculate the wind direction at the center of each segment by linearly interpolating
the wind direction to the height of the segment center. Circularity is
taken into account here.
4. Calculate the wind direction at hub height by linearly interpolating the wind
direction to the hub height.
5. Calculate the REWS as the cube root of the sum of the wind speed at each segment
center multiplied by the area-weight (area/total) of the segment and the cosine
of the difference between the wind direction at the segment center and the wind
direction at hub height.
Parameters
----------
wind_speed : xarray.DataArray
Wind speed on height levels.
wind_direction : xarray.DataArray
Wind direction on height levels.
hub_height : float
Turbine Hub height.
rotor_diameter : float
Turbine rotor diameter.
delta_z : float, optional
Height difference between segments of turbine spanned rotor area
(default: 1.0).
n_integrate : int, optional
Number of points to use for integration (default: 1001) of the area
of each segment.
Returns
-------
rews : xarray.DataArray
Rotor equivalent wind speed.
"""
if not isinstance(wind_speed, xr.DataArray):
raise TypeError("wind_speed must be a xarray.DataArray")
if not isinstance(wind_direction, xr.DataArray):
raise TypeError("wind_direction must be a xarray.DataArray")
if "height" not in wind_speed.dims:
raise ValueError("wind_speed must have a height dimension")
if "height" not in wind_direction.dims:
raise ValueError("wind_direction must have a height dimension")
hub_height = float(hub_height)
rotor_diameter = float(rotor_diameter)
delta_z = float(delta_z)
n_integrate = int(n_integrate)
rotor_radius = rotor_diameter / 2.0
zi = np.linspace(
hub_height - rotor_radius,
hub_height + rotor_radius,
int(np.round(rotor_diameter / delta_z)) + 1,
)
zc = (zi[1:] + zi[:-1]) / 2
zc = xr.DataArray(zc, dims=("height",), coords={"height": zc})
Ai = np.zeros_like(zc)
for i in range(len(zi) - 1):
zs = np.linspace(zi[i], zi[i + 1], n_integrate)
Ai[i] = trapezoid(2 * np.sqrt(rotor_radius**2 - (zs - hub_height) ** 2), zs)
# Area of rotor
A = np.pi * rotor_radius**2
# Area of rotor segment
Ai = xr.DataArray(Ai, dims=("height",), coords={"height": zc})
# wind speed and direction at segment center
ui = vinterp_wind_speed(
wind_speed, zc, method="linear", kwargs={"fill_value": "extrapolate"}
)
di = vinterp_wind_direction(
wind_direction, zc, method="linear", kwargs={"fill_value": "extrapolate"}
)
# wind direction at hub height
dh = vinterp_wind_direction(
wind_direction,
hub_height,
method="linear",
kwargs={"fill_value": "extrapolate"},
)
rews = ((1 / A) * Ai * (ui**3) * np.cos(np.deg2rad(di - dh))).sum(dim="height")
rews = xr.where(rews < 0, 0, rews)
rews = np.power(rews, 1.0 / 3.0)
rews = rews.expand_dims(height=[hub_height])
return rews
[docs]
def tswc_resample(
ds,
freq,
var_ws="wind_speed",
var_wd="wind_direction",
min_availability=0.5,
**kwargs,
):
"""Resample wind speed and direction to a given frequency.
Parameters
----------
ds : xarray.Dataset
Dataset with wind speed and direction.
freq : str
Resampling frequency.
var_ws : str, optional
Name of wind speed variable, by default "wind_speed".
var_wd : str, optional
Name of wind direction variable, by default "wind_direction".
Returns
-------
xarray.Dataset
Resampled dataset.
"""
ds = ds.copy()
def nan_mean(da):
return da.mean(dim="time").where(
da.notnull().sum(dim="time") >= len(da.time) * min_availability
)
ds["__U__"], ds["__V__"] = wind_vectors(ds[var_ws], ds[var_wd])
ws_resampled = ds[var_ws].resample(time=freq, **kwargs).map(nan_mean)
ds = ds.drop_vars([var_ws, var_wd])
ds = ds.resample(time=freq, **kwargs).map(nan_mean)
_, ds[var_wd] = wind_speed_and_direction(ds["__U__"], ds["__V__"])
ds[var_ws] = ws_resampled
ds = ds.drop_vars(["__U__", "__V__"])
return ds
[docs]
def shear_exponent(da):
"""
Compute the shear exponent from vertical wind speed profiles using finite
differences in log-space.
Parameters
----------
da : xarray.DataArray
Wind speed DataArray. The DataArray must contain
a height coordinate and horizontal coordinates ('west_east', 'south_north')
or be convertible to that structure. Wind speeds must be positive.
Returns
-------
xarray.DataArray
The shear exponent computed as (d ln u) / (d ln z) at log-space midpoints between input heights.
The returned DataArray has the same spatial structure as the input
and only contains the points where valid shear calculations could be made.
Raises
------
ValueError
If the input cannot be converted to a point/vertical structure by
windkit.spatial.to_point.
Notes
-----
- The function sorts profiles by (west_east, south_north, height) and
computes backward finite differences along the "point" axis.
- Midpoint heights are computed in log-space and assigned to the output
shear DataArray.
Examples
--------
>>> shear = shear_exponent(da)
"""
# convert to point structure
da = _spatial_stack(da)
# sort by x,y,height so that we can calculate forward differences between all points
da = da.sortby(["west_east", "south_north", "height"])
ln_u = np.log(da).diff("point", label="lower")
ln_z = (np.log(da["height"])).diff("point", label="lower")
shear = ln_u / ln_z
shear.name = "shear_exponent"
# find all points that are belonging to the same vertical profile with same x,y
same_xy = (da["west_east"].diff("point") == 0) & (
da["south_north"].diff("point") == 0
)
shear = shear.where(same_xy)
# calculate heights at mid-points in log space
logz = np.log(da["height"])[:-1] + 0.5 * ln_z
shear.coords["height"] = np.exp(logz)
shear = shear.assign_coords(same_xy["point"].drop_vars("height").coords)
# this is a bug in xarray where it did not keep attributes, only required for xarray<2026.01
# can be removed when support for older xarray versions is dropped
shear = shear.assign_attrs(da.attrs)
# drop invalid height that results from non-matching vertical profiles
shear = shear.where(same_xy.drop_vars("height"), drop=True)
# Remove the _pwio attributes as they are no longer valid for the reduced dimension
if "_pwio_stacked_index_names" in shear.attrs:
del shear.attrs["_pwio_stacked_index_names"]
# convert back to original structure
shear = _spatial_unstack(shear)
return shear
[docs]
def wind_veer(da):
"""
Calculate wind veer (change in wind direction with height).
Parameters
----------
da : xarray.DataArray
DataArray containing wind direction data. Must have coordinates for
'west_east', 'south_north', and 'height'.
Returns
-------
xarray.DataArray
DataArray containing the wind veer values in degrees per meter.
The structure matches the input structure (point, stacked_point,
raster, or cuboid).
"""
# convert to point structure
da = _spatial_stack(da)
# sort by x,y,height so that we can calculate forward differences between all points
da = da.sortby(["west_east", "south_north", "height"])
wd_diff = da.diff("point", label="lower")
# Adjust for circular nature of wind direction
wd_diff = (wd_diff + 180) % 360 - 180
z_diff = da["height"].diff("point", label="lower")
veer = wd_diff / z_diff
veer.name = "wind_veer"
# find all points that are belonging to the same vertical profile with same x,y
same_xy = (da["west_east"].diff("point") == 0) & (
da["south_north"].diff("point") == 0
)
veer = veer.where(same_xy)
# calculate heights at mid-points
new_z = da["height"][:-1] + 0.5 * z_diff
veer.coords["height"] = new_z
veer = veer.assign_coords(same_xy["point"].drop_vars("height").coords)
# this is a bug in xarray where it did not keep attributes, only required for xarray<2026.01
# can be removed when support for older xarray versions is dropped
veer = veer.assign_attrs(da.attrs)
# drop invalid height that results from non-matching vertical profiles
veer = veer.where(same_xy.drop_vars("height"), drop=True)
# Remove the _pwio attributes as they are no longer valid for the reduced dimension
if "_pwio_stacked_index_names" in veer.attrs:
del veer.attrs["_pwio_stacked_index_names"]
# convert back to original structure
veer = _spatial_unstack(veer)
return veer
