# (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 xarray as xr
from scipy.integrate import trapezoid
from windkit._errors import WindkitValidationError
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. The sector edges are the quantiles of all values in ``wd``
pooled together, so an input with several points gets one shared set of
edges rather than one set per point. Unlike the fixed width layout, whose
sector 0 is centred on north, quantile sector 0 starts at north. 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.
Raises
------
ValueError
If ``output_type`` is unknown or ``sectors`` is not positive.
WindkitValidationError
If ``quantiles=True`` and ``wd`` is an :py:class:`xarray.Dataset`, has
no finite values, or holds too few distinct directions to build
``sectors`` sectors with unique quantile edges.
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 _quantile_edges(wd_values, n_sectors):
valid = wd_values[np.isfinite(wd_values)]
n_distinct = np.unique(valid).size
if n_distinct == 0:
raise WindkitValidationError(
f"Cannot build {n_sectors} equal-probability sectors from "
"0 distinct wind directions: no finite values in wd."
)
# Quantiles interpolate between order statistics, so fewer distinct
# values than sectors can still give unique edges. Only the edges
# themselves decide.
edges = np.quantile(valid, np.linspace(0.0, 1.0, n_sectors + 1))
if np.unique(edges).size != n_sectors + 1:
raise WindkitValidationError(
f"Cannot build {n_sectors} equal-probability sectors from "
f"{n_distinct} distinct wind directions: the quantile "
"edges are not unique."
)
edges[0] = 0.0
edges[-1] = 360.0
return edges
def _wd_to_sector_quantiles(wd, edges):
# Left-closed bins on the interior edges, like the constant width layout.
sector = np.digitize(wd, edges[1:-1]).astype(np.float64)
sector[~np.isfinite(wd)] = np.nan
return sector
def _indices_to_centers(sector, centers):
out = np.full(sector.shape, np.nan)
valid = ~np.isnan(sector)
out[valid] = centers[sector[valid].astype(int)]
return out
if output_type not in ["centers", "indices"]:
raise ValueError("unkown output type. Possible values are 'centers','indices'")
if sectors < 1:
raise ValueError(f"'sectors' must be a positive integer, got {sectors}")
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 isinstance(wd, xr.Dataset):
raise WindkitValidationError(
"For quantiles=True, wd must be an xarray.DataArray or a numpy "
"array, not an xarray.Dataset"
)
edges = _quantile_edges(np.asarray(wd, dtype=np.float64).ravel(), sectors)
sector_indices = xr.apply_ufunc(
_wd_to_sector_quantiles, wd, kwargs={"edges": edges}
)
sector_coords = create_sector_coords(edges)
sector_centers = xr.apply_ufunc(
_indices_to_centers,
sector_indices,
kwargs={"centers": sector_coords.values},
)
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
_NO_PROFILE = (
"{name} differences along a vertical profile, so it needs at least two "
"heights at the same horizontal location. {detail}"
)
def _require_vertical_profile(da, name):
"""Reject an input with no vertical profile, before it is stacked.
Parameters
----------
da : xarray.DataArray
Input, in any spatial structure.
name : str
Name of the calling function, for the error message.
Raises
------
WindkitValidationError
If ``da`` has no ``height`` coordinate, or fewer than two heights.
"""
if "height" not in da.coords:
raise WindkitValidationError(
_NO_PROFILE.format(name=name, detail="The input has no height coordinate.")
)
if da["height"].size < 2:
raise WindkitValidationError(
_NO_PROFILE.format(
name=name, detail=f"The input carries {da['height'].size} height(s)."
)
)
def _usable_step(da, name):
"""Flag consecutive stacked points forming a usable profile step.
A usable step needs the same horizontal position and two distinct
heights, both compared for exact equality against ``da``'s own
coordinate values, with no tolerance applied. A NaN height counts as
distinct from any other value, including another NaN (``NaN != x`` is
always true), so a step with one is not filtered out here; the
resulting NaN difference propagates into the output silently rather
than raising, the same way ``same_xy`` already treated NaN positions
before this guard existed.
Parameters
----------
da : xarray.DataArray
Stacked and sorted input, on a flat ``point`` axis.
name : str
Name of the calling function, for the error message.
Returns
-------
xarray.DataArray
Boolean flag over ``point``, one shorter than ``da``'s point axis.
Raises
------
WindkitValidationError
If no horizontal location carries two distinct heights.
"""
same_xy = (da["west_east"].diff("point") == 0) & (
da["south_north"].diff("point") == 0
)
distinct_height = da["height"].diff("point") != 0
usable = same_xy & distinct_height
if not bool(usable.any()):
raise WindkitValidationError(
_NO_PROFILE.format(
name=name,
detail=(
f"The input has {da.sizes['point']} points, no two of them "
"at the same west_east/south_north with distinct heights."
),
)
)
return usable
[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
------
WindkitValidationError
If the input carries fewer than two distinct heights at the same
horizontal location, leaving no vertical profile to difference along.
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)
"""
_require_vertical_profile(da, "shear_exponent")
# 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"])
# find all consecutive point pairs forming a usable vertical profile step:
# same x,y and two distinct heights
usable = _usable_step(da, "shear_exponent")
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"
shear = shear.where(usable)
# 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(usable["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(usable.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).
Raises
------
WindkitValidationError
If the input carries fewer than two distinct heights at the same
horizontal location, leaving no vertical profile to difference along.
"""
_require_vertical_profile(da, "wind_veer")
# 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"])
# find all consecutive point pairs forming a usable vertical profile step:
# same x,y and two distinct heights
usable = _usable_step(da, "wind_veer")
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"
veer = veer.where(usable)
# calculate heights at mid-points
new_z = da["height"][:-1] + 0.5 * z_diff
veer.coords["height"] = new_z
veer = veer.assign_coords(usable["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(usable.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
