# (c) 2026 DTU Wind Energy
"""Synthetic mesoscale surface fields.
Stability-aware wind climate code that starts from raw mesoscale output needs
the surface fields the surface-layer temperature scale is derived from -
pressure, temperature, friction velocity, the heat fluxes, the boundary layer
height and the surface roughness. This module creates a synthetic version of
those, so that such code can be exercised without a mesoscale model run.
The derived stability variables themselves, ``temp_scale`` and
``ustar_over_pblh``, belong to a time series wind climate and are generated by
:func:`windkit.create_tswc` through its ``fields`` argument.
:func:`create_surface_fields` takes the same argument, for the case where the
two groups have to be correlated with each other. Both share the engine in
:mod:`windkit.wind_climate._correlated_fields`, so the values follow bounded
multi-harmonic sinusoids whose amplitude and phase differ per point, and
within a point they carry the cross-correlations of the reference run.
"""
__all__ = [
"create_surface_fields",
]
import numpy as np
import pandas as pd
import xarray as xr
from windkit.wind_climate._correlated_fields import (
_column_index,
_column_of,
_correlated_fields,
_seconds_since_start,
_stability_field_ranges,
)
from windkit.wind_climate.time_series_wind_climate import _create_tswc_core
from windkit.xarray_structures.empty import (
_native_batch,
_native_batch_horizontal,
_native_flat_coords,
_native_store,
)
from windkit.xarray_structures.metadata import (
_ALL_VARS_META,
_TS_ATTRS,
_update_history,
_update_var_attrs,
)
# One month at hourly resolution. Small enough to stay a cheap fixture, and
# long enough to clear the two-day span the correlated-field engine needs.
_DEFAULT_DATE_RANGE = pd.date_range(
"2001-01-01", "2001-02-01", freq="h", inclusive="left"
)
# Ranges for the mesoscale surface fields, used when the caller does not
# override them. Each entry is ``(low, high)`` in the variable's own units,
# taken as the 0.1st and 99.9th percentile of the same Cabauw reference and
# rounded. ZNT is absent from the reference, so its entry is a plausible
# range of land cover roughness instead.
_SURFACE_FIELD_RANGES = {
"PSFC": (97500.0, 104100.0),
"T2": (261.0, 304.0),
"UST": (0.01, 1.15),
"HFX": (-100.0, 265.0),
"LH": (-25.0, 440.0),
"PBLH": (19.0, 2100.0),
"Q2": (0.0011, 0.0152),
"ZNT": (5e-3, 0.1),
}
[docs]
def create_surface_fields(
output_locs,
date_range=None,
not_empty=True,
seed=9876538,
ranges=None,
elevation_range=(0.0, 200.0),
fields=None,
**tswc_kwargs,
):
"""Create synthetic mesoscale surface fields as a time series.
Produces the raw model output variables from which the surface-layer
temperature scale is derived, using the WRF naming convention, alongside
the wind time series from :func:`create_tswc`.
Time-varying variables: ``PSFC``, ``T2``, ``UST``, ``HFX``, ``LH``,
``PBLH``, ``Q2``, ``ZNT``, ``wind_speed``, ``wind_direction``.
Time-constant variables: ``LANDMASK``, ``elevation``.
The stability variables these fields feed, ``temp_scale`` and
``ustar_over_pblh``, are not part of the mesoscale set. They can be added
with the ``fields`` argument, which is resolved exactly as
:func:`create_tswc` resolves its own; asking for them here generates them
together with the surface fields, in one correlated draw, so that a point
holds the reference correlations between the two groups as well - the
``temp_scale``/``HFX`` pair among them.
Every point is over land: the values are drawn from an onshore reference,
so ``LANDMASK`` is 1 everywhere and the fields carry no offshore
behaviour. Code that has to exercise a land/sea
contrast should set ``LANDMASK`` and the surface fields of the sea points
itself.
These are all surface variables: every one of them declares
``height_dim_policy: "forbidden"``, so none carries a height dimension in
the returned dataset, whatever the structure of ``output_locs``. Each is
drawn once per ``(west_east, south_north)`` column, so no two columns hold
the same values — a fixture that repeated them could not detect code which
loses track of which location a value belongs to. A column is identified
by exact equality of its two coordinates, which holds for every structure
windkit builds. ``wind_speed`` and ``wind_direction`` come from
:func:`create_tswc` and keep their height dimension.
Within a point the fields carry the cross-correlations of the reference
site, including with ``wind_speed``.
Parameters
----------
output_locs : xarray.Dataset
Output geospatial information.
date_range : pandas.DatetimeIndex or None
Time range. If None, one month of 2001 at hourly resolution is used,
which is small enough to stay a cheap fixture. It must span more than
two days, whatever its resolution.
not_empty : bool
If True, the surface fields are filled with synthetic data. If False,
every variable is NaN, as :func:`create_tswc` does. Defaults to True.
seed : int
Seed for the random data, defaults to 9876538.
ranges : dict or None
Overrides for the ``(low, high)`` bounds of individual surface fields.
Keys not given fall back to the module defaults, and keys outside the
surface fields are rejected - the stability variables are bounded
by ``fields`` instead. Each entry must have ``low < high``;
neither the ordering nor the width of a pair is checked, so a reversed
pair flips the field about its mid-range and an equal pair makes it
constant.
elevation_range : tuple of float
Bounds of the terrain elevation [m], with ``low < high``. Defaults to
``(0.0, 200.0)``.
fields : str, sequence of str or None
Stability variables to generate alongside the surface fields, as in
:func:`create_tswc`: ``None`` for none of them, ``"stability"`` or
``"all"`` for every one, or the names to generate. Unlike passing them
to :func:`create_tswc` separately, these are drawn together with the
surface fields, so the two groups are correlated with each other.
**tswc_kwargs
Further keyword arguments passed to :func:`create_tswc`.
Returns
-------
ds : xarray.Dataset
Time series of mesoscale surface fields.
Raises
------
ValueError
If ``ranges`` holds a key that is not a surface field, if ``fields``
names something that is not a stability variable, or if ``date_range``
holds too few time steps or spans two days or less.
Notes
-----
Every variable added here is a surface variable, so all of them are stored
over the horizontal batch, never a separate height axis: float32,
C-contiguous, with ``time`` trailing that batch, the reverse of the
compiled core's declared argument order (windkit#807). The wind comes from
:func:`create_tswc` and is stored as that function documents. Coordinates
stay float64. Neither the dtype nor the dimension order is a supported API
contract: ``create_surface_fields`` is a test-data utility.
Examples
--------
>>> import pandas as pd
>>> import windkit as wk
>>> locs = wk.spatial.create_dataset(
... [10.0, 10.1], [55.0, 55.1], [100.0, 100.0], 4326, struct="point"
... ).drop_vars("output")
>>> ds = wk.create_surface_fields(locs)
>>> "PBLH" in ds and "LANDMASK" in ds
True
>>> ds.sizes["time"]
744
"""
unknown = sorted(set(ranges or {}) - set(_SURFACE_FIELD_RANGES))
if unknown:
raise ValueError(
f"Unknown surface field(s) {unknown}, cannot add to result. The keys "
f"of 'ranges' must be a subset of {sorted(_SURFACE_FIELD_RANGES)}."
)
# The stability variables come last, so that a call without ``fields``
# draws exactly what it drew before they could be asked for here.
field_ranges = {
**_SURFACE_FIELD_RANGES,
**(ranges or {}),
**_stability_field_ranges(fields),
}
if date_range is None:
date_range = _DEFAULT_DATE_RANGE
# The core rather than the public constructor, for its float64 wind: the
# dataset stores it float32 (windkit#807), and driving the fields off the
# rounded values would make them depend on the storage dtype.
ds, internals = _create_tswc_core(
output_locs,
date_range=date_range,
not_empty=not_empty,
seed=seed,
_return_internals=True,
**tswc_kwargs,
)
seconds = _seconds_since_start(pd.DatetimeIndex(ds["time"].values))
# Every variable added here is a surface variable, so all of them are
# built over the horizontal batch and never carry a height dimension.
# `create_tswc`'s own wind keeps the full batch, and its height.
dims, shape = _native_batch(output_locs)
horiz_dims, horiz_shape = _native_batch_horizontal(output_locs)
n_horiz = int(np.prod(horiz_shape, dtype=int))
if not not_empty:
# Mirrors create_tswc and create_wasp_site_effects: an empty dataset
# carries the full variable set, all of it NaN.
surface = {
"LANDMASK": np.full(n_horiz, np.nan),
"elevation": np.full(n_horiz, np.nan),
}
time_varying = {
name: np.full((len(seconds), n_horiz), np.nan) for name in field_ranges
}
return _finalize(ds, surface, time_varying, horiz_dims, horiz_shape)
rng = np.random.default_rng(seed + 2)
# Drawn once per horizontal column, then spread over the horizontal batch:
# a permutation where height is its own dimension, and a repeat over the
# heights above each column in a `point` structure, where it is not.
we, sn = _native_flat_coords(output_locs, dims, shape)
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)
)
n_col = len(base_of_column)
surface = {
"LANDMASK": np.ones(n_horiz),
"elevation": rng.uniform(*elevation_range, n_col)[column_of_location],
}
# The wind driving them is taken at the lowest height of each column, the
# closest thing to the surface wind the fields respond to.
generated = _correlated_fields(
field_ranges,
seconds,
n_col,
rng,
internals["wind_speed"][:, base_of_column],
)
time_varying = {
name: values[:, column_of_location] for name, values in generated.items()
}
return _finalize(ds, surface, time_varying, horiz_dims, horiz_shape)
def _finalize(ds, surface, time_varying, horiz_dims, horiz_shape):
"""Store the surface fields over the horizontal batch and attach metadata.
Nothing is unstacked: the height dimension every surface variable
declares ``"forbidden"`` is never introduced in the first place, because
these are built over a batch that has no height axis to begin with.
"""
for name, values in surface.items():
out_dims, array = _native_store(values, horiz_dims, horiz_shape, core_dims=())
ds[name] = xr.DataArray(array, dims=out_dims)
for name, values in time_varying.items():
out_dims, array = _native_store(
values, horiz_dims, horiz_shape, core_dims=("time",)
)
ds[name] = xr.DataArray(array, dims=out_dims)
meta = {
**_TS_ATTRS,
**{k: _ALL_VARS_META[k] for k in ds.data_vars if k in _ALL_VARS_META},
}
return _update_history(_update_var_attrs(ds, meta))
