"""Wind resources, the models of each case they sample, and the cases.
A resource gives the direction, speed, surface, stability and turbulence of a
site, and :meth:`WindResource.sample` draws :class:`SiteConditions`.
::
from oowl import resources as R
resource = R.WindResource(
direction=R.VonMises(seed=2),
mean_wind_speed=9.0,
roughness_length=R.Charnock(),
stability=R.KellyGryning(unstable_fraction=0.4),
turbulence=R.SurfaceLayer(coefficient=2.3),
mann=R.Mann(),
)
cases = resource.sample(1024, seed=1)
Every model argument is optional, and the defaults describe open water.
:class:`GWA` takes the directions and speeds of a site from the Global Wind
Atlas, and :class:`Dimitrov` is the coverage design for training data.
``R.profile`` holds the Monin-Obukhov profile functions the models use.
"""
from __future__ import annotations
import inspect
import queue
import threading
import time
import typing
import urllib.error
import urllib.request
from abc import ABC, abstractmethod
from dataclasses import KW_ONLY, dataclass, field, fields, replace
from functools import cached_property
from pathlib import Path
from typing import TYPE_CHECKING, Any, ClassVar, Sequence, cast
import numpy as np
if TYPE_CHECKING: # pragma: no cover
import matplotlib.axes
from numpy.typing import ArrayLike
from ._random import as_rng, resolve_seed
from ._repr import _ArgumentRepr, _describe, _FieldRepr
from ._site import profile
from ._site.base import SiteConditions, sector_angles
from ._site.coverage import Inflow, Roughness, Stability, Turbulence
from ._site.coverage import Mann as MannCoverage
from ._site.sampling import U_MAX, U_MIN, _checked_speed_range
from ._site.stability import Charnock, KellyGryning, friction_velocity
from ._site.turbulence import Mann, SurfaceLayer
from ._validation import (
_validate_finite,
_validate_integer,
_validate_positive,
_validate_real,
_validate_real_array,
_validate_unit_interval,
)
__all__ = [
"WindResource",
"GWA",
"Dimitrov",
"UniformDirection",
"SingleDirection",
"VonMises",
"Elliptical",
"SiteConditions",
"KellyGryning",
"Charnock",
"SurfaceLayer",
"Mann",
"profile",
]
def _validate_finite_pair(name: str, values) -> tuple[float, float]:
pair = _validate_real_array(name, values)
if pair.shape != (2,) or not np.all(np.isfinite(pair)):
raise ValueError(f"{name} must contain two finite values")
return float(pair[0]), float(pair[1])
def _validate_bins(name: str, bins) -> np.ndarray:
bins = _validate_real_array(name, bins)
if bins.ndim != 1 or len(bins) < 2:
raise ValueError(f"{name} must be a 1-D array with at least two edges")
if not np.all(np.isfinite(bins)) or not np.all(np.diff(bins) > 0):
raise ValueError(f"{name} must be finite and strictly increasing")
return bins
def _validate_bins_cover_range(
name: str,
bins: np.ndarray,
lower: float,
upper: float,
) -> None:
tol = 1e-12
if bins[0] > lower + tol or bins[-1] < upper - tol:
raise ValueError(
f"{name} must cover [{lower}, {upper}] when require_full_range=True"
)
class _WeibullSpeed:
"""Per-sector Weibull speeds, the distribution of a site.
``cdf`` and ``pdf`` take speeds and return one row per sector. ``ppf``
draws speeds truncated to ``[u_min, u_max]`` in the sectors given.
``time_share`` says whether the cases represent a share of time.
"""
time_share = True
def __init__(self, scale: np.ndarray, shape: np.ndarray):
self.scale = np.asarray(scale, dtype=float)
self.shape = np.asarray(shape, dtype=float)
# A large k overflows (u / a)^k to infinity above the scale, where the CDF
# is one and the density zero, so those overflows give the right limits.
def cdf(self, u, u_min: float, u_max: float) -> np.ndarray:
a, k = self.scale[:, None], self.shape[:, None]
x = np.clip(np.asarray(u, dtype=float)[None, :] / a, 0.0, None)
with np.errstate(over="ignore"):
return 1.0 - np.exp(-(x**k))
def pdf(self, u, u_min: float, u_max: float) -> np.ndarray:
a, k = self.scale[:, None], self.shape[:, None]
x = np.clip(np.asarray(u, dtype=float)[None, :] / a, 0.0, None)
# x^(k - 1) overflows only where exp(-x^k) is already zero.
with np.errstate(over="ignore", invalid="ignore"):
tail = np.exp(-(x**k))
density = (k / a) * x ** (k - 1.0) * tail
return np.where(tail > 0.0, density, 0.0)
def log_pdf(self, u) -> np.ndarray:
"""Return the log density at positive speeds, one row per sector."""
a, k = self.scale[:, None], self.shape[:, None]
x = np.asarray(u, dtype=float)[None, :] / a
with np.errstate(over="ignore"):
return np.log(k / a) + (k - 1.0) * np.log(x) - x**k
def ppf(self, q, sector, u_min: float, u_max: float) -> np.ndarray:
q = np.asarray(q, dtype=float)
a, k = self.scale[sector], self.shape[sector]
with np.errstate(over="ignore"):
f_min = 1.0 - np.exp(-(np.clip(u_min / a, 0.0, None) ** k))
f_max = 1.0 - np.exp(-(np.clip(u_max / a, 0.0, None) ** k))
f = f_min + q * (f_max - f_min)
speed = a * (-np.log1p(-np.clip(f, 0.0, 1.0 - 1e-15))) ** (1.0 / k)
# A range holding less than sqrt(eps) of a sector's time keeps under
# half the digits of its CDF difference, so its speeds are inverted in
# the cumulative hazard y = (u / a)^k instead, where
# y = y_min - log1p(q expm1(y_min - y_max)) stays accurate.
faint = f_max - f_min < np.sqrt(np.finfo(float).eps)
if np.any(faint):
a, k = a[faint], k[faint]
with np.errstate(over="ignore"):
y_min, y_max = (u_min / a) ** k, (u_max / a) ** k
y = y_min - np.log1p(q[faint] * np.expm1(y_min - y_max))
speed[faint] = a * y ** (1.0 / k)
return speed
def mean(self) -> np.ndarray:
"""Return each sector's mean speed, ``A Gamma(1 + 1/k)``."""
from scipy.special import gamma
return self.scale * gamma(1.0 + 1.0 / self.shape)
def second_moment(self) -> np.ndarray:
"""Return each sector's mean squared speed, ``A^2 Gamma(1 + 2/k)``."""
from scipy.special import gamma
return self.scale**2 * gamma(1.0 + 2.0 / self.shape)
class _TabulatedSpeed:
"""Per-sector speeds from a table over wind speed, as windIO stores them.
Each speed carries the probability of the cell around it, bounded by the
midpoints to its neighbors and by the first and last speed, and the
speed is uniform within a cell. ``lowest`` moves the first cell's lower
bound below the first speed, for counts that hold slower wind. A table of
one speed and no lower bound below it is a point mass, so every case draws
that speed.
"""
time_share = True
def __init__(self, speeds, masses, lowest: float | None = None):
speeds = np.atleast_1d(np.asarray(speeds, dtype=float))
masses = np.asarray(masses, dtype=float)
if speeds.ndim != 1 or len(speeds) < 1 or np.any(np.diff(speeds) <= 0):
raise ValueError(
"a probability table needs one or more increasing wind speeds"
)
if masses.ndim != 2 or masses.shape[1] != len(speeds):
raise ValueError("a probability table needs one column per wind speed")
if not np.all(np.isfinite(masses)) or np.any(masses < 0):
raise ValueError("probabilities must be finite and non-negative")
self.speeds = speeds
first = speeds[:1] if lowest is None else [min(float(lowest), speeds[0])]
self.edges = np.concatenate(
[first, 0.5 * (speeds[1:] + speeds[:-1]), speeds[-1:]]
)
totals = masses.sum(axis=1, keepdims=True)
# A sector without wind is never drawn, so any shape serves it.
masses = np.where(
totals > 0,
masses / np.where(totals > 0, totals, 1.0),
1.0 / masses.shape[1],
)
self.masses = masses
self.cumulative = np.concatenate(
[np.zeros((len(masses), 1)), np.cumsum(masses, axis=1)], axis=1
)
@property
def point_mass(self) -> bool:
"""Whether the table holds one speed only, in a cell of no width."""
return len(self.speeds) == 1 and self.edges[0] == self.edges[-1]
def cdf(self, u, u_min: float, u_max: float) -> np.ndarray:
u = np.asarray(u, dtype=float)
if self.point_mass:
# The sampled range is closed at u_min, so the CDF there is its
# left limit and a point mass on u_min still lies in the range.
reached = np.where(u == u_min, u > self.speeds[0], u >= self.speeds[0])
return np.repeat(reached[None, :].astype(float), len(self.masses), axis=0)
return np.array([np.interp(u, self.edges, row) for row in self.cumulative])
def pdf(self, u, u_min: float, u_max: float) -> np.ndarray:
u = np.asarray(u, dtype=float)
if self.point_mass:
return self._point_mass_pdf(u)
cell = np.clip(
np.searchsorted(self.edges, u, side="right") - 1, 0, len(self.speeds) - 1
)
inside = (u >= self.edges[0]) & (u <= self.edges[-1])
density = self.masses / np.diff(self.edges)[None, :]
return np.where(inside[None, :], density[:, cell], 0.0)
def _point_mass_pdf(self, u: np.ndarray) -> np.ndarray:
"""Return the point mass as a density on the quadrature nodes ``u``.
A point mass has no density, so it is spread over the two nodes around
its speed as the hat function there, divided by the trapezoid weights
of those nodes. The trapezoid rule over ``u`` then integrates it to one
and integrates its product with a function to that function's linear
interpolant at the speed.
"""
speed = self.speeds[0]
density = np.zeros((len(self.masses), len(u)))
if len(u) < 2 or speed < u[0] or speed > u[-1]:
return density
upper = int(np.clip(np.searchsorted(u, speed, side="right"), 1, len(u) - 1))
fraction = (speed - u[upper - 1]) / (u[upper] - u[upper - 1])
weights = np.empty(len(u))
weights[0], weights[-1] = 0.5 * (u[1] - u[0]), 0.5 * (u[-1] - u[-2])
weights[1:-1] = 0.5 * (u[2:] - u[:-2])
density[:, upper - 1] = (1.0 - fraction) / weights[upper - 1]
density[:, upper] = fraction / weights[upper]
return density
def ppf(self, q, sector, u_min: float, u_max: float) -> np.ndarray:
q = np.asarray(q, dtype=float)
if self.point_mass:
return np.full(q.shape, self.speeds[0])
sector = np.asarray(sector, dtype=int)
speeds = np.empty(q.shape)
for index in np.unique(sector):
rows = sector == index
low, high = np.interp([u_min, u_max], self.edges, self.cumulative[index])
target = low + q[rows] * (high - low)
speeds[rows] = np.interp(target, self.cumulative[index], self.edges)
return speeds
def mean(self) -> np.ndarray:
"""Return each sector's mean speed, the centers of its cells weighted."""
return self.masses @ (0.5 * (self.edges[:-1] + self.edges[1:]))
def second_moment(self) -> np.ndarray:
"""Return each sector's mean squared speed, uniform within each cell."""
low, high = self.edges[:-1], self.edges[1:]
return self.masses @ ((low**2 + low * high + high**2) / 3.0)
class _UniformSpeed:
"""Speeds spread evenly over the sampled range in every sector, a design."""
time_share = False
def __init__(self, n_sectors: int):
self.n_sectors = int(n_sectors)
def cdf(self, u, u_min: float, u_max: float) -> np.ndarray:
u = np.asarray(u, dtype=float)
share = np.clip((u - u_min) / (u_max - u_min), 0.0, 1.0)
return np.repeat(share[None, :], self.n_sectors, axis=0)
def pdf(self, u, u_min: float, u_max: float) -> np.ndarray:
u = np.asarray(u, dtype=float)
density = ((u >= u_min) & (u <= u_max)) / (u_max - u_min)
return np.repeat(density[None, :], self.n_sectors, axis=0)
def ppf(self, q, sector, u_min: float, u_max: float) -> np.ndarray:
return u_min + np.asarray(q, dtype=float) * (u_max - u_min)
def _mean_over_speed(distribution, wind_speed: np.ndarray, values) -> np.ndarray:
"""Return each sector's mean of ``values`` over ``wind_speed``, by quadrature.
``values`` has one row per sector and one column per speed, and any further
axes are averaged alike. The speeds are weighted by each sector's density
over them, normalized on the same nodes, so the weights sum to one however
coarsely the nodes resolve the density. A sector whose wind lies far from
the speeds has a Weibull density that underflows there, so its weights are
taken in logs, which keeps them peaked at the end nearest its wind. A
tabulated sector whose cells all lie outside the speeds, or between two of
them, has no density on them, so it takes the values at the speed nearest
its mean.
"""
pdf = distribution.pdf(wind_speed, wind_speed[0], wind_speed[-1])
# Below this peak the products with the values would turn subnormal.
faint = np.max(pdf, axis=1) < np.sqrt(np.finfo(float).tiny)
if np.any(faint) and isinstance(distribution, _WeibullSpeed):
log_pdf = distribution.log_pdf(wind_speed)[faint]
pdf[faint] = np.exp(log_pdf - log_pdf.max(axis=1, keepdims=True))
elif np.any(faint) and isinstance(distribution, _TabulatedSpeed):
offset = np.abs(wind_speed[None, :] - distribution.mean()[faint, None])
pdf[faint] = np.arange(len(wind_speed)) == np.argmin(offset, axis=1)[:, None]
pdf = pdf.reshape(pdf.shape + (1,) * (np.ndim(values) - 2))
return np.trapezoid(pdf * values, wind_speed, axis=1) / np.trapezoid(
pdf, wind_speed, axis=1
)
def _per_sector(name: str, value, n_sectors: int) -> np.ndarray:
"""Return positive finite ``value`` as one entry per sector."""
values = _validate_real_array(name, value)
if values.ndim == 0:
values = np.full(n_sectors, float(values), dtype=float)
if values.shape != (n_sectors,):
raise ValueError(f"{name} must be scalar or have length n_sectors")
if not np.all(np.isfinite(values)) or np.any(values <= 0):
raise ValueError(f"{name} must contain only positive finite values")
return values
class _Resource(ABC):
"""The chain of models every case is drawn along, shared by all resources.
A subclass sets the sectors and their probabilities, the speed
distribution and the models of each case, then calls :meth:`_setup_chain`.
"""
n_sectors: int
theta: np.ndarray
probability: np.ndarray
reference_height: float
roughness_length: float | Charnock | Roughness | Sequence[float | Charnock]
stability: KellyGryning | Stability
turbulence: SurfaceLayer | Turbulence
mann: Mann | MannCoverage | None
# Veer, yaw, tilt and air density, drawn only by a Dimitrov design.
_inflow: Inflow | None = None
_speed: "_WeibullSpeed | _UniformSpeed | _TabulatedSpeed"
def _setup_chain(self) -> None:
"""Check the surface against the height and each sector's TI level."""
self.reference_height = _validate_real(
"reference_height", self.reference_height
)
if not np.isfinite(self.reference_height) or self.reference_height <= 0:
raise ValueError("reference_height must be positive and finite")
# A Charnock z0 always lies below the height it is solved at.
for surface, _ in self._sector_surfaces():
if isinstance(surface, Charnock):
continue
roughest = surface.maximum if isinstance(surface, Roughness) else surface
if self.reference_height <= roughest:
raise ValueError("reference_height must lie above roughness_length")
# A set level scales each sector by :attr:`_turbulence_scale`, computed
# on first use as it averages the whole chain over the speed climate.
self._target_intensity = None
target = getattr(self.turbulence, "intensity", None)
if target is not None:
values = _validate_real_array("turbulence intensity", target)
if values.ndim == 0:
values = np.full(self.n_sectors, float(values))
if values.shape != (self.n_sectors,):
raise ValueError(
"turbulence intensity must be one value or one per sector"
)
self._target_intensity = values
@cached_property
def _turbulence_scale(self) -> np.ndarray:
"""Per-sector factor on ``sigma_u``, one unless the model sets a level."""
if self._target_intensity is None:
return np.ones(self.n_sectors)
return self._target_intensity / self._infer_sector_ti()
def __eq__(self, other: object) -> bool:
"""Whether ``other`` is a resource of the same class and record."""
if other.__class__ is not self.__class__:
return NotImplemented
return self.metadata() == other.metadata()
@abstractmethod
def metadata(self) -> dict:
"""Return a JSON-compatible record of the resource."""
@abstractmethod
def _speed_record(self, wind_speed: np.ndarray) -> dict:
"""Return the windIO fields of the direction and speed distribution."""
def _roughness_at(
self, wind_speed: ArrayLike, sector: ArrayLike | None = None
) -> np.ndarray:
"""Return the roughness length in meters at each ``reference_height`` speed.
``sector`` gives each speed's sector index, needed when the roughness
differs by sector.
"""
speed = np.asarray(wind_speed, dtype=float)
if not isinstance(self.roughness_length, tuple):
return self._surface_roughness(self.roughness_length, speed)
if sector is None:
raise ValueError("a roughness per sector needs the sector of each speed")
index = np.broadcast_to(np.asarray(sector, dtype=int), speed.shape)
lengths = np.empty(speed.shape)
for surface, rows in self._sector_surfaces():
cases = rows[index]
if np.any(cases):
lengths[cases] = self._surface_roughness(surface, speed[cases])
return lengths
def _surface_roughness(self, surface, speed: np.ndarray) -> np.ndarray:
"""Return the roughness of one surface at each speed.
A :class:`Dimitrov` roughness range is drawn per case, so its median
stands for it in summaries such as the TI rose and the windIO ``z0``.
"""
if isinstance(surface, Charnock):
return surface.roughness_length(speed, self.reference_height)
if isinstance(surface, Roughness):
return np.full(speed.shape, surface.median)
return np.full(speed.shape, surface)
def _sector_surfaces(self) -> list[tuple[Any, np.ndarray]]:
"""Return each distinct surface with a mask of the sectors it covers."""
roughness = self.roughness_length
if not isinstance(roughness, tuple):
return [(roughness, np.ones(self.n_sectors, dtype=bool))]
groups: list[tuple[Any, np.ndarray]] = []
for index, surface in enumerate(roughness):
for known, rows in groups:
if known == surface:
rows[index] = True
break
else:
rows = np.zeros(self.n_sectors, dtype=bool)
rows[index] = True
groups.append((surface, rows))
return groups
def _validate_probability(self) -> None:
if self.probability.shape != (self.n_sectors,):
raise ValueError("sector probability must have length n_sectors")
if not np.all(np.isfinite(self.probability)):
raise ValueError("sector probability must be finite")
if np.any(self.probability < 0):
raise ValueError("sector probability must be non-negative")
total = float(self.probability.sum())
if total <= 0:
raise ValueError("sector probability must contain positive mass")
# Dividing an already normalized array by a sum a few ulps from one
# would nudge it again on every reconstruction, so leave it as given.
if abs(total - 1.0) > 8 * np.finfo(float).eps * self.n_sectors:
self.probability = self.probability / total
@property
def sector_turbulence_intensity(self) -> np.ndarray:
"""Return each sector's mean turbulence intensity over 4 to 25 m/s.
The mean weighs each speed by its probability in the sector, so a
sector whose wind lies far outside the range takes the turbulence
intensity at the nearer end.
"""
if self._target_intensity is not None:
return self._target_intensity.copy()
return self._infer_sector_ti()
def sample(
self,
n_cases: int,
*,
seed: int | np.random.Generator | None = None,
u_min: float = U_MIN,
u_max: float = U_MAX,
) -> SiteConditions:
"""Return ``n_cases`` inflow cases drawn with a scrambled Sobol sequence.
Each case draws direction, wind speed, roughness, stability,
``sigma_u`` and the optional Mann and inflow variables in turn, each
given the ones before it. The cases are conditional on a wind speed in
``[u_min, u_max]`` in m/s, and ``metadata["in_range_fraction"]`` is the
share of time they represent, ``None`` for a :class:`Dimitrov` design.
A given ``seed`` replays the cases, and ``None`` draws fresh entropy
that the metadata records.
"""
from ._site.sampling import generate_site_conditions
return generate_site_conditions(
self,
n_cases,
seed=seed,
u_min=u_min,
u_max=u_max,
)
def plot(
self,
ax: "matplotlib.axes.Axes | None" = None,
*,
variable: typing.Literal["wind_speed", "turbulence_intensity"] = "wind_speed",
**kwargs: Any,
) -> "matplotlib.axes.Axes":
"""Plot the analytical wind resource and return its Matplotlib axes.
``variable="wind_speed"`` draws the joint direction-speed rose.
``variable="turbulence_intensity"`` draws the analytical TI rose. Bars
are percent of all time, or of a :class:`Dimitrov` design's cases, as
:meth:`wind_rose` and :meth:`ti_rose` give.
"""
from .plotting import _plot_ti_rose, _plot_wind_rose
if variable == "wind_speed":
return _plot_wind_rose(self, ax=ax, **kwargs)
if variable == "turbulence_intensity":
return _plot_ti_rose(self, ax=ax, **kwargs)
raise ValueError("variable must be 'wind_speed' or 'turbulence_intensity'")
@property
def ti_bins(self) -> np.ndarray:
"""Default turbulence-intensity bin edges used by :meth:`ti_rose`."""
return self._ti_bins(U_MIN, U_MAX)
def _ti_bins(self, u_min: float, u_max: float, n_bins: int = 8) -> np.ndarray:
"""Return ``n_bins`` equal bins over the reachable turbulence intensity."""
return np.linspace(*self._ti_span(u_min, u_max), n_bins + 1)
def _ti_span(self, u_min: float, u_max: float) -> tuple[float, float]:
"""Return the turbulence intensities the model can reach over a range."""
u = np.linspace(u_min, u_max, 2000)
lows, highs = [], []
for surface, rows in self._sector_surfaces():
low, high = self._turbulence_summary("intensity_range", u, surface=surface)
scale = self._turbulence_scale[rows]
lows.append(low.min() * scale.min())
highs.append(high.max() * scale.max())
return float(min(lows)), float(max(highs))
def _turbulence_summary(self, name: str, wind_speed, *args, surface):
"""Return a turbulence model's summary ``name`` over one surface."""
u = np.asarray(wind_speed, dtype=float)
return getattr(self.turbulence, name)(
u,
*args,
height=self.reference_height,
roughness_length=self._surface_roughness(surface, u),
stability=self.stability,
)
def _in_range_masses(self, u_min: float, u_max: float) -> np.ndarray:
"""Return each sector's probability of a speed in ``[u_min, u_max]``.
A mass lost to rounding counts as zero. A calm sector may have no wind
in the range, and it is never drawn.
"""
ends = self._speed.cdf(np.array([u_min, u_max]), u_min, u_max)
mass = ends[:, 1] - ends[:, 0]
return np.where(mass > 1e-15, mass, 0.0)
def _in_range_probability(
self, *, u_min: float = U_MIN, u_max: float = U_MAX
) -> tuple[np.ndarray, float]:
"""Sector probabilities given a speed in ``[u_min, u_max]``, and its share.
Returns ``(p_i m_i / M, M)``, where ``m_i`` is sector ``i``'s share of
time in the range and ``M`` the share over all sectors, the fraction of
time the sampled cases represent. Sectors that are often too calm or too
strong weigh less. A range that no probable sector reaches cannot be
sampled.
"""
u_min, u_max = _checked_speed_range(u_min, u_max)
weighted = self.probability * self._in_range_masses(u_min, u_max)
total = float(weighted.sum())
if total <= 0.0:
raise ValueError("wind-speed truncation range has zero probability")
return weighted / total, total
def wind_rose(
self,
speed_bins: ArrayLike | None = None,
*,
u_min: float = U_MIN,
u_max: float = U_MAX,
require_full_range: bool = True,
) -> np.ndarray:
"""Return the share of time in each direction sector and speed bin.
The table has shape ``(n_sectors, n_speed_bins)``, and entry ``(i, b)``
is ``P(sector i, speed in bin b)``, the share of all time, calm
included. It sums to the share of time in ``[u_min, u_max]`` in m/s,
so annual energy is ``8760 h`` times the sum of the table times the
power in each bin. ``u_min=0`` covers calm air, and ``1 - sum`` is then
the share above ``u_max``. Divide by the sum for the distribution given
a speed in the range. A :class:`Dimitrov` design has no share of time,
so its rose is its distribution over the range and sums to one. The
default bins split the range into six, and custom bins must cover it
unless ``require_full_range=False``.
"""
u_min, u_max = _checked_speed_range(u_min, u_max, calm=True)
if speed_bins is None:
speed_bins = np.linspace(u_min, u_max, 7)
bins = _validate_bins("speed_bins", speed_bins)
if require_full_range:
_validate_bins_cover_range("speed_bins", bins, u_min, u_max)
lows = np.clip(bins[:-1], u_min, u_max)
highs = np.clip(bins[1:], u_min, u_max)
widths = np.maximum(highs - lows, 0.0)
bin_mass = self._speed.cdf(highs, u_min, u_max) - self._speed.cdf(
lows, u_min, u_max
)
return self.probability[:, None] * np.where(
widths[None, :] > 0.0, bin_mass, 0.0
)
def ti_rose(
self,
ti_bins: ArrayLike | None = None,
*,
u_min: float = U_MIN,
u_max: float = U_MAX,
n_quadrature: int = 2000,
require_full_range: bool = True,
) -> np.ndarray:
"""Return the share of time in each direction sector and TI bin.
Entry ``(i, b)`` is ``P(sector i, TI in bin b, speed in [u_min,
u_max])``, the share of all time, integrated over ``n_quadrature``
speeds. A rose whose bins cover the reachable TI sums to the share of
time in the speed range, as :meth:`wind_rose` does, and its rows sum
to the rows of :meth:`wind_rose`, a :class:`Dimitrov` design's to one.
Turbulence intensity is undefined in
calm air, so ``u_min`` must be above 0. The default bins split the
reachable TI into eight, and custom bins must cover it unless
``require_full_range=False``.
"""
u_min, u_max = _checked_speed_range(u_min, u_max)
_validate_integer("n_quadrature", n_quadrature, minimum=2)
if ti_bins is None:
ti_bins = self._ti_bins(u_min, u_max)
bins = _validate_bins("ti_bins", ti_bins)
if require_full_range:
_validate_bins_cover_range("ti_bins", bins, *self._ti_span(u_min, u_max))
u = np.linspace(u_min, u_max, int(n_quadrature))
masses = self.probability * self._in_range_masses(u_min, u_max)
conditional = np.empty((self.n_sectors, len(u), len(bins) - 1))
for surface, covered in self._sector_surfaces():
# A sector scaled by s reaches a bin edge where the model reaches
# edge / s.
scales = np.unique(self._turbulence_scale[covered])
cdf = self._turbulence_summary(
"intensity_cdf",
u,
(bins[None, :] / scales[:, None]).ravel(),
surface=surface,
).reshape(len(u), len(scales), len(bins))
for position, value in enumerate(scales):
rows = covered & (self._turbulence_scale == value)
conditional[rows] = np.diff(cdf[:, position], axis=1)
# A fixed grid misses mass where the density falls off within a few
# nodes, so each sector's TI law is averaged with weights normalized
# on the nodes and then carries its exact share of time.
mean = _mean_over_speed(self._speed, u, conditional)
return np.where(masses[:, None] > 0.0, masses[:, None] * mean, 0.0)
def _infer_sector_ti(self) -> np.ndarray:
"""Return each sector's mean TI over 4 to 25 m/s from the unscaled model."""
# The conditional mean is smooth in the speed, so 200 nodes resolve it,
# and normalizing on the same nodes cancels the quadrature error.
u = np.linspace(U_MIN, U_MAX, 200)
return _mean_over_speed(self._speed, u, self._sector_mean_intensity(u))
def _sector_mean_intensity(self, wind_speed: np.ndarray) -> np.ndarray:
"""Return the unscaled mean TI, one row per sector and column per speed."""
table = np.empty((self.n_sectors, len(wind_speed)))
for surface, rows in self._sector_surfaces():
table[rows] = self._turbulence_summary(
"mean_intensity", wind_speed, surface=surface
)
return table
def save(self, path: str | Path, *, name: str | None = None) -> Path:
"""Write this resource to *path* as a windIO ``energy_resource`` document.
As :func:`oowl.save`, it returns the written path, and :func:`oowl.load`
reads it back.
"""
from ._windio import _save_resource
return _save_resource(self, Path(path), name)
def _to_windio_energy_resource(self, name: str = "oowl wind resource") -> dict:
"""Return a windIO ``energy_resource`` mapping.
``turbulence_intensity`` is the mean over the stability climate at each
direction and wind speed. Stability is written as ``z0``,
``reference_height`` and ``LMO``, the
median Monin-Obukhov length at each wind speed. A Charnock ``z0`` is
written per wind speed and loads back as the same Charnock model, and a
coverage roughness is written as its median. windIO
has no field for the stability climate, so :func:`oowl.load` restores
the default :class:`~oowl.resources.KellyGryning`.
"""
stability_speed = self._export_speeds()
per_sector = isinstance(self.roughness_length, tuple)
# One row per sector when the surface differs by sector, else one row.
sectors = np.arange(self.n_sectors if per_sector else 1)[:, None]
speeds = np.broadcast_to(stability_speed, (len(sectors), len(stability_speed)))
roughness = self._roughness_at(speeds, sectors if per_sector else None)
inverse_length = self.stability.sample(
speeds,
np.full(speeds.shape, 0.5),
height=self.reference_height,
roughness_length=roughness,
)
# Exact neutral stability has infinite Monin-Obukhov length. windIO
# documents are JSON-compatible, so the largest finite float stands for
# it, whose reciprocal is effectively zero.
neutral_lmo = np.finfo(float).max
monin_obukhov_length = np.full_like(inverse_length, neutral_lmo)
with np.errstate(over="ignore"):
np.divide(
1.0,
inverse_length,
out=monin_obukhov_length,
where=inverse_length != 0.0,
)
np.clip(
monin_obukhov_length,
-neutral_lmo,
neutral_lmo,
out=monin_obukhov_length,
)
return {
"name": name,
"wind_resource": {
**self._speed_record(stability_speed),
"turbulence_intensity": {
"data": (
self._turbulence_scale[:, None]
* self._sector_mean_intensity(stability_speed)
).tolist(),
"dims": ["wind_direction", "wind_speed"],
},
"wind_direction": self.theta.astype(float).tolist(),
"wind_speed": stability_speed.tolist(),
"LMO": (
{
"data": monin_obukhov_length.tolist(),
"dims": ["wind_direction", "wind_speed"],
}
if per_sector
else {
"data": monin_obukhov_length[0].tolist(),
"dims": ["wind_speed"],
}
),
"z0": (
{
"data": roughness.tolist(),
"dims": ["wind_direction", "wind_speed"],
}
if per_sector
else {"data": roughness[0].tolist(), "dims": ["wind_speed"]}
if isinstance(self.roughness_length, Charnock)
else {"data": float(roughness[0, 0]), "dims": []}
),
"reference_height": self.reference_height,
},
}
def _export_speeds(self) -> np.ndarray:
"""Return the wind speeds windIO fields over speed are written at."""
return np.linspace(U_MIN, U_MAX, 22)
def _surface_metadata(self) -> dict:
"""Return the reference height and roughness model as JSON."""
def record(surface):
return (
surface.metadata()
if isinstance(surface, (Charnock, Roughness))
else surface
)
roughness = self.roughness_length
return {
"reference_height": self.reference_height,
"roughness_length": (
[record(surface) for surface in roughness]
if isinstance(roughness, tuple)
else record(roughness)
),
}
[docs]
@dataclass(repr=False, eq=False, kw_only=True)
class WindResource(_Resource):
"""Wind resource of a site, its directions, speeds, surface and turbulence.
``direction`` shares the wind between ``n_sectors`` direction sectors,
:class:`UniformDirection` by default, or :class:`SingleDirection`,
:class:`VonMises` or :class:`Elliptical`. Each sector's speed is Weibull,
and roses follow analytically. ``mean_wind_speed`` is the exact
probability-weighted mean in m/s, 9 by default. ``weibull_a`` instead gives
the Weibull scale in m/s, one value or one per sector, and the mean follows
from it. ``weibull_k`` is one shape for all sectors or one per sector.
Each case is drawn along a chain of links. ``stability`` draws ``1/L``
from a :class:`~oowl.resources.KellyGryning` climate, ``turbulence`` draws
``sigma_u`` from :class:`~oowl.resources.SurfaceLayer` after it and also
gives the TI rose, and ``mann`` adds the :class:`~oowl.resources.Mann`
parameters. ``reference_height`` and ``roughness_length`` in meters define
the surface, and with the stability the log-law shear of each case. The
defaults describe open water. A :class:`~oowl.resources.Charnock`
roughness grows with the wind instead, and a sequence gives each sector
its own upwind surface, a length or a Charnock model, so a coast has
offshore and onshore sectors whose turbulence, shear and stability differ.
"""
direction: UniformDirection | SingleDirection | VonMises | Elliptical | None = (
field(default_factory=lambda: UniformDirection())
)
n_sectors: int = 12
mean_wind_speed: float | None = None
weibull_k: float | Sequence[float] | np.ndarray | None = 2.0
weibull_a: float | Sequence[float] | np.ndarray | None = None
stability: KellyGryning = field(default_factory=KellyGryning)
turbulence: SurfaceLayer = field(default_factory=SurfaceLayer)
reference_height: float = 100.0
roughness_length: float | Charnock | Sequence[float | Charnock] = 0.0002
mann: Mann | None = None
@property
def _weibull_k(self) -> np.ndarray:
"""Per-sector Weibull shape as an array."""
return np.asarray(self.weibull_k, dtype=float)
@property
def _weibull_a(self) -> np.ndarray:
"""Per-sector Weibull scale as an array."""
return np.asarray(self.weibull_a, dtype=float)
def __post_init__(self):
# A resource from data, such as GWA, carries its own directions.
if self.direction is not None and not isinstance(self.direction, _Direction):
raise TypeError(
"direction must be UniformDirection, SingleDirection, VonMises "
"or Elliptical"
)
self.n_sectors = _validate_integer("n_sectors", self.n_sectors)
self.theta = self._sector_angles(self.n_sectors)
self.probability = np.asarray(
self._sector_probabilities(self.n_sectors), dtype=float
)
self._validate_probability()
self.weibull_k = _per_sector("weibull_k", self.weibull_k, self.n_sectors)
# The speed climate is set by its mean or by the Weibull scale.
self._given_weibull_a = None
if self.weibull_a is not None:
if self.mean_wind_speed is not None:
raise ValueError("give mean_wind_speed or weibull_a, not both")
self._given_weibull_a = _per_sector(
"weibull_a", self.weibull_a, self.n_sectors
)
self.mean_wind_speed = float(
np.dot(
self.probability,
_WeibullSpeed(self._given_weibull_a, self._weibull_k).mean(),
)
)
else:
self.mean_wind_speed = _validate_positive(
"mean_wind_speed",
9.0 if self.mean_wind_speed is None else self.mean_wind_speed,
)
sector_mean_speed = np.asarray(self._sector_mean_speed(), dtype=float)
if sector_mean_speed.shape != (self.n_sectors,):
raise ValueError("sector mean speed must have length n_sectors")
if not np.all(np.isfinite(sector_mean_speed)) or np.any(sector_mean_speed <= 0):
raise ValueError("sector mean speed must be positive and finite")
weighted_mean = float(np.dot(self.probability, sector_mean_speed))
if not np.isclose(weighted_mean, self.mean_wind_speed, rtol=1e-10):
raise ValueError(
"probability-weighted sector means must equal mean_wind_speed"
)
self._sector_mean_speeds = sector_mean_speed
from scipy.special import gamma
self.weibull_a = (
self._given_weibull_a.copy()
if self._given_weibull_a is not None
else sector_mean_speed / gamma(1.0 + 1.0 / self._weibull_k)
)
self._speed = self._speed_distribution()
from ._site.stability import _resolve_model, _resolve_roughness
from ._site.turbulence import _resolve_turbulence
given = self.roughness_length
if isinstance(given, (list, tuple, np.ndarray)):
surfaces = tuple(_resolve_roughness(surface) for surface in list(given))
if len(surfaces) != self.n_sectors:
raise ValueError("roughness_length needs one surface per sector")
self.roughness_length = surfaces
else:
self.roughness_length = _resolve_roughness(self.roughness_length)
self.stability = _resolve_model(self.stability)
self.turbulence = _resolve_turbulence(self.turbulence)
if self.mann is not None and not isinstance(self.mann, Mann):
raise TypeError("mann must be None or a Mann model")
self._setup_chain()
@property
def _pinned_seed(self) -> bool:
"""Whether the direction model draws nothing or was seeded."""
return getattr(self.direction, "_pinned_seed", True)
def _reseeded(self, seed: int) -> "WindResource":
"""Return this resource, or one whose unseeded direction is redrawn."""
if self._pinned_seed:
return self
direction = cast(VonMises, self.direction)
given_a = self._given_weibull_a
return replace(
self,
direction=direction._reseeded(seed),
mean_wind_speed=None if given_a is not None else self.mean_wind_speed,
weibull_a=given_a,
)
def __repr__(self) -> str:
"""Return the call that rebuilds the resource, every argument shown.
Only arguments the class accepts are shown, so a resource that derives
its speeds, such as :class:`GWA`, does not print them.
"""
init = type(self).__init__
if init is not WindResource.__init__:
names: list[str] = []
for name, parameter in inspect.signature(init).parameters.items():
if parameter.kind is inspect.Parameter.VAR_KEYWORD:
names += [item.name for item in fields(WindResource) if item.init]
elif name != "self":
names.append(name)
else:
names = [item.name for item in fields(WindResource) if item.init]
parts = []
for name in dict.fromkeys(names):
if not hasattr(self, name):
continue
value = getattr(self, name)
# The speed climate was given as its mean or as the Weibull scale,
# and only that one is an argument.
if name == "mean_wind_speed" and self._given_weibull_a is not None:
value = None
if name == "weibull_a":
value = self._given_weibull_a
# Per-sector values in full, so the text rebuilds the resource.
if isinstance(value, np.ndarray):
value = float(value[0]) if np.all(value == value[0]) else tuple(value)
parts.append(f"{name}={_describe(value)}")
return f"{type(self).__name__}({', '.join(parts)})"
def _sector_angles(self, n_sectors: int) -> np.ndarray:
if self.direction is None:
return sector_angles(n_sectors)
return self.direction.sector_angles(n_sectors)
def _sector_probabilities(self, n_sectors: int) -> np.ndarray:
if self.direction is None:
raise ValueError("this resource carries no direction model")
return self.direction.probabilities(n_sectors)
def _sector_mean_speed(self) -> np.ndarray:
if self._given_weibull_a is not None:
return _WeibullSpeed(self._given_weibull_a, self._weibull_k).mean()
return np.full(self.n_sectors, cast(float, self.mean_wind_speed))
def _speed_distribution(self) -> _WeibullSpeed | _TabulatedSpeed:
"""Return the distribution of each sector's speed, set before the chain."""
return _WeibullSpeed(self._weibull_a, self._weibull_k)
@property
def sector_mean_wind_speed(self) -> np.ndarray:
"""Return the conditional mean wind speed in every sector, in m/s."""
return self._sector_mean_speeds.copy()
@property
def wind_speed_standard_deviation(self) -> float:
"""Return the standard deviation in m/s of the speed over all time."""
second_moment = self._speed_distribution().second_moment()
variance = float(np.dot(self.probability, second_moment))
variance -= cast(float, self.mean_wind_speed) ** 2
return float(np.sqrt(max(variance, 0.0)))
def _speed_record(self, wind_speed: np.ndarray) -> dict:
"""Return the windIO fields of the direction and speed distribution."""
return {
"sector_probability": {
"data": self.probability.astype(float).tolist(),
"dims": ["wind_direction"],
},
"weibull_a": {
"data": self._weibull_a.astype(float).tolist(),
"dims": ["wind_direction"],
},
"weibull_k": {
"data": self._weibull_k.astype(float).tolist(),
"dims": ["wind_direction"],
},
}
class _WindIOResource(WindResource):
"""A resource read from windIO data, Weibull, a probability table or a rose.
``speed`` is the distribution of each sector. ``turbulence`` is the file's
``turbulence_intensity`` field, or ``None``, ``wind_speed`` its speed
coordinate and ``table`` a loaded probability field, kept for saving back.
"""
def __init__(
self,
theta,
probability,
speed,
turbulence,
wind_speed,
stability,
reference_height,
roughness_length,
table=None,
):
self._windio_turbulence = (turbulence, wind_speed)
self._windio_table = table
self._windio_theta = np.asarray(theta, dtype=float)
self._windio_probability = np.asarray(probability, dtype=float)
n_sectors = len(self._windio_theta)
if n_sectors < 1 or self._windio_probability.shape != (n_sectors,):
raise ValueError("windIO resource arrays must have equal non-zero length")
if not np.all(np.isfinite(self._windio_theta)):
raise ValueError("windIO wind directions must be finite")
if not np.all(np.isfinite(self._windio_probability)) or np.any(
self._windio_probability < 0
):
raise ValueError(
"windIO sector probabilities must be finite and non-negative"
)
probability_sum = float(self._windio_probability.sum())
if probability_sum <= 0:
raise ValueError("windIO sector probabilities must contain positive mass")
weibull = isinstance(speed, _WeibullSpeed)
if weibull and (
speed.scale.shape != (n_sectors,)
or not np.all(np.isfinite(speed.scale))
or np.any(speed.scale <= 0)
):
raise ValueError("windIO Weibull A must be positive and finite")
if weibull and (
speed.shape.shape != (n_sectors,)
or not np.all(np.isfinite(speed.shape))
or np.any(speed.shape <= 0)
):
raise ValueError("windIO Weibull k must be positive and finite")
self._windio_speed = speed
mean = float(np.dot(self._windio_probability / probability_sum, speed.mean()))
super().__init__(
direction=None,
n_sectors=n_sectors,
mean_wind_speed=mean,
weibull_k=speed.shape if weibull else 2.0,
turbulence=SurfaceLayer(
intensity=None
if turbulence is None
else np.asarray(
_windio_turbulence_intensity(turbulence, wind_speed, speed)
).tolist()
),
stability=stability,
reference_height=reference_height,
roughness_length=roughness_length,
)
self.theta = self._windio_theta.copy()
# The file's own scale, as rebuilding it from the mean would round it.
self.weibull_a = speed.scale.copy() if weibull else None
if not weibull:
self.weibull_k = None
def __repr__(self) -> str:
# Loaded data has no constructor call, so summarize what it holds.
return (
f"windIO resource(n_sectors={self.n_sectors}, "
f"mean_wind_speed={self.mean_wind_speed:.6g})"
)
def _sector_probabilities(self, n_sectors: int) -> np.ndarray:
return self._windio_probability
def _sector_mean_speed(self) -> np.ndarray:
return self._windio_speed.mean()
def _speed_distribution(self) -> _WeibullSpeed | _TabulatedSpeed:
# Each sector's turbulence level averages over its speeds, so the
# file's distribution must be in place before the chain is set up.
return self._windio_speed
def _export_speeds(self) -> np.ndarray:
if isinstance(self._windio_speed, _TabulatedSpeed):
return self._windio_speed.speeds.copy()
return super()._export_speeds()
def _speed_record(self, wind_speed: np.ndarray) -> dict:
if not isinstance(self._windio_speed, _TabulatedSpeed):
return super()._speed_record(wind_speed)
return {
"probability": {
"data": (
self.probability[:, None] * self._windio_speed.masses
).tolist(),
"dims": ["wind_direction", "wind_speed"],
}
}
def metadata(self) -> dict:
"""Return a record of the loaded resource, type ``windIO``."""
speed = self._windio_speed
record = {
"type": "windIO",
"wind_direction": self.theta.astype(float).tolist(),
"n_sectors": self.n_sectors,
"mean_wind_speed": cast(float, self.mean_wind_speed),
"sector_probability": self.probability.astype(float).tolist(),
"turbulence": self.turbulence.metadata(),
"surface": self._surface_metadata(),
"stability": self.stability.metadata(),
}
if isinstance(speed, _TabulatedSpeed):
record["wind_speed"] = speed.speeds.tolist()
record["speed_probability"] = speed.masses.tolist()
else:
record["weibull_a"] = self._weibull_a.astype(float).tolist()
record["weibull_k"] = self._weibull_k.astype(float).tolist()
return record
def _to_windio_energy_resource(self, name: str = "oowl wind resource") -> dict:
"""Return the windIO mapping, with the turbulence intensity as loaded.
The loaded field is kept unless it runs over other wind speeds than
the ones written, so a loaded file saves back unchanged.
"""
import copy
document = super()._to_windio_energy_resource(name)
if self._windio_table is not None:
document["wind_resource"]["probability"] = copy.deepcopy(self._windio_table)
field, speed = self._windio_turbulence
dims = field.get("dims", []) if isinstance(field, dict) else []
written = document["wind_resource"]
if field is not None and (
"wind_speed" not in dims
or np.array_equal(_windio_data(speed), written["wind_speed"])
):
written["turbulence_intensity"] = copy.deepcopy(field)
return document
def _windio_data(value):
"""Return windIO data, unwrapping a ``{"data": ...}`` wrapper."""
return value.get("data") if isinstance(value, dict) else value
# Dimensions of the windIO forms OOWL does not read, a resource mapped over
# the site, given per turbine or varying with height.
_UNSUPPORTED_WINDIO_DIMENSIONS = ("x", "y", "wind_turbine", "height")
def _from_windio_energy_resource(energy_resource):
"""Return the resource a windIO ``energy_resource`` mapping describes.
Per-sector Weibull parameters and a probability table over wind direction
and speed load as the distributions they state, each tabulated speed
carrying the cell around it. A table over wind direction alone with one
wind speed is a rose whose every case blows at that speed. A time series
is counted into 12 sectors and 1 m/s speed bins over the speeds observed,
and loads as that table. ``None`` stands for a missing resource and for
OOWL's own placeholder, a probability without dimensions. A resource over
``x`` and ``y``, per turbine or per height raises :class:`ValueError`.
"""
if energy_resource is None:
return None
if not isinstance(energy_resource, dict) or not isinstance(
energy_resource.get("wind_resource"), dict
):
raise TypeError("energy_resource must be a mapping holding a wind_resource")
wind_resource = energy_resource["wind_resource"]
probability_field = wind_resource.get("probability")
# A site saved without a resource carries a probability of one over no
# dimensions, which describes no wind.
scalar = np.ndim(_windio_data(probability_field)) == 0
if isinstance(probability_field, dict) and scalar:
return None
unsupported = [
name for name in _UNSUPPORTED_WINDIO_DIMENSIONS if name in wind_resource
]
if unsupported:
raise ValueError(
f"windIO wind_resource over {unsupported} is not supported, OOWL reads "
"one resource for the whole site at its reference_height"
)
if "wind_direction" not in wind_resource:
raise ValueError("windIO wind_resource is missing wind_direction")
def data(name):
return _windio_data(wind_resource[name])
table = None
wind_speed = wind_resource.get("wind_speed")
if wind_speed is not None and "time" not in wind_resource:
wind_speed = np.atleast_1d(np.asarray(_windio_data(wind_speed), dtype=float))
if "time" in wind_resource:
if wind_speed is None:
raise ValueError("windIO time series is missing wind_speed")
theta, probability, wind_speed, values, slowest = _windio_time_series(
data("time"), data("wind_direction"), _windio_data(wind_speed)
)
speed = _TabulatedSpeed(wind_speed, values, lowest=slowest)
elif all(
name in wind_resource
for name in ("sector_probability", "weibull_a", "weibull_k")
):
theta = data("wind_direction")
probability = data("sector_probability")
speed = _WeibullSpeed(data("weibull_a"), data("weibull_k"))
elif isinstance(probability_field, dict) and wind_speed is not None:
theta = data("wind_direction")
table = probability_field
values = np.asarray(table.get("data"), dtype=float)
dims = list(table.get("dims", []))
if dims == ["wind_direction"] and len(wind_speed) == 1:
values = values[:, None]
elif dims == ["wind_speed", "wind_direction"]:
values = values.T
elif dims != ["wind_direction", "wind_speed"]:
raise ValueError(f"windIO probability over {dims} is not supported")
probability = values.sum(axis=1)
speed = _TabulatedSpeed(wind_speed, values)
else:
raise ValueError(
"windIO wind_resource needs sector Weibull parameters, a probability "
"table over wind_direction and wind_speed, or a time series"
)
# windIO carries no stability climate, so the default model is restored.
stability = KellyGryning()
reference_height = float(wind_resource.get("reference_height", 100.0))
roughness_length = (
_windio_roughness(wind_resource["z0"], wind_speed, reference_height)
if "z0" in wind_resource
else 0.0002
)
return _WindIOResource(
theta,
probability,
speed,
wind_resource.get("turbulence_intensity"),
wind_speed,
stability,
reference_height,
roughness_length,
table=table,
)
def _windio_time_series(time, wind_direction, wind_speed):
"""Count a windIO time series into a rose and a table over speed.
Returns the sector angles, each sector's share of the steps, the table
speeds, the counts per sector and speed and the slowest step. The steps
fall into 12 sectors of 30 degrees and into 1 m/s bins centered on the
whole m/s from the one at or below the slowest step to the one at or above
the fastest. The table's last cell reaches only to its whole m/s, where no
step lies beyond. The table starts at 1 m/s at the lowest, as the
turbulence models have no 0 m/s, and calm steps below 1 m/s count to that
cell, which reaches down to the slowest step.
"""
n_sectors = 12
direction = np.asarray(wind_direction, dtype=float)
speed = np.asarray(wind_speed, dtype=float)
if direction.ndim != 1 or direction.shape != speed.shape:
raise ValueError(
"windIO time series wind_direction and wind_speed must be one value "
"per time"
)
if len(np.atleast_1d(time)) != len(speed):
raise ValueError("windIO time series time must have one entry per step")
if not np.all(np.isfinite(direction)) or not np.all(np.isfinite(speed)):
raise ValueError(
"windIO time series wind_direction and wind_speed must be finite"
)
if np.any(speed < 0):
raise ValueError("windIO time series wind_speed must be non-negative")
sector = np.round(direction / (360.0 / n_sectors)).astype(int) % n_sectors
lowest = max(np.floor(speed.min()), 1.0)
speeds = np.arange(lowest, max(np.ceil(speed.max()), lowest) + 1.0)
edges = np.append(speeds - 0.5, speeds[-1] + 0.5)
edges[0] = 0.0
counts, _, _ = np.histogram2d(
sector, speed, bins=[np.arange(n_sectors + 1) - 0.5, edges]
)
probability = counts.sum(axis=1) / len(speed)
return sector_angles(n_sectors), probability, speeds, counts, float(speed.min())
def _windio_turbulence_intensity(field, wind_speed, distribution):
"""Return a windIO turbulence intensity as one value or one per sector.
A turbulence intensity over wind speed is averaged over each sector's
speeds between 4 and 25 m/s, the mean a resource's target sets.
"""
if not isinstance(field, dict):
return field
values = np.asarray(field.get("data"), dtype=float)
dims = list(field.get("dims", []))
if "wind_speed" not in dims:
if dims not in ([], ["wind_direction"]):
raise ValueError(f"turbulence_intensity over {dims} is not supported")
return values
if dims == ["wind_speed"]:
values = values[None, :]
elif dims == ["wind_speed", "wind_direction"]:
values = values.T
elif dims != ["wind_direction", "wind_speed"]:
raise ValueError(f"turbulence_intensity over {dims} is not supported")
speed = np.asarray(_windio_data(wind_speed), dtype=float)
if speed.shape != values.shape[-1:]:
raise ValueError("turbulence_intensity needs one wind_speed per column")
grid = np.linspace(U_MIN, U_MAX, 200)
curves = np.array([np.interp(grid, speed, row) for row in values])
return _mean_over_speed(distribution, grid, curves)
def _windio_roughness(z0, wind_speed, reference_height):
"""Return a windIO ``z0`` as a length, a Charnock model, or one per sector.
A ``z0`` that changes with wind speed is inverted for
``A_c = g z0 / u_star^2`` with the neutral ``u_star`` of that speed, the
Charnock constant it was written with. A ``z0`` over wind direction gives
each sector its own surface.
"""
from ._site.stability import _GRAVITY
if not isinstance(z0, dict):
return float(z0)
values = np.asarray(z0.get("data"), dtype=float)
dims = list(z0.get("dims", []))
if dims == ["wind_speed", "wind_direction"]:
values, dims = values.T, ["wind_direction", "wind_speed"]
if dims not in (
[],
["wind_speed"],
["wind_direction"],
["wind_direction", "wind_speed"],
):
raise ValueError(f"z0 over {dims} is not supported")
if "wind_speed" in dims and wind_speed is None:
raise ValueError("z0 over wind_speed needs a wind_speed coordinate")
def surface(row):
row = np.atleast_1d(row)
if np.ptp(row) == 0.0:
return float(row[0])
u_star = friction_velocity(wind_speed, reference_height, row)
return Charnock(constant=float(np.mean(_GRAVITY * row / u_star**2)))
if "wind_direction" not in dims:
return surface(values)
return tuple(surface(row) for row in values)
class _Direction(_FieldRepr):
"""How a resource shares the wind between its direction sectors.
``sector_angles`` gives the sector centers and ``probabilities`` each
sector's share of time.
"""
def sector_angles(self, n_sectors: int) -> np.ndarray:
"""Return the sector centers in degrees, ``360 i / n`` by default."""
return sector_angles(n_sectors)
def probabilities(self, n_sectors: int) -> np.ndarray:
"""Return each sector's share of time, summing to one."""
raise NotImplementedError
def metadata(self) -> dict:
"""Return a JSON-compatible record of the direction model."""
from ._sampling_specs import _json_safe
return {
"type": type(self).__name__,
**{
item.name: _json_safe(getattr(self, item.name))
for item in fields(cast(Any, self))
},
}
[docs]
@dataclass(frozen=True, repr=False)
class SingleDirection(_Direction):
"""All probability in one sector, centered on ``direction``.
``direction`` is in degrees, 0 north and 90 east. The sectors turn so one
is centered on it, and every case blows from exactly that direction.
"""
#: Direction the wind blows from, in degrees.
direction: float
def __post_init__(self):
direction = _validate_finite("direction", self.direction)
object.__setattr__(self, "direction", direction % 360.0)
[docs]
def sector_angles(self, n_sectors: int) -> np.ndarray:
return sector_angles(n_sectors) + self.direction % (360.0 / n_sectors)
[docs]
def probabilities(self, n_sectors: int) -> np.ndarray:
n_sectors = _validate_integer("n_sectors", n_sectors)
offset = self.direction % (360.0 / n_sectors)
index = int(round((self.direction - offset) * n_sectors / 360.0)) % n_sectors
probability = np.zeros(n_sectors, dtype=float)
probability[index] = 1.0
return probability
[docs]
@dataclass(frozen=True, repr=False)
class VonMises(_Direction):
"""A mixture of von Mises distributions on a uniform background.
The wind direction model of Carta et al. (2008),
https://doi.org/10.1016/j.enconman.2007.10.017. Each mode has a mean
direction in degrees, a concentration, larger for a narrower mode, and a
weight. ``n_modes`` modes are drawn from ``seed`` unless ``modes`` gives
them as ``(direction, concentration, weight)`` tuples. Without a seed, a
sampler redraws the modes for each scenario from its own seed. A drawn
mode takes its direction uniformly, its concentration from
``concentration_range`` and its weight from 0.5 to 1.5. The weights are
scaled to share ``1 - background_fraction``, and the background spreads
the rest evenly, a mode of zero concentration.
"""
#: Number of modes drawn from ``seed``.
n_modes: int = 2
_: KW_ONLY
#: Modes as ``(direction, concentration, weight)``, replacing drawn ones.
modes: Sequence[tuple[float, float, float]] | None = None
#: Share of probability spread evenly over the compass.
background_fraction: float = 0.20
#: Range of the drawn concentrations.
concentration_range: tuple[float, float] = (3.0, 12.0)
#: Seed of the drawn modes.
seed: int | None = None
# Set per instance in __post_init__. A ClassVar stays out of the fields.
_pinned_seed: ClassVar[bool] = True
def __post_init__(self):
background = _validate_unit_interval(
"background_fraction", self.background_fraction
)
concentration = _validate_finite_pair(
"concentration_range", self.concentration_range
)
if concentration[0] <= 0 or concentration[1] < concentration[0]:
raise ValueError("concentration_range must be positive and increasing")
# Unseeded drawn modes are redrawn per scenario, see _reseeded.
object.__setattr__(
self, "_pinned_seed", self.seed is not None or self.modes is not None
)
for name, value in (
("background_fraction", background),
("concentration_range", concentration),
# Given modes never draw, so an unset seed stays None and the repr
# and metadata of equal models match.
(
"seed",
None
if self.modes is not None and self.seed is None
else resolve_seed(self.seed),
),
):
object.__setattr__(self, name, value)
if self.modes is None:
count = _validate_integer("n_modes", self.n_modes)
object.__setattr__(self, "n_modes", count)
rng = as_rng(self.seed)
modes = tuple(
(
float(rng.uniform(0.0, 360.0)),
float(rng.uniform(*concentration)),
float(rng.uniform(0.5, 1.5)),
)
for _ in range(int(count))
)
else:
parsed = []
for item in self.modes:
if len(item) != 3:
raise ValueError(
"modes must contain (direction, concentration, weight) tuples"
)
direction, kappa, weight = (
_validate_real("mode " + name, value)
for name, value in zip(
("direction", "concentration", "weight"), item
)
)
if not np.isfinite(direction):
raise ValueError("mode direction must be finite")
if not np.isfinite(weight) or weight < 0:
raise ValueError("mode weight must be finite and >= 0")
parsed.append(
(
direction % 360.0,
_validate_positive("mode concentration", kappa),
weight,
)
)
if not parsed:
raise ValueError("modes must contain at least one tuple")
modes = tuple(parsed)
object.__setattr__(self, "modes", modes)
object.__setattr__(self, "n_modes", len(modes))
weights = np.asarray([mode[2] for mode in modes], dtype=float)
if weights.sum() <= 0:
raise ValueError("at least one mode weight must be positive")
object.__setattr__(self, "_modes", modes)
object.__setattr__(
self, "_weights", (1.0 - background) * weights / weights.sum()
)
def _reseeded(self, seed: int) -> "VonMises":
"""Return this model, or an unseeded one redrawn from ``seed``."""
return self if self._pinned_seed else replace(self, seed=int(seed))
def _density(self, theta_rad) -> np.ndarray:
from scipy.special import i0e
theta_rad = np.asarray(theta_rad, dtype=float)
density = np.full(theta_rad.shape, self.background_fraction / (2.0 * np.pi))
state = cast(Any, self)
for (direction, kappa, _), weight in zip(state._modes, state._weights):
mu = np.deg2rad(direction)
density += (
weight
* np.exp(kappa * (np.cos(theta_rad - mu) - 1.0))
/ (2.0 * np.pi * i0e(kappa))
)
return density
[docs]
def probabilities(self, n_sectors: int) -> np.ndarray:
n_sectors = _validate_integer("n_sectors", n_sectors)
delta = 2.0 * np.pi / n_sectors
n_points = max(64, int(np.ceil(4096 / n_sectors)))
offsets = (np.arange(n_points, dtype=float) + 0.5) / n_points
# One row of midpoints per sector.
theta = np.arange(n_sectors)[:, None] * delta - 0.5 * delta + offsets * delta
probability = self._density(theta).mean(axis=1) * delta
return probability / probability.sum()
def _elliptical_cdf(theta_rad, a):
"""CDF of the unit-area ellipse direction distribution."""
theta = np.asarray(theta_rad, dtype=float)
theta_mod = theta % (2.0 * np.pi)
quadrant = np.clip(np.floor(theta_mod / (0.5 * np.pi)).astype(int), 0, 3)
q1 = np.where(
quadrant == 0,
theta_mod,
np.where(
quadrant == 1,
np.pi - theta_mod,
np.where(quadrant == 2, theta_mod - np.pi, 2.0 * np.pi - theta_mod),
),
)
q1 = np.clip(q1, 0.0, 0.5 * np.pi)
at_boundary = q1 >= 0.5 * np.pi - 1e-14
safe_q1 = np.where(at_boundary, 0.0, q1)
first_quadrant = np.where(
at_boundary,
0.25,
np.arctan(np.pi * a**2 * np.tan(safe_q1)) / (2.0 * np.pi),
)
return np.where(
quadrant == 0,
first_quadrant,
np.where(
quadrant == 1,
0.5 - first_quadrant,
np.where(quadrant == 2, 0.5 + first_quadrant, 1.0 - first_quadrant),
),
)
def _periodic_elliptical_cdf(theta_rad, a):
"""Unwrapped CDF of the periodic elliptical direction distribution."""
theta = np.asarray(theta_rad, dtype=float)
return _elliptical_cdf(theta, a) + np.floor(theta / (2.0 * np.pi))
def _elliptical_front_cdf(theta_rad, a):
"""Unwrapped integral of elliptical density restricted to the front half."""
theta = np.asarray(theta_rad, dtype=float)
period = 2.0 * np.pi
cycle = np.floor((theta + 0.5 * np.pi) / period)
front_lower = -0.5 * np.pi + cycle * period
front_upper = front_lower + np.pi
within_cycle = _periodic_elliptical_cdf(
np.minimum(theta, front_upper), a
) - _periodic_elliptical_cdf(front_lower, a)
return 0.5 * cycle + within_cycle
[docs]
@dataclass(frozen=True, repr=False, kw_only=True)
class Elliptical(_Direction):
"""Hart's three-parameter elliptical wind rose.
Described in https://doi.org/10.5194/wes-10-1821-2025. Directions follow
an ellipse of unit area whose ``semi_axis`` lies along the prevailing
axis. The default ``1 / sqrt(pi)`` is a circle, uniform when ``folding``
is 0, a longer one concentrates the rose along the axis and a shorter one
across it. ``folding`` moves probability from the opposite half towards
the axis. ``direction`` is the prevailing direction in degrees, 0 north
and 90 east.
"""
#: Semi-axis of the unit-area ellipse along the prevailing axis.
semi_axis: float = float(1.0 / np.sqrt(np.pi))
#: Share moved from the opposite half, in [0, 1].
folding: float = 0.0
#: Prevailing direction in degrees.
direction: float = 0.0
def __post_init__(self):
semi_axis = _validate_positive("semi_axis", self.semi_axis)
object.__setattr__(self, "semi_axis", semi_axis)
folding = _validate_unit_interval("folding", self.folding)
direction = _validate_finite("direction", self.direction)
object.__setattr__(self, "folding", folding)
object.__setattr__(self, "direction", direction % 360.0)
[docs]
def probabilities(self, n_sectors: int) -> np.ndarray:
n_sectors = _validate_integer("n_sectors", n_sectors)
width = 2.0 * np.pi / n_sectors
theta = np.deg2rad(sector_angles(n_sectors)) - np.deg2rad(self.direction)
lower, upper = theta - 0.5 * width, theta + 0.5 * width
base_mass = _periodic_elliptical_cdf(upper, self.semi_axis) - (
_periodic_elliptical_cdf(lower, self.semi_axis)
)
front_mass = _elliptical_front_cdf(
upper, self.semi_axis
) - _elliptical_front_cdf(lower, self.semi_axis)
probabilities = (1.0 - self.folding) * base_mass
probabilities += 2.0 * self.folding * front_mass
total = probabilities.sum()
if not np.isfinite(total) or total <= 0:
raise ValueError("elliptical direction model produced no probability")
return probabilities / total
# Roughness WAsP gives the water class, listed as 0 in libfiles.
_WATER = 0.0002
[docs]
class GWA(WindResource):
"""Wind resource from the Global Wind Atlas (GWA) Global Wind Climate (GWC).
Downloads the GWC libfile for a ``latitude``/``longitude`` point from the
Global Wind Atlas v4 API and interpolates its per-sector Weibull ``A`` and
``k`` linearly in ``ln z0`` to ``roughness_length`` and linearly in height
to ``reference_height``. The sector frequencies, which the libfile gives
per roughness class only, are interpolated in ``ln z0``. The same surface
values define the stability model. Adapted from PyWake's
``py_wake.site.xrsite.GlobalWindAtlasSite``.
One roughness length per atlas sector gives each sector the atlas climate
at its own roughness, for a site whose upwind surface changes with
direction, such as a coast. The frequencies are then renormalized. The
atlas climates assume that roughness holds far upwind in each sector.
"""
#: Global Wind Atlas v4 GWC libfile endpoint.
GWA_URL = "https://globalwindatlas.info/api/gwa/custom/Lib"
#: Browser ``Referer`` required by the v4 API.
REQUEST_HEADERS = {
"User-Agent": "oowl.GWA",
"Referer": "https://globalwindatlas.info/",
}
#: Maximum response body size accepted before parsing.
MAX_RESPONSE_BYTES = 10 * 1024 * 1024
#: Bytes read from the response at a time.
READ_CHUNK_BYTES = 64 * 1024
def __init__(
self,
latitude: float,
longitude: float,
roughness_length: float | Sequence[float],
reference_height: float,
*,
stability: KellyGryning = KellyGryning(),
turbulence: SurfaceLayer = SurfaceLayer(),
mann: Mann | None = None,
timeout: float = 30.0,
libfile_path: str | Path | None = None,
):
"""
Parameters
----------
latitude, longitude : float
Latitude in ``[-90, 90]`` and longitude in ``[-180, 180]`` degrees.
roughness_length : float or sequence of float
Roughness length in meters, or one per atlas sector.
reference_height : float
Height above ground in meters.
stability : KellyGryning
Stability climate of the site.
turbulence : SurfaceLayer
Turbulence model.
mann : Mann, optional
Mann turbulence parameters, ``None`` to omit them.
timeout : float
Maximum seconds to wait for the Global Wind Atlas API. Ignored when
``libfile_path`` is supplied.
libfile_path : str or pathlib.Path, optional
Previously downloaded GWC libfile to read instead of the network.
"""
self.latitude = _validate_real("latitude", latitude)
self.longitude = _validate_real("longitude", longitude)
lengths = _validate_real_array("roughness_length", roughness_length)
reference_height = _validate_real("reference_height", reference_height)
if not np.isfinite(self.latitude) or not -90.0 <= self.latitude <= 90.0:
raise ValueError("latitude must be finite and in [-90, 90]")
if not np.isfinite(self.longitude) or not -180.0 <= self.longitude <= 180.0:
raise ValueError("longitude must be finite and in [-180, 180]")
if (
lengths.ndim > 1
or lengths.size == 0
or not np.all(np.isfinite(lengths))
or np.any(lengths <= 0.0)
):
raise ValueError(
"roughness_length must be positive and finite, one length or one "
"per sector"
)
if not np.isfinite(reference_height) or reference_height <= lengths.max():
raise ValueError(
"reference_height must be finite and above roughness_length"
)
self.timeout = _validate_positive("timeout", timeout)
self.libfile_path = None if libfile_path is None else Path(libfile_path)
if self.libfile_path is None:
gwc = self._read_gwc(self.latitude, self.longitude, timeout=self.timeout)
else:
gwc = self._read_gwc_file(self.libfile_path)
native = gwc["frequency"].shape[1]
if lengths.ndim == 1 and len(lengths) != native:
raise ValueError(
f"roughness_length needs one length per atlas sector, {native} here"
)
a, k, freq = self._interp_gwc(gwc, lengths, reference_height)
self._gwa_freq = freq
self._gwa_mean = _WeibullSpeed(a, k).mean()
mean_wind_speed = float(np.dot(freq, self._gwa_mean))
super().__init__(
direction=None,
n_sectors=len(freq),
mean_wind_speed=mean_wind_speed,
weibull_k=k,
stability=stability,
turbulence=turbulence,
reference_height=reference_height,
roughness_length=(
float(lengths) if lengths.ndim == 0 else tuple(lengths.tolist())
),
mann=mann,
)
@classmethod
def _read_gwc(cls, lat: float, long: float, timeout: float = 30.0) -> dict:
"""Download and parse a size-bounded GWC libfile within a total deadline."""
timeout = _validate_positive("timeout", timeout)
deadline = time.monotonic() + timeout
url = f"{cls.GWA_URL}?lat={lat}&long={long}&areaname=oowl"
request = urllib.request.Request(url, headers=cls.REQUEST_HEADERS)
result: queue.Queue[tuple[str, object]] = queue.Queue(maxsize=1)
response_holder = []
cancelled = threading.Event()
def download() -> None:
try:
with urllib.request.urlopen(
request,
context=cls._ssl_context(),
timeout=timeout,
) as response:
response_holder.append(response)
if cancelled.is_set():
return
payload = bytearray()
reader = getattr(response, "read1", response.read)
while True:
if cancelled.is_set():
return
remaining = cls.MAX_RESPONSE_BYTES + 1 - len(payload)
chunk = reader(min(cls.READ_CHUNK_BYTES, remaining))
if not chunk:
break
payload.extend(chunk)
if len(payload) > cls.MAX_RESPONSE_BYTES:
raise RuntimeError(
"Global Wind Atlas response exceeds "
f"{cls.MAX_RESPONSE_BYTES} bytes"
)
result.put(("payload", bytes(payload)))
except Exception as error:
result.put(("error", error))
worker = threading.Thread(target=download, daemon=True)
worker.start()
worker.join(max(0.0, deadline - time.monotonic()))
if worker.is_alive():
cancelled.set()
if response_holder:
try:
response_holder[0].close()
except Exception:
pass
raise RuntimeError(
"could not download Global Wind Atlas data for "
f"latitude={lat}, longitude={long}: "
f"total deadline of {timeout:g}s exceeded"
)
kind, value = result.get_nowait()
if kind == "error":
error = value if isinstance(value, BaseException) else RuntimeError(value)
raise RuntimeError(
"could not download Global Wind Atlas data for "
f"latitude={lat}, longitude={long}: {error}"
) from error
if not isinstance(value, (bytes, bytearray)):
raise RuntimeError("Global Wind Atlas download returned no payload")
try:
lines = bytes(value).decode("utf-8").strip().splitlines()
except UnicodeDecodeError as error:
raise RuntimeError(
"Global Wind Atlas returned a libfile that is not valid UTF-8"
) from error
return cls._parse_gwc(lines)
@classmethod
def _read_gwc_file(cls, path: str | Path) -> dict:
"""Read and parse a local GWC libfile without network access."""
path = Path(path)
try:
lines = path.read_text(encoding="utf-8").strip().splitlines()
except OSError as error:
raise RuntimeError(f"could not read GWC libfile {path}: {error}") from error
return cls._parse_gwc(lines)
@staticmethod
def _ssl_context():
"""Return an SSL context, using certifi's CA bundle when installed."""
import ssl
try:
import certifi
return ssl.create_default_context(cafile=certifi.where())
except ImportError: # pragma: no cover - depends on the environment
return ssl.create_default_context()
@staticmethod
def _parse_gwc(lines: Sequence[str]) -> dict:
"""Parse GWC libfile lines into roughness and height grids of A, k, freq.
After a title line, one line holds ``nrough nhgt nsec`` and the next two
the roughnesses and heights. Then follow ``(nrough, 2 * nhgt + 1, nsec)``
values, the frequencies and then A and k at each height.
"""
nrough, nhgt, nsec = map(int, lines[1].split())
roughnesses = np.array(lines[2].split(), dtype=float)
heights = np.array(lines[3].split(), dtype=float)
data = np.array([line.split() for line in lines[4:]], dtype=float).reshape(
(nrough, nhgt * 2 + 1, nsec)
)
freq = data[:, 0] / data[:, 0].sum(axis=1, keepdims=True)
return {
"roughnesses": roughnesses,
"heights": heights,
"weibull_a": data[:, 1::2], # (nrough, nhgt, nsec)
"weibull_k": data[:, 2::2], # (nrough, nhgt, nsec)
"frequency": freq, # (nrough, nsec)
}
@staticmethod
def _interp_axis(values, coords, target: float, axis: int) -> np.ndarray:
"""Interpolate ``values`` linearly along ``axis`` at ``target``.
``target`` is clamped to the ends of ``coords``, like ``numpy.interp``.
"""
values = np.asarray(values, dtype=float)
coords = np.asarray(coords, dtype=float)
target = float(np.clip(target, coords[0], coords[-1]))
idx = int(np.clip(np.searchsorted(coords, target), 1, len(coords) - 1))
lo, hi = coords[idx - 1], coords[idx]
weight = 0.0 if hi == lo else (target - lo) / (hi - lo)
low = np.take(values, idx - 1, axis=axis)
high = np.take(values, idx, axis=axis)
return low * (1.0 - weight) + high * weight
@classmethod
def _interp_gwc(cls, gwc: dict, roughness, height: float):
"""Interpolate A and k in roughness then height, and freq in roughness.
Roughness enters the log law as ``ln z0``, so the climates are
interpolated linearly in ``ln z0``. Libfiles list water as ``0``, which
WAsP takes as ``0.0002`` m. ``roughness`` is one length or one per
sector, and each sector takes the climate at its own.
"""
rough = np.log(np.maximum(np.asarray(gwc["roughnesses"], float), _WATER))
n_sectors = gwc["frequency"].shape[1]
targets = np.log(
np.maximum(
np.broadcast_to(np.asarray(roughness, float), (n_sectors,)), _WATER
)
)
heights = gwc["heights"]
a, k, freq = np.empty(n_sectors), np.empty(n_sectors), np.empty(n_sectors)
for target in np.unique(targets):
sectors = targets == target
for out, name in ((a, "weibull_a"), (k, "weibull_k")):
values = cls._interp_axis(gwc[name], rough, target, axis=0)
out[sectors] = cls._interp_axis(values, heights, height, axis=0)[
sectors
]
freq[sectors] = cls._interp_axis(gwc["frequency"], rough, target, axis=0)[
sectors
]
return a, k, freq / freq.sum()
def _sector_probabilities(self, n_sectors: int) -> np.ndarray:
return self._gwa_freq
def _sector_mean_speed(self) -> np.ndarray:
return self._gwa_mean
[docs]
class Dimitrov(_ArgumentRepr, _Resource):
"""The coverage design of Dimitrov et al. (2018), for training data.
Cases spread evenly over bounds that hold nearly all measured site
conditions, so a surrogate trained on them seldom extrapolates
(https://doi.org/10.5194/wes-3-767-2018). Check that the conditions it is
applied to lie inside them. The cases are not weighted by how often they
occur, so do not use the design for the energy, loads or statistics of a
real site.
Every case draws each variable within bounds that depend on the variables
before it, uniformly unless stated below. The wind speed at ``reference_height`` in meters
spreads over the sampled range, ``u_min`` to ``u_max`` of
:meth:`~WindResource.sample`, and the direction over ``n_sectors``
sectors. ``sigma_u`` lies between ``minimum_intensity U`` and
``reference_intensity (6.8 + 0.75 U + 3 (10 / U)^2)``, capped at
``maximum_intensity U`` when given. The Mann length scale lies between
``max(7.5, 15 |alpha|^-2/3)`` and 275 m and the anisotropy between 1 and 5.
Yaw and tilt lie within 10 degrees and air density between 1.1 and
1.35 kg/m^3. Veer across the rotor scales with its size, so it is drawn
only with ``rotor_diameter`` D in meters, between ``-0.1 D (5 / U)`` and
``min(60 sin|latitude|, D (5 / U)^2)`` degrees. These bounds are those of
Dimitrov et al.
Their power-law shear bounds are replaced by the log law, which OOWL
bounds instead. The roughness length is log-uniform over
``roughness_range`` in meters. ``1/L`` lies between
``-kappa g heat_flux / (T0 u*^3)`` and
``kappa g temperature_scale / (T0 u*^2)``, with ``u*`` over the smoothest
surface, so ``heat_flux`` in K m/s bounds unstable air and
``temperature_scale`` in K stable air. The shear of each case follows from
its roughness and stability.
"""
roughness_length: Roughness
stability: Stability
turbulence: Turbulence
mann: MannCoverage
def __init__(
self,
*,
minimum_intensity: float = 0.025,
reference_intensity: float = 0.18,
maximum_intensity: float | None = None,
roughness_range: tuple[float, float] = (1e-5, 1.0),
heat_flux: float = 0.3,
temperature_scale: float = 0.3,
rotor_diameter: float | None = None,
latitude: float = 50.0,
n_sectors: int = 12,
reference_height: float = 100.0,
):
smoothest, roughest = _validate_finite_pair("roughness_range", roughness_range)
if not 0.0 < smoothest < roughest:
raise ValueError("roughness_range must be positive and increasing")
self.rotor_diameter = (
None
if rotor_diameter is None
else _validate_positive("rotor_diameter", rotor_diameter)
)
self.roughness_range = (smoothest, roughest)
turbulence = Turbulence(
minimum_intensity, reference_intensity, maximum_intensity
)
stability = Stability(
heat_flux, temperature_scale, minimum_roughness_length=smoothest
)
inflow = Inflow(self.rotor_diameter, latitude)
self.latitude = inflow.latitude
self.minimum_intensity = turbulence.minimum_intensity
self.reference_intensity = turbulence.reference_intensity
self.maximum_intensity = turbulence.maximum_intensity
self.heat_flux = stability.heat_flux
self.temperature_scale = stability.temperature_scale
self.n_sectors = _validate_integer("n_sectors", n_sectors)
directions = UniformDirection()
self.theta = directions.sector_angles(self.n_sectors)
self.probability = directions.probabilities(self.n_sectors)
self.reference_height = reference_height
self.roughness_length = Roughness(smoothest, roughest)
self.stability = stability
self.turbulence = turbulence
self.mann = MannCoverage()
self._inflow = inflow
self._speed = _UniformSpeed(self.n_sectors)
self._setup_chain()
stability._checked_height(self.reference_height)
if (
self.rotor_diameter is not None
and 0.5 * self.rotor_diameter >= self.reference_height
):
raise ValueError("the rotor radius must lie below reference_height")
def _speed_record(self, wind_speed: np.ndarray) -> dict:
"""Return the design's speeds as a direction by speed probability table.
Each speed carries the probability of the cell around it, bounded by
the midpoints to its neighbors and by the first and last speed.
"""
edges = np.concatenate(
[wind_speed[:1], 0.5 * (wind_speed[1:] + wind_speed[:-1]), wind_speed[-1:]]
)
cells = np.diff(edges) / (edges[-1] - edges[0])
return {
"probability": {
"data": (self.probability[:, None] * cells[None, :]).tolist(),
"dims": ["wind_direction", "wind_speed"],
}
}
