"""Base classes and containers for synthetic site conditions.
Direction sectors are evenly spaced at ``theta_i = 360 * i / n`` degrees.
"""
from __future__ import annotations
import json
from collections.abc import Sequence
from dataclasses import KW_ONLY, dataclass, field, fields
from pathlib import Path
from typing import TYPE_CHECKING, Any
import numpy as np
from .._validation import _validate_integer
from .stability import _checked_height, _profile
if TYPE_CHECKING: # pragma: no cover
import matplotlib.axes
import matplotlib.figure
__all__ = ["SiteConditions", "sector_angles"]
def sector_angles(n_sectors: int) -> np.ndarray:
"""Return the ``n_sectors`` evenly spaced sector centers ``360 * i / n``."""
_validate_integer("n_sectors", n_sectors)
return 360.0 * np.arange(int(n_sectors), dtype=float) / int(n_sectors)
# Optional per-case Mann parameters, filled by a mann link.
_MANN_FIELDS = ("mann_length_scale", "mann_anisotropy")
# Optional per-case inflow variables drawn by the Dimitrov design.
_INFLOW_FIELDS = (
"wind_veer",
"yaw_misalignment",
"inflow_tilt",
"air_density",
)
#: Version of the ``.npz`` layout :meth:`SiteConditions.save` writes.
_SAVED_VERSION = 1
[docs]
@dataclass(eq=False)
class SiteConditions:
"""Sampled inflow cases and their sampling metadata.
Each case holds a wind-direction sector center, wind speed, ``sigma_u``,
``TI = sigma_u / U``, inverse Obukhov length, the site's
``reference_height`` repeated per case, and the ``roughness_length`` each
case sees, which varies with wind speed for a Charnock roughness and is drawn
per case by :class:`~oowl.resources.Dimitrov`. The wind
profile is the log law of that roughness and stability, so
:meth:`wind_speed_at` gives the speed at any height and
``shear_exponent`` is its local power-law exponent at the reference
height. ``weight`` is each case's probability mass given a speed in the
sampled range and sums to one. ``metadata["in_range_fraction"]`` is that
range's share of time.
A ``mann`` link fills ``mann_length_scale`` in meters and
``mann_anisotropy``. A :class:`~oowl.resources.Dimitrov` design fills
``yaw_misalignment`` and ``inflow_tilt`` in degrees and ``air_density`` in
``kg/m^3``, and with a rotor diameter ``wind_veer`` across the rotor in
degrees. Otherwise they are ``None``.
"""
wind_direction: np.ndarray
wind_speed: np.ndarray
sigma_u: np.ndarray
turbulence_intensity: np.ndarray
inverse_obukhov_length: np.ndarray
shear_exponent: np.ndarray
reference_height: np.ndarray
roughness_length: np.ndarray
weight: np.ndarray
_: KW_ONLY
metadata: dict = field(default_factory=dict)
mann_length_scale: np.ndarray | None = None
mann_anisotropy: np.ndarray | None = None
wind_veer: np.ndarray | None = None
yaw_misalignment: np.ndarray | None = None
inflow_tilt: np.ndarray | None = None
air_density: np.ndarray | None = None
def __repr__(self) -> str:
return f"SiteConditions(n_cases={len(self.wind_speed)})"
def __eq__(self, other: object) -> bool:
"""Whether ``other`` holds the same cases, arrays and metadata."""
if other.__class__ is not self.__class__:
return NotImplemented
def same(left, right) -> bool:
# An optional field may be an array on one side and None on the other.
if isinstance(left, np.ndarray) or isinstance(right, np.ndarray):
return (
isinstance(left, np.ndarray)
and isinstance(right, np.ndarray)
and np.array_equal(left, right)
)
return bool(left == right)
return all(
same(getattr(self, item.name), getattr(other, item.name))
for item in fields(self)
)
def __post_init__(self):
arrays = {
name: np.asarray(getattr(self, name), dtype=float)
for name in (
"wind_direction",
"wind_speed",
"sigma_u",
"turbulence_intensity",
"inverse_obukhov_length",
"shear_exponent",
"reference_height",
"roughness_length",
"weight",
*_MANN_FIELDS,
*_INFLOW_FIELDS,
)
if getattr(self, name) is not None
}
# A Mann link fills both, while the inflow fields come without veer
# when no rotor diameter is given.
if len({name in arrays for name in _MANN_FIELDS}) > 1:
raise ValueError(f"{' and '.join(_MANN_FIELDS)} must be given together")
n = arrays["wind_direction"].shape
if any(arr.shape != n or arr.ndim != 1 for arr in arrays.values()):
raise ValueError("all site-condition arrays must be 1-D and equal length")
for name, arr in arrays.items():
if not np.all(np.isfinite(arr)):
raise ValueError(f"{name} must contain only finite values")
setattr(self, name, arr)
if np.any(self.weight < 0):
raise ValueError("weight must be non-negative")
if not np.isclose(float(self.weight.sum()), 1.0, rtol=1e-10, atol=1e-12):
raise ValueError("weight must sum to one")
if np.any(self.wind_speed <= 0):
raise ValueError("wind_speed must be positive")
if np.any(self.sigma_u < 0):
raise ValueError("sigma_u must be non-negative")
if np.any(self.roughness_length <= 0):
raise ValueError("roughness_length must be positive")
if np.any(self.reference_height != self.reference_height[0]):
raise ValueError("reference_height must be constant across one site")
if np.any(self.reference_height <= self.roughness_length):
raise ValueError("reference_height must lie above roughness_length")
expected_ti = self.sigma_u / self.wind_speed
if not np.allclose(
self.turbulence_intensity,
expected_ti,
rtol=1e-10,
atol=1e-12,
):
raise ValueError("turbulence_intensity must equal sigma_u / wind_speed")
if not isinstance(self.metadata, dict):
raise TypeError("metadata must be a dictionary")
try:
json.dumps(self.metadata, allow_nan=False)
except (TypeError, ValueError) as error:
raise ValueError("metadata must be finite and JSON-serializable") from error
@property
def n_cases(self) -> int:
"""Number of sampled inflow cases."""
return int(len(self.wind_direction))
[docs]
def stability_parameter(self, height: float) -> np.ndarray:
"""Return ``z/L`` at ``height`` meters for every case, positive when stable."""
return _checked_height(height) * self.inverse_obukhov_length
[docs]
def wind_speed_at(self, height: float) -> np.ndarray:
"""Return every case's wind speed in ``m/s`` at ``height`` meters.
The Monin-Obukhov log law of the case, scaled to its speed at the
reference height,
.. math::
U(z) = U(z_r) \\frac{\\ln(z/z_0) - \\Psi(z/L) + \\Psi(z_0/L)}
{\\ln(z_r/z_0) - \\Psi(z_r/L) + \\Psi(z_0/L)},
with the stable gradient held beyond
:data:`~oowl.resources.profile.STABLE_LIMIT`. ``height`` must lie above every
case's roughness length.
"""
height = _checked_height(height)
reference = float(self.reference_height[0])
at_height = _profile(height, self.roughness_length, self.inverse_obukhov_length)
at_reference = _profile(
reference, self.roughness_length, self.inverse_obukhov_length
)
return self.wind_speed * at_height / at_reference
[docs]
def plot(
self,
ax: "Sequence[matplotlib.axes.Axes] | np.ndarray | None" = None,
**kwargs: Any,
) -> np.ndarray:
"""Plot the wind rose, turbulence, stability and wind-speed panels.
Draws four panels, so ``ax`` is a sequence of four Axes, the first
polar, and a new figure is created without it. Returns the array of the
four Axes.
"""
from ..plotting import _plot_site_conditions
return _plot_site_conditions(self, ax=ax, **kwargs)
[docs]
def save(self, path: str | Path) -> Path:
"""Write the cases to ``path``, a NumPy ``.npz`` file, and return the path.
Each array is stored under its field name, the optional ones only
when set, and the metadata as JSON text under ``metadata``, so
:meth:`load` rebuilds equal conditions. Raises ``ValueError`` for
another suffix.
"""
path = Path(path)
if path.suffix.lower() != ".npz":
raise ValueError(
f"site conditions are saved as a .npz file, not {path.name}"
)
arrays = {
item.name: getattr(self, item.name)
for item in fields(self)
if item.name != "metadata" and getattr(self, item.name) is not None
}
path.parent.mkdir(parents=True, exist_ok=True)
np.savez(
path,
schema_version=np.asarray(_SAVED_VERSION),
metadata=np.asarray(json.dumps(self.metadata)),
**arrays,
)
return path
[docs]
@classmethod
def load(cls, path: str | Path) -> SiteConditions:
"""Read the conditions :meth:`save` wrote to ``path``.
Raises ``ValueError`` for a file of another version or one missing a
case array, and the checks of the constructor apply to what it holds.
"""
path = Path(path)
with np.load(path, allow_pickle=False) as data:
stored = dict(data)
if int(stored.pop("schema_version", -1)) != _SAVED_VERSION:
raise ValueError(f"{path.name} is not site conditions OOWL saved")
known = {item.name for item in fields(cls)}
unknown = sorted(set(stored) - known)
if unknown:
raise ValueError(f"{path.name} holds unknown arrays {unknown}")
required = [
item.name
for item in fields(cls)
if item.name not in stored
and item.name != "metadata"
and item.kw_only is False
]
if required:
raise ValueError(f"{path.name} is missing {required}")
metadata = json.loads(str(stored.pop("metadata", "{}")))
return cls(**stored, metadata=metadata)
