# (c) 2026 DTU Wind Energy
"""Time series wind climate module
A time series wind climate is defined by dataset with a time series
``wind speed`` and ``wind direction``.
A valid time series wind climate therefore has a dimension ``time``.
Also it must have one of the valid :ref:`geospatial_structures`. This module contains
functions that operate on time series wind climates.
This includes the ability to create time series datasets from files and from
existing data.
"""
__all__ = [
"validate_tswc",
"is_tswc",
"create_tswc",
"read_tswc",
"tswc_to_file",
"tswc_from_dataframe",
"create_tswc_pair",
]
import collections
import logging
import re
import warnings
import numpy as np
import pandas as pd
import xarray as xr
from windkit.xarray_structures._schema import (
DataVar,
DatasetSchema,
make_is_valid,
make_validator,
validate,
validate_input,
)
from windkit.xarray_structures.empty import (
_copy_chunks,
_native_batch,
_native_batch_horizontal,
_native_flat_coords,
_native_frame,
_native_store,
)
from windkit.xarray_structures.metadata import (
_ALL_VARS_META,
_TS_ATTRS,
_update_history,
_update_var_attrs,
)
from windkit.xarray_structures.sector import create_sector_coords
from windkit.xarray_structures.wsbin import create_wsbin_coords
from windkit.spatial import to_stacked_point
from windkit.spatial._crs import set_crs
from windkit.utils import _infer_file_format
from windkit.wind import shear_extrapolate, wd_to_sector
from windkit.wind_climate._correlated_fields import (
_column_index,
_column_of,
_correlated_fields,
_seconds_since_start,
_stability_field_ranges,
)
logger = logging.getLogger(__name__)
SUPPORTED_TSWC_FILE_FORMATS_READ = ["csv", "txt", "nc", "yml", "yaml"]
SUPPORTED_TSWC_FILE_FORMATS_WRITE = ["yml", "yaml"]
WS = "wind_speed"
WD = "wind_direction"
DIM_TIME = "time"
DIM_SECTOR = "sector"
DIM_WSBIN = "wsbin"
_DEFAULT_DATE_RANGE = pd.date_range(
"2001-01-01", "2002-01-01", freq="10min", inclusive="left"
)
# Windpro date format tokens → Python strftime codes, pre-sorted longest-first
# to avoid partial matches during sequential replacement (e.g. "MM" before "M").
# Note: windpro uses "hh" for 24-hour clock (unlike Java/.NET where hh = 12-hour),
# confirmed by real windpro exports containing hours up to 23.
_WINDPRO_TO_STRFTIME = (
("yyyy", "%Y"),
("MM", "%m"),
("dd", "%d"),
("HH", "%H"),
("hh", "%H"),
("mm", "%M"),
("ss", "%S"),
("yy", "%y"),
)
_TSWC_SCHEMA = DatasetSchema(
name="tswc",
data_vars={WS: DataVar(dims=(DIM_TIME,)), WD: DataVar(dims=(DIM_TIME,))},
dims=(DIM_TIME,),
coords=("time",),
)
validate_tswc = make_validator(_TSWC_SCHEMA)
is_tswc = make_is_valid(_TSWC_SCHEMA)
[docs]
def create_tswc(
output_locs,
date_range=None,
not_empty=True,
seed=9876538,
weibull_A=8.0,
weibull_k=2.1,
speed_tau=14400.0,
direction_mean=240.0,
direction_kappa=1.5,
direction_tau=21600.0,
spatial_correlation_length=10_000.0,
fields=None,
):
"""Create time series wind climate dataset.
When ``not_empty=True`` (default), the data variables are filled with
temporally and spatially correlated synthetic wind data using an AR(1)
process in Gaussian space transformed to a Weibull distribution for speed,
and a circular AR(1) with von Mises innovations for direction.
Parameters
----------
output_locs : xarray.Dataset
Output geospatial information.
date_range : pandas.DatetimeIndex or None
Time range as a pandas DateTimeIndex. If None is passed, a default
range of one year (2001) at 10 min resolution is created.
not_empty : bool
If True, the empty dataset is filled with realistic synthetic
wind data. Defaults to True.
seed : int
Seed for the random data, defaults to 9876538.
weibull_A : float
Weibull scale parameter [m/s]. Defaults to 8.0.
weibull_k : float
Weibull shape parameter. Defaults to 2.1.
speed_tau : float
Autocorrelation e-folding timescale for speed [s]. Defaults to 14400 (4 h).
direction_mean : float
Prevailing wind direction [degrees]. Defaults to 240.
direction_kappa : float
Von Mises concentration around prevailing direction. Defaults to 1.5.
direction_tau : float
Autocorrelation timescale for direction [s]. Defaults to 21600 (6 h).
spatial_correlation_length : float
Spatial e-folding length scale for inter-point correlation [same unit
as coordinates]. Defaults to 10000.
fields : str, sequence of str or None
Stability variables to add alongside the wind: ``"stability"`` or
``"all"`` for both ``temp_scale`` (:math:`T_*`) and
``ustar_over_pblh`` (:math:`u_*/z_i`), or a subset of those names.
Defaults to None, which adds neither. The values are bounded by the
0.1st to 99.9th percentile of an onshore mesoscale reference and carry
that reference's correlations with each other and with the generated
wind speed. Both are surface-layer scales, declared
``height_dim_policy: "forbidden"``, so they are drawn once per
horizontal column - driven by the wind at its lowest height - and
carry no height dimension, while the wind keeps its own.
``temp_scale`` follows the convention in which
:math:`T_*` has the sign of the surface heat flux, so positive is
unstable and negative stable
(https://doi.org/10.1007/s10546-023-00803-3).
Returns
-------
ds : xarray.Dataset
Time series wind climate dataset.
Raises
------
ValueError
If ``fields`` names something other than a stability variable, or if
the time range spans 2 days or less, which is too short to generate
correlated stability variables.
Notes
-----
``wind_speed`` and ``wind_direction`` are stored float32, C-contiguous,
with ``time`` trailing the full spatial batch; the surface-layer scales
are stored the same way over the horizontal batch. That is the reverse of
the compiled core's declared argument order, so PyWAsP's compiled boundary
can take a Fortran-ordered view instead of converting a copy
(windkit#807). Coordinates stay float64: the core declares the location
arguments of ``predict_*_points`` double precision, and float32 spacing at
5e5 m is 0.031 m. Neither the dtype nor the dimension order is a supported
API contract: ``create_tswc`` is a test-data utility.
Examples
--------
Create a single-point TSWC with default parameters:
>>> import pandas as pd
>>> from windkit.spatial import create_dataset
>>> import windkit as wk
>>> locs = create_dataset(10.0, 55.0, 100.0, 4326).drop_vars("output")
>>> dates = pd.date_range("2020-01-01", periods=8760, freq="h")
>>> tswc = wk.create_tswc(locs, date_range=dates)
>>> tswc
<xarray.Dataset> ...
Create a multi-point TSWC with custom Weibull parameters:
>>> import numpy as np
>>> locs = create_dataset(
... [10.0, 10.1], [55.0, 55.1], [100.0, 100.0],
... 4326, struct="point",
... ).drop_vars("output")
>>> tswc = wk.create_tswc(locs, date_range=dates, weibull_A=9.5, weibull_k=2.3)
Create an empty TSWC (NaN-filled):
>>> tswc_empty = wk.create_tswc(locs, date_range=dates, not_empty=False)
>>> np.all(np.isnan(tswc_empty["wind_speed"].values))
True
Create a TSWC carrying the stability variables:
>>> tswc = wk.create_tswc(locs, date_range=dates, fields="stability")
>>> sorted(tswc.data_vars)
['temp_scale', 'ustar_over_pblh', 'wind_direction', 'wind_speed']
"""
return _create_tswc_core(
output_locs,
date_range=date_range,
not_empty=not_empty,
seed=seed,
weibull_A=weibull_A,
weibull_k=weibull_k,
speed_tau=speed_tau,
direction_mean=direction_mean,
direction_kappa=direction_kappa,
direction_tau=direction_tau,
spatial_correlation_length=spatial_correlation_length,
fields=fields,
_return_internals=False,
)
def _create_tswc_core(
output_locs,
date_range=None,
not_empty=True,
seed=9876538,
weibull_A=8.0,
weibull_k=2.1,
speed_tau=14400.0,
direction_mean=240.0,
direction_kappa=1.5,
direction_tau=21600.0,
spatial_correlation_length=10_000.0,
fields=None,
_return_internals=False,
):
"""Internal helper for create_tswc with optional access to its internals.
With ``_return_internals``, also returns the pre-storage arrays the two
in-house callers need: the latent speed and direction fields
:func:`create_tswc_pair` derives its second series from, and the float64
wind speed :func:`create_surface_fields` drives its fields with. All three
are ``(n_timesteps, n_locations)``, laid out over the flattened batch
:func:`_native_batch` gives, and none of them can be recovered from the
returned dataset, whose variables are stored float32 in a different
dimension order (windkit#807).
"""
from windkit.wind_climate._synthetic_time_series import (
_generate_latent_direction,
_generate_latent_speed,
_latent_to_weibull,
_spatial_correlation_cholesky,
)
# Resolved before anything is built, so a misspelled field fails fast.
field_ranges = _stability_field_ranges(fields)
if date_range is None:
time_values = _DEFAULT_DATE_RANGE
elif type(date_range) is not pd.DatetimeIndex:
raise TypeError("date_range must be of type pandas.DatetimeIndex")
else:
time_values = date_range
# The wind carries height, the surface-layer scales do not, so the two
# groups are built over different batches. Both are laid out over
# `output_locs`'s own dimensions, so a scalar height or a reduced
# horizontal axis stays reduced.
dims, shape = _native_batch(output_locs)
horiz_dims, horiz_shape = _native_batch_horizontal(output_locs)
n_pt = int(np.prod(shape, dtype=int))
n_horiz = int(np.prod(horiz_shape, dtype=int))
n_timesteps = len(time_values)
z_latent = None
theta_latent = None
generated = {}
if not_empty:
rng = np.random.default_rng(seed)
we, sn = _native_flat_coords(output_locs, dims, shape)
chol_L = _spatial_correlation_cholesky(we, sn, spatial_correlation_length)
dt = np.median(np.diff(time_values)).astype("timedelta64[s]").astype(float)
z_latent = _generate_latent_speed(n_timesteps, n_pt, dt, speed_tau, chol_L, rng)
ws = _latent_to_weibull(z_latent, weibull_A, weibull_k)
theta_latent = _generate_latent_direction(
n_timesteps,
n_pt,
dt,
direction_mean,
direction_kappa,
direction_tau,
chol_L,
rng,
)
wd = np.rad2deg(theta_latent) % 360.0
if field_ranges:
# Surface-layer scales have no vertical variation, so they are
# drawn once per horizontal column - driven by the wind at its
# lowest height - and then spread over the horizontal batch.
# Where height is its own dimension that spread is a permutation;
# in a `point` structure, where height is folded into the flat
# axis, it repeats each column's value over the heights above it.
height = _native_flat_coords(output_locs, dims, shape, ("height",))[0]
_, base_of_column = _column_index(we, sn, height)
column_of_location = _column_of(
*_native_flat_coords(output_locs, horiz_dims, horiz_shape)
)
stability = _correlated_fields(
field_ranges,
_seconds_since_start(pd.DatetimeIndex(time_values)),
len(base_of_column),
np.random.default_rng(seed + 1),
ws[:, base_of_column],
)
generated = {
name: values[:, column_of_location]
for name, values in stability.items()
}
else:
ws = np.full((n_timesteps, n_pt), np.nan)
wd = np.full((n_timesteps, n_pt), np.nan)
generated = {
name: np.full((n_timesteps, n_horiz), np.nan) for name in field_ranges
}
ds = _native_frame(output_locs).assign_coords(time=time_values)
for name, values in (("wind_speed", ws), ("wind_direction", wd)):
out_dims, array = _native_store(values, dims, shape, core_dims=("time",))
ds[name] = xr.DataArray(array, dims=out_dims)
for name, values in generated.items():
out_dims, array = _native_store(
values, horiz_dims, horiz_shape, core_dims=("time",)
)
ds[name] = xr.DataArray(array, dims=out_dims)
stability_meta = {k: _ALL_VARS_META[k] for k in field_ranges}
ds = _update_var_attrs(
_copy_chunks(output_locs, ds),
{**_TS_ATTRS, **stability_meta},
)
ds = _update_history(ds)
if _return_internals:
return ds, {
"z_latent": z_latent,
"theta_latent": theta_latent,
"wind_speed": ws,
}
return ds
def _pack_tswc_dataset(output_locs, date_range, speed_array, direction_array):
"""Build a TSWC dataset from pre-computed speed and direction arrays.
Parameters
----------
output_locs : xarray.Dataset
Output geospatial information.
date_range : pandas.DatetimeIndex
Time coordinate values.
speed_array : numpy.ndarray
Wind speed array, shape ``(n_timesteps, n_points)``.
direction_array : numpy.ndarray
Wind direction array, shape ``(n_timesteps, n_points)``.
Returns
-------
ds : xarray.Dataset
Time series wind climate dataset.
"""
dims, shape = _native_batch(output_locs)
ds = _native_frame(output_locs).assign_coords(time=date_range)
for name, values in (
("wind_speed", speed_array),
("wind_direction", direction_array),
):
out_dims, array = _native_store(values, dims, shape, core_dims=("time",))
ds[name] = xr.DataArray(array, dims=out_dims)
ds = _update_var_attrs(_copy_chunks(output_locs, ds), _TS_ATTRS)
return _update_history(ds)
[docs]
def create_tswc_pair(
output_locs,
date_range=None,
seed=9876538,
weibull_A=(8.0, 7.2),
weibull_k=(2.1, 2.3),
target_r2=0.80,
smooth_sigma=1.5,
speed_tau=14400.0,
direction_mean=240.0,
direction_kappa=1.5,
direction_tau=21600.0,
direction_bias=5.0,
direction_noise_std=15.0,
spatial_correlation_length=10_000.0,
):
"""Create a pair of correlated synthetic TSWC datasets.
Generates a first TSWC with realistic temporal and spatial correlations,
then derives a second TSWC that is correlated with the first but has
different Weibull parameters (speed bias) and a directional bias.
Parameters
----------
output_locs : xarray.Dataset
Output geospatial information.
date_range : pandas.DatetimeIndex or None
Time range. If None, defaults to one year (2001) at 10 min resolution.
seed : int
Random seed for reproducibility.
weibull_A : tuple of float
Weibull scale parameters [m/s] for the first and second dataset.
weibull_k : tuple of float
Weibull shape parameters for the first and second dataset.
target_r2 : float
Target Pearson r-squared between the two datasets' speeds.
smooth_sigma : float
Gaussian filter width in time steps for latent smoothing.
speed_tau : float
Autocorrelation e-folding timescale for speed [s].
direction_mean : float
Prevailing wind direction [degrees].
direction_kappa : float
Von Mises concentration around prevailing direction.
direction_tau : float
Autocorrelation timescale for direction [s].
direction_bias : float
Constant directional bias of the second dataset relative to the
first [degrees].
direction_noise_std : float
Std of Gaussian angular noise added to the second dataset's
direction [degrees].
spatial_correlation_length : float
Spatial e-folding length scale for inter-point correlation within
each dataset. Only relevant when ``output_locs`` contains multiple
points.
Returns
-------
tswc_1 : xarray.Dataset
First TSWC dataset.
tswc_2 : xarray.Dataset
Second TSWC dataset.
Notes
-----
Both datasets are stored the way :func:`create_tswc` stores its own
(windkit#807).
Examples
--------
Create a measurement/reference pair for MCP testing:
>>> import pandas as pd
>>> from windkit.spatial import create_dataset
>>> import windkit as wk
>>> locs = create_dataset(10.0, 55.0, 100.0, 4326).drop_vars("output")
>>> dates = pd.date_range("2020-01-01", periods=8760, freq="h")
>>> meas, ref = wk.create_tswc_pair(locs, date_range=dates)
Customize Weibull parameters and target correlation:
>>> meas, ref = wk.create_tswc_pair(
... locs,
... date_range=dates,
... weibull_A=(9.0, 7.5),
... weibull_k=(2.0, 2.4),
... target_r2=0.85,
... )
Both outputs are valid TSWC datasets:
>>> wk.validate_tswc(meas)
>>> wk.validate_tswc(ref)
"""
from windkit.wind_climate._synthetic_time_series import _generate_reference
A_1, A_2 = weibull_A
k_1, k_2 = weibull_k
ds_1, internals = _create_tswc_core(
output_locs,
date_range=date_range,
not_empty=True,
seed=seed,
weibull_A=A_1,
weibull_k=k_1,
speed_tau=speed_tau,
direction_mean=direction_mean,
direction_kappa=direction_kappa,
direction_tau=direction_tau,
spatial_correlation_length=spatial_correlation_length,
_return_internals=True,
)
rng = np.random.default_rng(seed + 1)
if date_range is None:
time_values = _DEFAULT_DATE_RANGE
else:
time_values = date_range
dt = np.median(np.diff(time_values)).astype("timedelta64[s]").astype(float)
speed_2, direction_2 = _generate_reference(
internals["z_latent"],
internals["theta_latent"],
ref_A=A_2,
ref_k=k_2,
target_r2=target_r2,
smooth_sigma=smooth_sigma,
ref_direction_bias=direction_bias,
ref_direction_noise_std=direction_noise_std,
dt=dt,
speed_tau=speed_tau,
rng=rng,
)
if date_range is None:
date_range = _DEFAULT_DATE_RANGE
ds_2 = _pack_tswc_dataset(output_locs, date_range, speed_2, direction_2)
return ds_1, ds_2
def _read_ts_windpro_txt(fpath):
"""Parses windpro format txt file into a dataset.
Parameters
----------
fpath : [str]
[file path to be parsed]
Returns
-------
xarray.Dataset
"""
def _is_float(value):
try:
float(value)
return True
except ValueError:
return False
# parse parameters from windpro header
header_line = None
date_format = None
lng, lat = None, None
disp_height = -1
with open(fpath, "r") as file:
for i, line in enumerate(file):
if "Geographical Coordinates" in line:
parts = line.split()
for j, p in enumerate(parts):
if _is_float(p) and parts[j - 1] == "Longitude:":
lng = float(p)
if _is_float(p) and parts[j - 1] == "Latitude:":
lat = float(p)
if "Displacement height" in line:
parts = line.split()
for p in parts:
if _is_float(p):
disp_height = float(p)
if "Date time format" in line:
fmt_str = line.split(":", 1)[-1].strip()
date_format = fmt_str
for token, code in _WINDPRO_TO_STRFTIME:
date_format = date_format.replace(token, code)
if "TimeStamp" in line:
header_line = i
break
if disp_height > 0:
warnings.warn(
"Displacement height cannot be used in WindKit. Set it up via the map instead."
)
if lng is None or lat is None:
raise ValueError("Couldn't parse coordinates from windpro file header")
if header_line is None:
raise ValueError("Could not find 'TimeStamp' header in windpro file")
if date_format is None:
date_format = "%d-%m-%Y %H:%M"
warnings.warn(
"Could not find 'Date time format' in windpro header, "
f"falling back to '{date_format}'."
)
# skip_blank_lines=False ensures header_line index matches the actual line
# number in the file (pandas would shift it otherwise by skipping blank lines).
# skiprows skips the units row immediately after the header.
ts = pd.read_csv(
fpath,
delimiter="\t{2}|\t",
skip_blank_lines=False,
header=header_line,
skiprows=[header_line + 1],
engine="python",
)
try:
ts["TimeStamp"] = pd.to_datetime(ts["TimeStamp"], format=date_format)
except ValueError as e:
raise ValueError(
f"Could not parse timestamps with format '{date_format}'. "
f"The windpro file may use a different date format. Original error: {e}"
) from e
# parse height from the wind speed/direction column
height = None
for col in ts.columns:
if "Mean wind speed" in col:
matches = re.findall(r"[0-9]+(?:\.[0-9]+)?m", col)
if matches:
height = float(matches[0].replace("m", ""))
ts = ts.rename({col: "ws"}, axis="columns")
if "Wind direction" in col:
ts = ts.rename({col: "wd"}, axis="columns")
if "ws" not in ts.columns or "wd" not in ts.columns:
raise ValueError(
"Could not find 'Mean wind speed' and/or 'Wind direction' columns"
)
if height is None:
raise ValueError("Could not parse measurement height from column names")
ts = ts[~ts.ws.str.contains("-")]
ts = ts[ts["ws"].notna()]
ts = ts[ts["wd"].notna()]
ts["ws"] = ts["ws"].astype(float)
ts["wd"] = ts["wd"].astype(float)
ts_ds = xr.Dataset(
{
"wind_speed": (["time"], ts["ws"]),
"wind_direction": (["time"], ts["wd"]),
},
coords={
"time": ("time", ts["TimeStamp"]),
"south_north": lat,
"west_east": lng,
"height": height,
"crs": 0,
},
)
set_crs(ts_ds, 4326)
_update_var_attrs(ts_ds, {**_TS_ATTRS})
# validate the type schema only (scalar coordinates, no spatial dims yet)
validate(ts_ds, _TSWC_SCHEMA)
return ts_ds
def _read_csv(
file_name,
west_east=0.0,
south_north=0.0,
crs=4326,
time_col=0,
height_to_columns=None,
pandas_kwargs=None,
):
"""
Reads a csv file into a time series wind climate xarray.Dataset. The file must have one time
entry per row, a column with a time stamp and at least one wind speed and one wind direction. It
allows to create a dataset for several heights.
Parameters
----------
file_name : str
file path to a csv file with wind speed and wind direction measurements for different timestamps.
west_east: float, optional
west east locaton of the measurement. Defaults to 0.0.
south_north: float, optional
south north location of the measurement. Defaults to 0.0.
crs : int, dict, str or pyproj.crs.CRS, optional
Value to initialize `pyproj.crs.CRS`. Defaults to 4326.
time_col: int, str
column position (integer) or header (str) where the timestamp is located. it can be overriden by
using `pandas.read_csv` kwargs. Defaults to 0 (first column in the file).
height_to_columns: dict
dictionary to map the wind speed and directions to its corresponding height. The key is a float
with the height, and the value is a tuple (str,str) with the header for the wind speed and the
header for the wind direction, respectively. If the parameter is `None`, the columns are inferred
from the column names in the files. The function will find wind speeds for different heights and
after that will look for wind direction columns, matching them to the closest height.
Examples of autodetected header formats:
- ws_10, ws_10_mean, ws10, WS10 (wind speed at 10 m)
- windagl10, windagl_10, windagl_10_mean (wind speed at 10 m)
- wd_15, wd_15_mean, w15, WD15 (wind direction at 15m)
- wdiragl15, wdiragl_15, wdiragl_15_mean (wind direction at 15 m)
pandas_kwargs: None, dict, optional
Optional arguments that are forwarded to `pandas.read_csv` for customizing its behavior.
Returns
-------
da: xarray.Dataset
Time series wind climate dataset with variables 'wind_speed' and 'wind_direction'
and with a coordinate and dimension 'time'.
Raises
------
RuntimeError
If the time column cannot be parsed or if the wind speed and wind direction columns cannot
be detected.
"""
if pandas_kwargs is None:
pandas_kwargs = {}
default_pandas_kwargs = {
"index_col": time_col,
}
pandas_kwargs = {**default_pandas_kwargs, **pandas_kwargs}
pd_df = pd.read_csv(file_name, **pandas_kwargs)
if not isinstance(pd_df.index, pd.DatetimeIndex) and pd.api.types.is_string_dtype(
pd_df.index
):
try:
pd_df.index = pd.to_datetime(pd_df.index)
except (ValueError, pd.errors.ParserError) as e:
logger.debug("Could not auto-parse datetime index: %s", e)
return tswc_from_dataframe(pd_df, west_east, south_north, crs, height_to_columns)
[docs]
def tswc_from_dataframe(
pd_df,
west_east,
south_north,
crs,
height_to_columns=None,
):
"""
transforms a pandas.DataFrame into a time series wind climate xarray.Dataset. The dataframe must have
an index with time format and at least one wind speed and one wind direction. It allows to create a
dataset for several heights.
Parameters
----------
pd_df : pandas.DataFrame
pandas dataframe with wind speed and wind direction measurements for different timestamps and
heights.
west_east: float
west east locaton of the measurement
south_north: float
south north location of the measurement
crs : int, dict, str or pyproj.crs.CRS
Value to initialize `pyproj.crs.CRS`
height_to_columns: dict
dictionary to map the wind speed and directions to its corresponding height. The key is a float
with the height, and the value is a tuple (str,str) with the header for the wind speed and the
header for the wind direction, respectively. If the parameter is `None`, the columns are inferred
from the column names in the dataframe. The function will find wind speeds for different heights
and after that will look for wind direction columns, matching them to the closest height.
Examples of autodetected header formats:
- ws_10, ws_10_mean, ws10, WS10 (wind speed at 10 m)
- windagl10, windagl_10, windagl_10_mean (wind speed at 10 m)
- wd_15, wd_15_mean, w15, WD15 (wind direction at 15m)
- wdiragl15, wdiragl_15, wdiragl_15_mean (wind direction at 15 m)
Returns
-------
da: xarray.Dataset
Time series wind climate dataset with variables 'wind_speed' and 'wind_direction'
and with a coordinate and dimension 'time'.
Raises
------
RuntimeError
If it fails to autodetect the columns
"""
# Check if index is datetime
if not isinstance(pd_df.index, pd.DatetimeIndex):
raise RuntimeError(
"The dataframe index is not of type 'datetime'. Please provide a pandas.DataFrame with the time as index."
)
if height_to_columns is None:
try:
height_to_columns = _headers_to_dict(pd_df)
rows = ["{:<5} {:<12} {:<12}".format("h", "Wind speed", "Wind dir")]
rows += [
"{:<5} {:<12} {:<12}".format(k, v[0], v[1])
for k, v in height_to_columns.items()
]
logger.info("Columns detected\n%s", "\n".join(rows))
except Exception as err:
raise RuntimeError(
str(err)
+ "\nColumns could not be detected automatically. Provide a height_to_columns dictionary."
)
ds_pieces = []
for k, v in height_to_columns.items():
ws = xr.DataArray(pd_df[v[0]], dims=["time"])
wd = xr.DataArray(pd_df[v[1]], dims=["time"])
ds_piece = xr.Dataset({"wind_speed": ws, "wind_direction": wd}).assign_coords(
height=k,
)
ds_pieces.append(ds_piece)
ds = xr.concat(
ds_pieces, dim="height", data_vars="all", compat="equals", join="outer"
)
ds = ds.assign_coords(
{
"west_east": west_east,
"south_north": south_north,
}
)
ds = ds.transpose("time", ...)
ds = set_crs(ds, crs)
ds = _update_history(ds)
return to_stacked_point(_update_var_attrs(ds, {**_TS_ATTRS}))
def _headers_to_dict(df):
"""
Tries to detect the wind speed and wind direction columns on a pandas dataframe
and builds a dictionary
Parameters
----------
df : pandas.DataFrame
dataframe with wind data
Returns
-------
dict: Dictionary where the key is the height (float) and the value is a tuple (string,string)
with the wind speed column header and the wind direction column header, or None if the headers
format does not allow autodetection.
"""
headers_list = df.columns
vel_columns, dir_columns = _parse_heights_from_headers(headers_list)
if len(vel_columns) == 0:
raise RuntimeError("No wind speed columns were found.")
if len(dir_columns) == 0:
raise RuntimeError("No wind direction columns were found.")
v_list = np.array([x[0] for x in vel_columns])
d_list = np.array([x[0] for x in dir_columns])
duplicate_heights = [
x for x, count in collections.Counter(v_list).items() if count > 1
]
if len(duplicate_heights) > 0:
duplicate_heights_err_msg = ""
for val in duplicate_heights:
duplicate_heights_err_msg += (
f"There are duplicate entries for height {val}\n"
)
raise RuntimeError(duplicate_heights_err_msg.rstrip())
indices = _closest_dir_index(v_list, d_list)
final_dict = {}
for i, val in enumerate(vel_columns):
final_dict.update({val[0]: (val[1], dir_columns[indices[i]][1])})
return final_dict
def _parse_heights_from_headers(header_list):
"""Detect the wind speed and wind direction columns from a header list
Parameters
----------
header_list : list of string
list with each header
Returns
-------
list: list with tuples (float,string) with the height and the string header of the wind speed
list: list with tuples (float,string) with the height and the string header of the wind direction
"""
response_h = []
response_d = []
velocity_patterns_list = [
r"(?:(?:windagl)|(?:ws))_*(?P<height>\d*\.*\d+).*(?:mean)*",
r"a(?P<height>\d*\.*\d+)(:?(:?|:?T0deg))_wind_speed_mean",
]
direction_patterns_list = [
r"(?:(?:wdiragl)|(?:wd))_*(?P<height>\d*\.*\d+).*(?:mean)*",
r"d(?P<height>\d*\.*\d+)(:?(:?|:?T0deg))_wind_direction_mean",
]
while len(velocity_patterns_list) != 0:
velocity_pattern = velocity_patterns_list.pop(0)
direction_pattern = direction_patterns_list.pop(0)
for val in header_list:
match_vel = re.match(velocity_pattern, val, re.IGNORECASE)
match_dir = re.match(direction_pattern, val, re.IGNORECASE)
if match_vel is not None:
height_vel = match_vel.group("height")
response_h.append((float(height_vel), val))
if match_dir is not None:
height_dir = match_dir.group("height")
response_d.append((float(height_dir), val))
return response_h, response_d
def _closest_dir_index(vel_list, dir_list):
"""
returns a list with the indices with the closest value of wind
direction for a given wind velocity
Parameters
----------
vel_list : numpy.array
array with heights where the velocity was measured
dir_list : numpy.array
array with height where the direction was measured
Returns
-------
list : list
list with the indices in dir_list corresponding to vel_list
"""
resp = []
for val in vel_list:
resp.append(np.argmin(abs(val - dir_list)))
return resp
def _read_windiofile(windiofile, crs=None, level="plant", **kwargs):
"""Read a windIO YAML file into a time series wind climate dataset.
Parameters
----------
windiofile : str or pathlib.Path
Path to windIO YAML file.
crs : int, dict, str or CRS
Coordinate reference system. Required if not present in windIO file.
level : str
The windIO hierarchy level of the input file: 'plant', 'site',
'energy_resource', or 'wind_resource'. Defaults to 'plant'.
**kwargs
Additional keyword arguments (unused, for API consistency).
Returns
-------
xarray.Dataset
Time series wind climate dataset.
Raises
------
ValueError
If CRS is not provided and not present in windIO file.
ValueError
If level is not valid.
"""
from windkit.io.windio.translate import tswc_from_windio
from windkit.io.windio.utils import _read_windio_climate
return _read_windio_climate(windiofile, tswc_from_windio, crs=crs, level=level)
[docs]
def read_tswc(filename, *, crs=None, file_format="infer", **kwargs):
"""
Reads a time series wind climate file into a xarray.Dataset.
Parameters
----------
filename : str
path to the file to be read.
crs : int, dict, str or CRS, None, optional
Coordinate reference system. Required for windIO YAML files.
file_format : str
format of the file to be read. If 'infer' is passed, the function will try to infer the format
from the file extension. Defaults to 'infer'.
**kwargs : dict
Additional arguments that are forwarded to the specific file reader.
For windIO YAML files, the 'level' argument can be provided here to specify the windIO hierarchy level
of the input file: 'plant', 'site', 'energy_resource', or 'wind_resource'. Defaults to 'plant'.
Returns
-------
xarray.Dataset
Time series wind climate dataset.
"""
if file_format == "infer":
file_format = _infer_file_format(filename)
if file_format == "csv":
ds = _read_csv(filename, **kwargs)
elif file_format == "txt":
ds = _read_ts_windpro_txt(filename, **kwargs)
elif file_format == "nc":
ds = xr.open_dataset(filename, **kwargs)
elif file_format == "yaml":
ds = _read_windiofile(filename, crs=crs, **kwargs)
else:
raise NotImplementedError(f"File format {file_format} is not supported.")
ds = _update_var_attrs(ds, _TS_ATTRS)
return _update_history(ds)
def _to_windiofile(tswc, windiofile, level="plant", **kwargs):
"""Write a time series wind climate dataset to windIO YAML format.
Parameters
----------
tswc : xarray.Dataset
Time series wind climate dataset.
windiofile : str or pathlib.Path
Path to output windIO YAML file.
level : str
Export level: 'plant', 'site', 'energy_resource', or 'wind_resource'.
Defaults to 'plant'.
**kwargs
Additional keyword arguments (unused, for API consistency).
"""
from windkit.io.windio.translate import tswc_to_windio
from windkit.io.windio.utils import _write_windio_climate
_write_windio_climate(tswc, windiofile, tswc_to_windio, level=level)
[docs]
@validate_input(_TSWC_SCHEMA)
def tswc_to_file(tswc, filename, *, file_format="infer", **kwargs):
"""Write time series wind climate to file.
Parameters
----------
tswc : xarray.Dataset
Time series wind climate dataset.
filename : str or pathlib.Path
Path to the file to write the time series wind climate to.
file_format : str
The file format of the output file. If 'infer', the file format will be
inferred from the file extension. Supported formats are 'yaml'.
Defaults to 'infer'.
**kwargs : dict
Additional keyword arguments to pass to the writing function.
For windIO YAML files, the 'level' argument can be provided here to specify the windIO hierarchy level
of the output file: 'plant', 'site', 'energy_resource', or 'wind_resource'. Defaults to 'plant'.
"""
if file_format == "infer":
file_format = _infer_file_format(filename)
if file_format in ["yml", "yaml"]:
_to_windiofile(tswc, filename, **kwargs)
else:
raise ValueError(f"Unsupported file format for TSWC: {file_format}")
def _sample_groups(values, wd, n_sectors):
"""Group of each sample along the last axis.
The group is the sample's wind direction sector, or 0 when ``n_sectors``
is None. Samples with a missing value, or a missing direction when
grouping by sector, get -1.
"""
usable = ~np.isnan(values)
groups = np.full(values.shape, -1, dtype=np.int32)
if n_sectors is None:
groups[usable] = 0
return groups
sector, _ = wd_to_sector(wd, n_sectors, output_type="indices")
usable &= ~np.isnan(sector)
np.copyto(groups, sector, casting="unsafe", where=usable)
return groups
def _group_index(groups, n_groups):
"""Flat index into a (batch, group) array of every usable sample.
Also returns the mask of usable samples, as a (batch, time) array.
"""
n_batch = int(np.prod(groups.shape[:-1], dtype=int))
groups = groups.reshape(n_batch, -1)
usable = groups >= 0
offsets = np.arange(n_batch, dtype=np.intp)[:, np.newaxis] * n_groups
flat = np.broadcast_to(offsets, groups.shape)[usable]
flat += groups[usable]
return flat, usable
def _ratio(num, den):
"""Divide, giving NaN where ``den`` is zero (a group without samples)."""
out = np.full(np.broadcast_shapes(num.shape, den.shape), np.nan)
return np.divide(num, den, out=out, where=den > 0)
def _sector_mean(values, wd, *, n_sectors):
"""Mean of ``values`` over the last axis in each sector."""
groups = _sample_groups(values, wd, n_sectors)
flat, usable = _group_index(groups, n_sectors)
shape = (*groups.shape[:-1], n_sectors)
size = int(np.prod(shape, dtype=int))
kept = np.reshape(values, usable.shape)[usable]
sums = np.bincount(flat, weights=kept, minlength=size)
return _ratio(sums, np.bincount(flat, minlength=size)).reshape(shape)
def _sector_freq_gt_mean(ws, wd, *, n_sectors):
"""Fraction of samples above the mean of their sector."""
groups = _sample_groups(ws, wd, n_sectors)
flat, usable = _group_index(groups, n_sectors)
shape = (*groups.shape[:-1], n_sectors)
size = int(np.prod(shape, dtype=int))
kept = np.reshape(ws, usable.shape)[usable]
counts = np.bincount(flat, minlength=size)
means = _ratio(np.bincount(flat, weights=kept, minlength=size), counts)
above = np.bincount(flat, weights=kept > means[flat], minlength=size)
return _ratio(above, counts).reshape(shape)
def _grouped_cdf(ws, wd=None, *, wsceil, n_sectors=None):
"""Fraction of samples below each of ``wsceil``, per sector if ``n_sectors``."""
groups = _sample_groups(ws, wd, n_sectors)
n_groups = n_sectors or 1
n_bins = wsceil.size
flat, usable = _group_index(groups, n_groups)
# A sample counts towards the CDF from the first ceiling above it on
first_bin = np.searchsorted(
wsceil, np.reshape(ws, usable.shape)[usable], side="right"
)
flat *= n_bins + 1
flat += first_bin
shape = (*groups.shape[:-1], n_groups, n_bins + 1)
counts = np.bincount(flat, minlength=int(np.prod(shape, dtype=int)))
counts = counts.reshape(shape)
below = np.cumsum(counts[..., :n_bins], axis=-1)
return _ratio(below, counts.sum(axis=-1, keepdims=True))
def _grouped_stat(kernel, values, tswc, n_sectors, extra_dims=None, **kwargs):
"""Apply a grouped kernel over time, per sector unless ``n_sectors`` is None.
Lazy for dask-backed input, whose time chunks are merged into one. Recent
dask splits the other dimensions to keep the blocks' size; older versions
do not, so chunk those dimensions too if a block would not fit in memory.
"""
extra_dims = extra_dims or {}
args = [values] if n_sectors is None else [values, tswc[WD]]
result = xr.apply_ufunc(
kernel,
*args,
kwargs={"n_sectors": n_sectors, **kwargs},
input_core_dims=[[DIM_TIME]] * len(args),
output_core_dims=[[DIM_SECTOR, *extra_dims]],
dask="parallelized",
output_dtypes=[np.float64],
dask_gufunc_kwargs={
"output_sizes": {DIM_SECTOR: n_sectors or 1, **extra_dims},
"allow_rechunk": True,
},
)
if n_sectors is None:
return result.isel({DIM_SECTOR: 0}, drop=True)
return result.assign_coords(create_sector_coords(n_sectors).coords)
def _time_mean(values, tswc, bysector, n_sectors):
"""Mean of ``values`` over time, per sector if ``bysector``.
Accumulated in float64. Samples with a missing value, or a missing wind
direction when ``bysector`` is True, are skipped. A sector without samples
gives NaN.
"""
if bysector:
return _grouped_stat(_sector_mean, values, tswc, n_sectors)
return values.mean(dim=DIM_TIME, dtype=np.float64)
def _mean_ws_moment(tswc, moment=1, *, bysector=False, n_sectors=12):
"""Calculate the mean wind speed moment of a time series wind climate."""
return _time_mean(tswc[WS] ** moment, tswc, bysector, n_sectors)
def _ws_cdf(tswc, *, bysector=False, ws_bins=None, n_sectors=12):
"""Calculate the wind speed CDF of a time series wind climate.
The CDF at each ``wsceil`` is the fraction of samples strictly below it.
Samples are skipped as in :func:`_time_mean`.
"""
wsbins = create_wsbin_coords() if ws_bins is None else create_wsbin_coords(ws_bins)
wsceil = wsbins["wsceil"].values
result = _grouped_stat(
_grouped_cdf,
tswc[WS],
tswc,
n_sectors if bysector else None,
extra_dims={DIM_WSBIN: wsceil.size},
wsceil=wsceil,
)
return result.assign_coords(wsbins.coords)
def _ws_freq_gt_mean(tswc, *, bysector=False, n_sectors=12):
"""Calculate the frequency of wind speeds greater than the mean wind speed.
Samples are skipped as in :func:`_time_mean`.
"""
if bysector:
return _grouped_stat(_sector_freq_gt_mean, tswc[WS], tswc, n_sectors)
ws = tswc[WS]
n_valid = ws.notnull().sum(dim=DIM_TIME)
n_above = (ws > ws.mean(dim=DIM_TIME, dtype=np.float64)).sum(dim=DIM_TIME)
return n_above / n_valid.where(n_valid > 0)
def _mean_wind_speed(tswc, *, bysector=False, n_sectors=12):
"""Calculate the mean wind speed of a time series wind climate."""
return _mean_ws_moment(tswc, 1, bysector=bysector, n_sectors=n_sectors)
def _mean_power_density(tswc, *, air_density=1.225, bysector=False, n_sectors=12):
"""Calculate the mean power density of a time series wind climate."""
return _time_mean(0.5 * air_density * tswc[WS] ** 3, tswc, bysector, n_sectors)
def _extrapolate_to_height(tswc, height, shear_exponent=0.143):
"""Extrapolates wind speed and wind direction to a given height using power law for wind speed"""
ws_extrapolated = shear_extrapolate(
tswc["wind_speed"], height, shear_exponent=shear_exponent
)
wd_extrapolated = tswc["wind_direction"].interp(
height=height, method="nearest", kwargs={"fill_value": "extrapolate"}
)
result = xr.Dataset(
{
"wind_speed": ws_extrapolated,
"wind_direction": wd_extrapolated,
}
).assign_attrs(tswc.attrs)
result = result.transpose(*result.wind_direction.dims, ...)
return result
