"""Binned Wind Climate Module
Extends windkit Binned Wind Climate module to add a function for creating a weibull
distribution out of the binned wind climate using the WAsP method. It contains a much
faster fortran based binning routine that can be used for example in combination with
mesoscale model output, where speed is crucial.
It also contains functions for calculating stability information efficiently using
a binning algorithm which computes the conditional mean stability based on the wind
histogram in a single pass without storing the whole time series of wind speed.
"""
__all__ = [
"bwc_from_tswc",
"calc_temp_scale",
"stability_histogram",
"create_histogram_z0",
"mean_z0",
"baroclinicity_histogram",
"bwc_resample_like",
"bwc_resample_sectors",
"bwc_resample_wsbins_like",
]
import logging
import numpy as np
import xarray as xr
from scipy import interpolate
from windkit.wind_climate.binned_wind_climate import (
_create_time_attributes,
_validate_bwc_wrapper_factory,
_bwc_from_counts,
)
from windkit.wind_climate.time_series_wind_climate import (
_validate_tswc_wrapper_factory,
WS,
WD,
DIM_TIME,
)
from windkit import create_sector_coords
from windkit import create_wsbin_coords
from windkit.spatial import _spatial_stack, _spatial_unstack, to_point, set_crs
from windkit.xarray_structures.metadata import (
_BWC_ATTRS,
_HIS_ATTRS,
_HIS_STAB_ATTRS,
_TS_STAB_ATTRS,
_MEAN_STAB_ATTRS,
_HIS_BARO_ATTRS,
_MEAN_BARO_ATTRS,
_update_var_attrs,
_ALL_VARS_META,
_update_history,
)
from pywasp.core import Rvea0288, Rvea0326 # pylint:disable=no-name-in-module
logger = logging.getLogger(__name__)
# fortran calls
_CREATE_HIST = Rvea0288.histograms.add_histogram_points
_CALC_TSCALE = Rvea0326.calc_points.tstar_points
_NORM_HIST = Rvea0288.histograms.normalize_histogram_points
_COND_HIST = Rvea0288.histograms.conditional_hist
_MEAN_HIST = Rvea0288.histograms.mean_from_hist
_STD_HIST = Rvea0288.histograms.std_from_hist
_RESAMPLE_HIST = Rvea0288.histograms.histogram_resample_sectors
_RESAMPLE_HIST_WSBINS = Rvea0288.histograms.histogram_resample_wsbins
_RESAMPLE_HIST_WSBINS_AK = Rvea0288.histograms.histogram_resample_wsbins_ak
_DGDZ = Rvea0288.baroclin.compute_dgdz
_DGDZ_TIME = Rvea0288.baroclin.compute_dgdz_t
# constants
RTD = 180 / np.pi
[docs]
@_validate_tswc_wrapper_factory(run_extra_checks=False)
def bwc_from_tswc(
ds,
hist=None,
normalize=True,
revert_to_original=True,
wsbin_width=1,
n_wsbins=30,
n_sectors=12,
):
"""Add timeseries to histogram
Converts a time series with wind speed & direction to derive
summed histogram of the wind vector.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
ds: xr.Dataset
PyWAsP Timeseries dataset. Can have any of the pywasp
spatial data structures
hist: xr.Dataset, optional
Histogram with dimensions point, wsbin, sector
containg a count
Default is None, creating a histogram using wsbin_width, n_wsbins, and n_sectors
normalize: boolean
Normalize histogram for each sector, storing the frequency per sector in "wdfreq"
variable, Default True
revert_to_original: boolean
Return the histogram input in the same spatial struction as ds? If false, keeps
the data is kept in stacked point, which is more efficient when doing large
calculations for numerical wind atlases.
wsbin_width: float
width of wind speed bins
n_wsbins: int
Number of wind speed bins
n_sectors: int
Number of sectors (wind direction bins)
Returns
-------
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
with the values from ds-timeseries added.
"""
logger.debug(
f"Convert time series ({ds[WS].sizes[DIM_TIME]} samples) into histogram"
)
if hist is None:
hist = _init_histogram_wind(ds, wsbin_width, n_wsbins, n_sectors)
dss = _spatial_stack(ds, copy=False).transpose(..., DIM_TIME, "point")
add = np.full(dss[WS].shape, 1.0, order="F", dtype=np.float32)
if "point" not in hist.dims:
hist = _spatial_stack(hist, copy=False)
hist = hist.transpose(..., "wsbin", "sector", "point")
wv_count, count_bins_exceeded = _create_hist(
dss[WS].fillna(-9999.0),
dss[WD].fillna(-9999.0),
0,
hist.wsbin,
hist.sector,
add,
hist["wv_count"].values,
hist["count_bins_exceeded"].values,
)
nr_fails = count_bins_exceeded[count_bins_exceeded > 0].shape[0]
if nr_fails > 0:
logger.warning(
f"""The wind speeds you provided
exceeded the range used for binning ({np.max(hist.wsbin.values)})
in {nr_fails} points"""
)
var_names = ("wv_count", "count_bins_exceeded")
for name in var_names:
hist[name] = (hist[name].dims, eval(name), _HIS_ATTRS[name])
_update_var_attrs(hist[name], _HIS_ATTRS)
hist[name].attrs["_pwio_data_is_2d"] = False
dic_attrs = _create_time_attributes(ds, hist=hist)
hist = hist.assign_attrs(dic_attrs).drop_vars("point", errors="ignore")
if normalize:
hist = _normalize_hist(hist)
if revert_to_original:
hist = _spatial_unstack(hist)
hist = _update_var_attrs(hist, _BWC_ATTRS)
return _update_history(hist)
def _create_hist(wspd, wdir, isuv, wsbin, wdsec, add, wv_count, count_bins_exceeded):
"""creating histogram from wind speed and wind dir
Wrapper around fortran call that loops over all points
and generates 2d histograms from wind speed and direction.
If the wind speed or wind direction is outside of the range
used for binning it adds 1 to `count_bins_exceeded` for each point
Parameters
----------
wspd: array of float
Array with wind speeds dimensions (`time`, `point`)
wdir: array of float
Array with meteorological wind direction dimensions (`time`, `point`)
isuv: int {0,1}
Are the input values wind speed / wind direction (0)
or u / v (1)?
wsbin: array of float
center of wind speed bins used for creating histogram
wdsec: array of float
center of wind direction bins used for creating histogram
add: array of float
Array with values that will be added to histogram at each point
wv_count: array of float
Array with histogram with dimensions #wsbin, #wdsec, #point
count_bins_exceeded: array of int
Array with dimension point
Returns
-------
wv_count: array of float
Array with histogram with dimensions #wsbin, #wdsec, #point
count_bins_exceeded: array of int
Array with dimension point
"""
_CREATE_HIST(wspd, wdir, isuv, wsbin, wdsec, add, wv_count, count_bins_exceeded)
return (wv_count, count_bins_exceeded)
[docs]
def calc_temp_scale(ds, revert_to_original=True, enhmix=1.0, use_humidity=False):
"""Calculate temperature scale
This calculates a surface layer temperature scale that can
be used to model the wind profile. Using the enhmix can be used
to control mixing over land, because ERA5 has a large bias
in $u_*$ over land. If you supply enhmix you will have to
add a landmask because it is only applied over land. If use_humidity
is true the humidity is used for calculating the air density.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
ds : xarray.Dataset
Dataset containing variables
["PSFC", "T2", "UST", "HFX", "LH", "PBLH"]
Can have any pywasp spatial structure
revert_to_original: boolean
Return the histogram input in original format from the
input :class:`xarray.Dataset`? If false, keeps the data in stacked point
format which is more efficient when doing large
calculations for numerical wind atlases.
enhmix : float
Decrease mixing over land with this factor. If different from
1.0, requires a land mask variable "LANDMASK" since it is only
applied over land.
use_humidity: bool
If True, use the humdity for calculating the air density. If
true, requires variable "elevation" (for the terrain elevation)
and "Q2" to be present on the dataset. By default False
Returns
-------
hist : xr.DataArray
DataArray 'temp_scale' with same dimensions as input
"""
dss = _spatial_stack(ds, copy=revert_to_original).transpose(..., DIM_TIME, "point")
req_vars = ["PSFC", "T2", "UST", "HFX", "LH", "PBLH"]
if enhmix != 1.0:
req_vars = req_vars + ["LANDMASK"]
if use_humidity:
req_vars = req_vars + ["elevation", "Q2"]
if "LANDMASK" not in req_vars:
logger.debug("Adding dummy variable LANDMASK")
# land mask not required in dataset, so we just create a dummy array
dss["LANDMASK"] = xr.full_like(dss["point"], 0)
if "Q2" not in req_vars:
# humidity not required in dataset, so we just create a dummy array
logger.debug("Adding dummy variables Q2 and elevation")
dss["elevation"] = xr.full_like(dss["point"], 0)
dss["Q2"] = xr.full_like(dss["UST"], 0)
try:
logger.debug("Calculating temperature scale")
tstar = _CALC_TSCALE(
dss["PSFC"],
dss["T2"],
dss["UST"],
dss["HFX"],
dss["LH"],
dss["PBLH"],
dss["elevation"],
dss["Q2"],
dss["LANDMASK"],
enhmix,
int(use_humidity),
)
except KeyError as e:
raise KeyError(f"Required variable {e} not found in dataset!")
dss["temp_scale"] = (dss.UST.dims, tstar)
dss = dss[["temp_scale"]]
_update_var_attrs(dss, _TS_STAB_ATTRS)
if revert_to_original:
res = _spatial_unstack(dss)["temp_scale"]
else:
res = dss["temp_scale"]
return _update_history(res)
[docs]
def stability_histogram(
ds,
hist=None,
finalize=True,
revert_to_original=True,
wsbin_width=1,
n_wsbins=40,
n_sectors=12,
percentile=0.5,
wv_count=None,
landmask=None,
):
"""Add timeseries to existing histogram
Uses the stability parameters "temp_scale" and "ustar_over_pblh".
See the paper "Using observed and modelled heat fluxes for
improved extrapolation of wind distributions" for the definitions
of these variables. The temp_scale variable is defined as $T_*$
in that paper, whereas the "ustar_over_pblh" should be defined as
$u_*/pblh$ as input for this function. This is done to give priority
to the smallest pblh that are most important for the wind profile.
The output is reversed ($**-1$) to give the input format as required
by WAsP ($pblh/u_*$).
The time dimension is copied so that one can see what was the
period that was used to generate the histogram.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
ds : xarray.Dataset
Dataset containing variables ['wind_speed', 'wind_direction',
'temp_scale', 'ustar_over_pblh']
Can have any pywasp spatial structure
hist : xarray.Dataset
Histogram with dimensions point, wsbin, sector
finalize : boolean
Convert the stabilility histogram into a mean and
standard deviation?
revert_to_original : boolean
Return the histogram input in original format from the
input ``ds``? If false, keeps the data in stacked point
format which is more efficient when doing large
calculations for numerical wind atlases.
wsbin_width : float
width of wind speed bins
n_wsbins : int
Number of wind speed bins
n_sectors : int
Number of sectors (wind direction bins)
kwargs :
Other keyword arguments passed on to _finalize_meso
Returns
-------
hist : xarray.Dataset
Histogram with dimensions point, wsbin, sector
with the values from ds-timeseries added.
"""
if hist is None:
hist = _init_histogram_stab(ds, wsbin_width, n_wsbins, n_sectors)
ds_sub = ds[[WS, WD, "temp_scale", "ustar_over_pblh"]]
dss = _spatial_stack(ds_sub, copy=revert_to_original).transpose(
..., DIM_TIME, "point"
)
hist = _spatial_stack(hist, copy=revert_to_original).transpose(
..., "wsbin", "sector", "point"
)
# error arrays are inout in fortran, so we must create copies
# to not add the error counter 3 times
ierr_sum_squared_temp_scale = hist["count_bins_exceeded"].copy().values
ierr_sum_pblh = hist["count_bins_exceeded"].copy().values
sum_temp_scale, count_bins_exceeded = _create_hist(
dss[WS],
dss[WD],
0,
hist.wsbin,
hist.sector,
dss["temp_scale"],
hist["sum_temp_scale"].values,
hist["count_bins_exceeded"].values,
)
sum_squared_temp_scale, count_bins_exceeded = _create_hist(
dss[WS],
dss[WD],
0,
hist.wsbin,
hist.sector,
dss["temp_scale"] ** 2,
hist["sum_squared_temp_scale"].values,
ierr_sum_squared_temp_scale,
)
sum_pblh, count_bins_exceeded = _create_hist(
dss[WS],
dss[WD],
0,
hist.wsbin,
hist.sector,
dss["ustar_over_pblh"],
hist["sum_pblh"].values,
ierr_sum_pblh,
)
var_names = (
"sum_temp_scale",
"sum_squared_temp_scale",
"sum_pblh",
"count_bins_exceeded",
)
for name in var_names:
hist[name] = (hist[name].dims, eval(name))
_update_var_attrs(hist[name], _HIS_STAB_ATTRS)
hist[name].attrs["_pwio_data_is_2d"] = True
dic_attrs = _create_time_attributes(ds_sub, hist=hist)
hist = hist.assign_attrs(dic_attrs)
if finalize:
if wv_count is None:
wv_count = bwc_from_tswc(
dss,
normalize=False,
revert_to_original=False,
wsbin_width=wsbin_width,
n_wsbins=n_wsbins,
n_sectors=n_sectors,
)
else:
wv_count = bwc_from_tswc(
dss, hist=wv_count, normalize=False, revert_to_original=False
)
res = _finalize_meso(
hist, wv_count["wv_count"], landmask, percentile=percentile
)
return _update_history(res)
if revert_to_original:
res = _spatial_unstack(hist)
return _update_history(res)
else:
res = hist
return _update_history(res)
def _init_histogram_stab(ds, wsbin_width=1, nws=40, nwd=12):
"""Initialize histogram
Initialize histogram with the variables needed for
a stability histogram. The dimensions of the spatial
structure are the same as the input ds and extended with
dimensions nws, nwd.
Parameters
----------
ds: xr.Dataset
Can have any of the pywasp spatial data structures
var_dic: dict
Dictonairy contain the variables that will be added to
the xarray container, including meta-data according to
the CF conventions
wsbin: float
width of wind speed bins
wdsec: float
width of wind direction bins
nws: int
Number of wind speed bins
nwd: int
Number of wind direction bins
Returns
-------
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
"""
hist = _spatial_stack(ds, copy=False)[["point"]]
var_names = ("sum_temp_scale", "sum_squared_temp_scale", "sum_pblh")
for var in var_names:
hist[var] = _init_histogram(ds, var, nwd=nwd, nws=nws, wsbin_width=wsbin_width)
hist["count_bins_exceeded"] = xr.zeros_like(hist["point"], dtype=np.int32)
return hist
def _init_histogram(ds, name, nwd=12, nws=None, wsbin_width=None):
"""Initialize generic histogram
Initialize histogram with the variables needed for
a wind speed and wind direction dependent histogram. The dimensions
of the spatial structure are the same as the input ds and extended
with dimensions nws, nwd.
Parameters
----------
ds: `xarray.Dataset`
Can have any of the pywasp spatial data structures
name: str
Name of the `xarray.DataArray` that is created
nwd: int
Number of wind direction bins
nws: int
Number of wind speed bins
wsbin_width: float
width of wind speed bins
Returns
-------
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
"""
all_dims = (_spatial_stack(ds, copy=False)["point"],)
all_dims = all_dims + (create_sector_coords(nwd),)
if wsbin_width is not None and nws is not None:
all_dims = all_dims + (create_wsbin_coords(nws, width=wsbin_width),)
template = xr.broadcast(*all_dims)
# combined point, ws_bin, wd_bin array filled with zeroes
hist = xr.zeros_like(template[0], dtype="float32")
# make sure we copy auxilary coordinates from sector and wind speed bins
for i in template[1:]:
hist = hist.assign_coords(i.coords)
hist.name = name
try:
_update_var_attrs(hist, {**_ALL_VARS_META})
except KeyError:
logging.warning("No metadata could be added to the variable!")
return hist
def _init_histogram_wind(ds, wsbin_width, nws, nwd):
"""Initialize histogram
Initialize histogram with the variables needed for
a wind vector count histogram. The dimensions of the spatial
structure are the same as the input ds and extended with
dimensions nws, nwd.
Parameters
----------
ds: xr.Dataset
Can have any of the pywasp spatial data structures
wsbin_width: float
width of wind speed bins
nws: int
Number of wind speed bins
nwd: int
Number of wind direction bins
Returns
-------
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
"""
hist = _spatial_stack(ds, copy=False)[["point"]]
hist["wv_count"] = _init_histogram(
ds, "wv_count", nwd=nwd, nws=nws, wsbin_width=wsbin_width
)
hist["count_bins_exceeded"] = xr.zeros_like(hist["point"], dtype=np.int32)
return hist
def _stab_mean_std(hist, wv_count, percentile=0.5):
"""Compute mean or standard deviation from histogram
Compute mean or standard deviation from histogram. The
variables are binned by sector and wind speed bin and
a mean or std. dev. is calulcated based on the highest
percentile of wind speeds for each sector.
Parameters
----------
hist: xr.Dataset
Dataset with variables `sum_temp_scale`, `sum_squared_temp_scale`
and `sum_pblh` and dimensions `wsbin`, `sector`
Returns
-------
clim: xr.Dataset
Dataset with mean mesoclimate variables per sector
"""
hist["wv_count"] = wv_count
hist = _spatial_stack(hist, copy=False)
cond_hist, freq = _COND_HIST(hist.sum_pblh, hist.wv_count, percentile)
pblhmean = _MEAN_HIST(cond_hist, freq)
cond_hist, freq = _COND_HIST(hist.sum_temp_scale, hist.wv_count, percentile)
cond_hist2, freq = _COND_HIST(
hist.sum_squared_temp_scale, hist.wv_count, percentile
)
tstarstd = _STD_HIST(cond_hist, cond_hist2, freq)
tstarmean = _MEAN_HIST(cond_hist, freq)
clim = hist[["point", "sector"]]
clim["mean_pblh_scale"] = (("sector", "point"), 1 / pblhmean)
clim["rms_temp_scale"] = (("sector", "point"), tstarstd)
clim["mean_temp_scale"] = (("sector", "point"), tstarmean)
return _spatial_unstack(clim)
def create_histogram_z0(
ds,
hist=None,
finalize=True,
revert_to_original=True,
wsbin_width=1,
n_wsbins=40,
n_sectors=12,
percentile=0.5,
wv_count=None,
landmask=None,
):
"""Add timeseries to existing histogram
Uses the stability parameters to derive
summed histogram of the surface layer temperature scale,
squared surface layer temperature scale, boundary
layer height and logarithm of roughness length. The time
dimension is copied so that one can see what was the
period that was used to generate the histogram.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
ds: xr.Dataset
Dataset containing variables
['wind_speed', 'wind_direction', 'LANDMASK', 'ZNT']
Can have any pywasp spatial structure
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
finalize: boolean
Convert the stabilility histogram into a mean and
standard deviation?
revert_to_original: boolean
Return the histogram input in original format from the
input `ds`? If false, keeps the data in stacked point
format which is more efficient when doing large
calculations for numerical wind atlases.
wsbin_width: float
width of wind speed bins
n_wsbins: int
Number of wind speed bins
n_sectors: int
Number of sectors (wind direction bins)
kwargs:
Other keyword arguments passed on to _finalize_meso
Returns
-------
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
with the values from ds-timeseries added.
"""
if hist is None:
hist = _spatial_stack(ds, copy=False)[["point"]]
hist["z0sum"] = _init_histogram(
ds[["ZNT"]], "z0sum", n_sectors, n_wsbins, wsbin_width
)
hist["count_bins_exceeded"] = xr.zeros_like(hist["point"], dtype=np.int32)
ds_sub = ds[[WS, WD, "ZNT"]]
dss = _spatial_stack(ds_sub, copy=False).transpose(..., DIM_TIME, "point")
hist = _spatial_stack(hist, copy=False).transpose(..., "wsbin", "sector", "point")
z0sum, count_bins_exceeded = _create_hist(
dss[WS],
dss[WD],
0,
hist.wsbin,
hist.sector,
np.log(dss["ZNT"]),
hist["z0sum"].values,
hist["count_bins_exceeded"].values,
)
hist["z0sum"] = (("wsbin", "sector", "point"), z0sum)
hist["count_bins_exceeded"] = (("point"), count_bins_exceeded)
dic_attrs = _create_time_attributes(ds_sub, hist=hist)
hist = hist.assign_attrs(dic_attrs).drop_vars("point", errors="ignore")
if finalize:
if wv_count is None:
wv_count = bwc_from_tswc(
dss,
normalize=False,
revert_to_original=False,
wsbin_width=wsbin_width,
n_wsbins=n_wsbins,
n_sectors=n_sectors,
)
else:
wv_count = bwc_from_tswc(
dss, hist=wv_count, normalize=False, revert_to_original=False
)
res = _z0_mean(hist, wv_count["wv_count"], landmask, percentile=percentile)
return _update_history(res)
if revert_to_original:
res = _spatial_unstack(hist)
return _update_history(res)
else:
return _update_history(hist)
def _z0_mean(hist, wv_count, landmask, percentile=0.5):
"""Compute mean or standard deviation from histogram
Compute mean or standard deviation from histogram. The
variables are binned by sector and wind speed bin and
a mean or std. dev. is calulcated based on the highest
percentile of wind speeds for each sector.
Parameters
----------
hist: xr.Dataset
Dataset with variables `sum_temp_scale`, `sum_squared_temp_scale`
and `sum_pblh` and dimensions `wsbin`, `sector`
Returns
-------
clim: xr.Dataset
Dataset with mean mesoclimate variables per sector
"""
hist["wv_count"] = wv_count
hist["landmask"] = landmask
hist = _spatial_stack(hist, copy=False)
cond_hist, wd_count = _COND_HIST(hist["z0sum"], hist["wv_count"], percentile)
z0mean = _MEAN_HIST(cond_hist, wd_count)
hist["logz0sec"] = (("sector", "point"), z0mean)
hist["wd_count"] = (("sector", "point"), wd_count)
hist["mean_z0"] = mean_z0(hist)
hist["mean_z0"] = xr.where(hist["landmask"] == 0, 0.0, hist["mean_z0"])
return _spatial_unstack(hist)
def mean_z0(hist):
"""Compute mean roughness
Compute mean roughness from log(z0) per wsbin and
sector
.. warning::
This function is experimental and its signature may change.
Parameters
----------
hist: xr.Dataset
Dataset with variables `logz0sec` and `wd_count`
and dimensions `wsbin`, `sector`
Returns
-------
ds: xr.DataArray
Dataarray with mean z0 per sector
"""
mean_z0 = hist["logz0sec"] * hist["wd_count"] / hist["wd_count"].sum(dim=["sector"])
return np.exp(mean_z0.sum(dim=["sector"]))
def _normalize_hist(hist_point):
"""Convert frequency to probability-density histogram
Convert frequency to probability density histogram in the
WAsP core structure, i.e. with an additional variables
containing wind speed sector frequencies.
Parameters
----------
hist: xr.Dataset
Dataset with the frequency of occurence in each bin
Returns
-------
hist: xr.Dataset
Binned wind climate with probability density in each bin
"""
wsfreq, wdfreq = _NORM_HIST(hist_point["wv_count"])
var_names = ("wsfreq", "wdfreq")
var_dims = (("wsbin", "sector", "point"), ("sector", "point"))
for name, dims in zip(var_names, var_dims):
hist_point[name] = (dims, eval(name))
hist = hist_point.drop_vars(["wv_count", "count_bins_exceeded"], errors="ignore")
return hist
def _finalize_meso(total_stab_histogram, wv_count, landmask=None, percentile=0.5):
"""Finalize cumulative mesoscale climate
Computes conditional mean mesoscale climate from cumulative
histogram with mesoscale climate variables. The mean
and standard deviation of heat flux are calculated based on
the values occuring at the 50 percent highest wind speeds
(percentile=0.5). These fields are calculated seperately
over land and sea.
Parameters
----------
total_stab_histogram: xr.Dataset
Dataset with the sum of mesoscale variables
wrf: xr.Dataset
Dataset (e.g. `wrfout`) with the WRF variables
`LANDMASK`, `HGT`, `LU_INDEX`, `COSALPHA`,
`SINALPHA`
Returns
-------
ds: xr.Dataset
Dataset with mean or standard deviation of the
histograms
"""
hist = _stab_mean_std(total_stab_histogram, wv_count, percentile=percentile)
if landmask is not None:
hist = _nearest_landsea(hist, landmask)
_update_var_attrs(hist, _MEAN_STAB_ATTRS)
return hist
def _nearest_landsea(ds, landmask):
"""Calculate mean values over sea and land
To use stability setting over land and sea seperately
we calculate two fields from one: the land and sea fields
are generated using Scipy's nearest neighbour interpolation
Parameters
----------
ds: xr.Dataset
Dataset with the variables given in `ori_vars`
Returns
-------
ds: xr.Dataset
Dataset with the variables given in `new_vars`
"""
ori_vars = [
"mean_temp_scale",
"rms_temp_scale",
"mean_pblh_scale",
"mean_temp_scale",
"rms_temp_scale",
"mean_pblh_scale",
]
new_vars = [
"mean_temp_scale_land",
"rms_temp_scale_land",
"mean_pblh_scale_land",
"mean_temp_scale_sea",
"rms_temp_scale_sea",
"mean_pblh_scale_sea",
]
for var, ori_var in zip(new_vars, ori_vars):
list_secs = []
for sec in ds.sector.values:
# era5 has a fraction for the landmask, if the ratio
# land sea is larger than 0.5 we consider it land
if "_land" in var:
landflag = 1
else:
landflag = 0
tarr = (
ds[ori_var]
.squeeze()
.sel(sector=sec)
.where(np.round(landmask) == landflag)
)
array = np.ma.masked_invalid(tarr)
xx, yy = np.meshgrid(ds.west_east.values, ds.south_north.values)
x1 = xx[~array.mask]
y1 = yy[~array.mask]
newarr = array[~array.mask]
GD1 = interpolate.griddata(
(x1, y1), newarr.ravel(), (xx, yy), method="nearest"
)
list_secs.append(GD1)
ds[var] = (("sector", "south_north", "west_east"), list_secs)
return ds[new_vars]
def _init_histogram_baro(ds, wsbin_width=1, nws=40, nwd=12):
"""Initialize histogram
Initialize histogram with the variables needed for
a stability histogram. The dimensions of the spatial
structure are the same as the input ds and extended with
dimensions nws, nwd.
Parameters
----------
ds: xr.Dataset
Can have any of the pywasp spatial data structures
var_dic: dict
Dictonairy contain the variables that will be added to
the xarray container, including meta-data according to
the CF conventions
wsbin: float
width of wind speed bins
wdsec: float
width of wind direction bins
nws: int
Number of wind speed bins
nwd: int
Number of wind direction bins
Returns
-------
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
"""
bin_centers_ws = create_wsbin_coords(bins=nws, width=wsbin_width)
bin_centers_wd = create_sector_coords(nwd)
hist = _spatial_stack(ds, copy=False)[["point"]]
# stabdim = hist.sizes["point"]
# sum_dugdz = np.zeros((nws, nwd, stabdim), order="F", dtype=np.float32)
# sum_dvgdz = np.zeros((nws, nwd, stabdim), order="F", dtype=np.float32)
# count_bins_exceeded = np.zeros(stabdim, order="F", dtype=np.int32)
hist = hist.assign_coords({**bin_centers_wd.coords, **bin_centers_ws.coords})
var_names = (
"sum_dugdz",
"sum_dvgdz",
"count_bins_exceeded",
)
var_dims = (("wsbin", "sector", "point"),) * 2 + ("point",)
for name, dims in zip(var_names, var_dims):
hist[name] = (dims, eval(name), _HIS_BARO_ATTRS[name])
return hist
def _finalize_baro(hist, percentile=0.5):
"""Compute mean or standard deviation from histogram
Compute mean or standard deviation from histogram. The
variables are binned by sector and wind speed bin and
a mean or std. dev. is calulcated based on the highest
percentile of wind speeds for each sector.
Parameters
----------
hist: xr.Dataset
Dataset with variables `tstarsum`, `tstar2sum`
and `pblhsum` and dimensions `wsbin`, `sector`
Returns
-------
clim: xr.Dataset
Dataset with mean mesoclimate variables per sector
"""
hist = _spatial_stack(hist, copy=False)
cond_hist, freq = _COND_HIST(hist["sum_dugdz"], hist["wv_count"], percentile)
mean_dugdz = _MEAN_HIST(cond_hist, freq)
cond_hist, freq = _COND_HIST(hist["sum_dvgdz"], hist["wv_count"], percentile)
mean_dvgdz = _MEAN_HIST(cond_hist, freq)
clim = hist[["point", "sector"]]
clim["mean_dgdz_dir_north"] = (
("sector", "point"),
np.arctan2(-mean_dvgdz, -mean_dugdz) * RTD,
)
clim["mean_dgdz_dir_south"] = (
("sector", "point"),
np.arctan2(mean_dvgdz, mean_dugdz) * RTD,
)
clim["mean_dgdz_dir"] = xr.where(
clim["south_north"] < 0,
clim["mean_dgdz_dir_south"],
clim["mean_dgdz_dir_north"],
)
clim["mean_dgdz"] = (
("sector", "point"),
np.sqrt(mean_dugdz**2 + mean_dvgdz**2),
)
out = _spatial_unstack(
clim.drop_vars(["mean_dgdz_dir_north", "mean_dgdz_dir_south"])
)
_update_var_attrs(out, _MEAN_BARO_ATTRS)
return out
def baroclinicity_histogram(
ds, hist=None, finalize=True, revert_to_original=True, **kw
):
"""Add timeseries to existing histogram
Uses the baroclinicity parameters to derive
summed histogram of the geostrophic shear components.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
ds: xr.Dataset
Dataset containing variables
['wind_speed', 'wind_direction', 'PSFC', 'T2', 'UST',
'HFX', 'LH', 'PBLH', 'LANDMASK', 'ZNT']
Can have any pywasp spatial structure
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
finalize: boolean
Convert the stabilility histogram into a mean and
standard deviation?
revert_to_original: boolean
Return the histogram input in original format from the
input `ds`? If false, keeps the data in stacked point
format which is more efficient when doing large
calculations for numerical wind atlases.
**kwargs
Arbitrary keyword arguments passed on to
initialization of histogram
Returns
-------
hist: xr.Dataset
Histogram with dimensions point, wsbin, sector
with the values from ds-timeseries added.
"""
if hist is None:
hist = _init_histogram_baro(ds, **kw)
ds_sub = ds[[WS, WD, "dugdz", "dvgdz"]]
dss = _spatial_stack(ds_sub, copy=False).transpose(..., DIM_TIME, "point")
# add = np.full(dss[WS].shape, 1.0, order="F", dtype=np.float32)
hist = _spatial_stack(hist, copy=False).transpose(..., "wsbin", "sector", "point")
# error arrays are inout in fortran, so we must create copies
# to not add the error counter 2 times
ierr_sum_dvgdz = hist["count_bins_exceeded"].copy().values
sum_dugdz, count_bins_exceeded = _create_hist(
dss[WS],
dss[WD],
0,
hist.wsbin,
hist.sector,
dss["dugdz"],
hist["sum_dugdz"].values,
hist["count_bins_exceeded"].values,
)
sum_dvgdz, count_bins_exceeded = _create_hist(
dss[WS],
dss[WD],
0,
hist.wsbin,
hist.sector,
dss["dvgdz"],
hist["sum_dvgdz"].values,
ierr_sum_dvgdz,
)
var_names = (
"sum_dugdz",
"sum_dvgdz",
"count_bins_exceeded",
)
for name in var_names:
hist[name] = (hist[name].dims, eval(name))
_update_var_attrs(hist[name], _HIS_BARO_ATTRS)
hist[name].attrs["_pwio_data_is_2d"] = False
dic_attrs = _create_time_attributes(ds_sub, hist=hist)
hist = hist.assign_attrs(dic_attrs)
if finalize:
return _finalize_baro(hist)
if revert_to_original:
return _spatial_unstack(hist)
else:
return hist
[docs]
@_validate_bwc_wrapper_factory(run_extra_checks=False)
def bwc_resample_like(source, target):
"""
Resamples a histogram to different sector or wind speed bin structure.
Parameters
----------
source: xarray.Dataset
Binned wind climate xr.Dataset object that we want to resample
target: xarray.Dataset
Binned wind climate xr.Dataset object structure that we want to resample to
Returns
-------
xarray.Dataset
Binned Wind Climate
"""
result = bwc_resample_sectors(source, target.sizes["sector"])
result = bwc_resample_wsbins_like(result, target)
return result
[docs]
@_validate_bwc_wrapper_factory(run_extra_checks=False)
def bwc_resample_sectors(source, n_sectors=12, offset=None):
"""
Resamples a histogram to different sector structure.
Parameters
----------
source: xarray.Dataset
Binned wind climate xr.Dataset object that we want to resample
target: xarray.Dataset
Binned wind climate xr.Dataset object structure that we want to resample to
offset: xarray.DataArray or float
Offset in degrees with which the histogram are rotated (positive
for clockwise rotation). If it is a DataArray, needs to have the same spatial shape as
the binned wind climate.
Returns
-------
xarray.Dataset
Binned Wind Climate
"""
if source.sizes["sector"] == n_sectors and offset is None:
return source
source = _spatial_stack(source)
# copy of old attributes so we can unstack to the original structure
oldattrs = {key: val for key, val in source.attrs.items() if "_pwio" in key}
if offset is None:
offset = xr.full_like(source.point, 0.0)
elif isinstance(offset, float | int):
offset = xr.full_like(source.point, offset)
else:
offset = to_point(offset)
def _resample(wv_freq, doff, n_sectors):
return _RESAMPLE_HIST(wv_freq, doff, n_sectors)
wv_freq = xr.apply_ufunc(
_resample,
(source.wsfreq * source.wdfreq).transpose(..., "wsbin", "sector"),
offset,
input_core_dims=[["wsbin", "sector"], []],
output_core_dims=[["wsbin", "sector"]],
exclude_dims={"sector"},
dask="allowed",
vectorize=True,
kwargs={"n_sectors": n_sectors},
keep_attrs=True,
)
wdcenters = create_sector_coords(n_sectors)
wv_freq = wv_freq.assign_coords({**wdcenters.coords})
source = _bwc_from_counts(wv_freq)
# Update metadata
source.attrs["title"] = "Binned wind climate"
# apply old attributes so that we can unstack
source = source.assign_attrs(oldattrs)
res = _spatial_unstack(source)
return _update_history(res)
[docs]
@_validate_bwc_wrapper_factory(run_extra_checks=False)
def bwc_resample_wsbins_like(source, target, fit_weibull=True, errors="raise"):
"""
Resamples a histogram to different wind speed bin structure.
Parameters
----------
source: xarray.Dataset
Binned wind climate xr.Dataset object that we want to resample
target: xarray.Dataset
Binned wind climate xr.Dataset object structure that we want to resample to
fit_weibull: boolean, optional
Default True, when resampling the histogram find the sectorwise A and k and
use these to resample the histogram. If false, use the A and k provided in
the dataset instead.
errors: {'raise', 'ignore'}, optional
If 'raise' (default), raises a ValueError error if probability mass
is transfered beyond the highest wind speed bin.
If "ignore", all points including those that lose probability mass
are returned.
Returns
-------
xarray.Dataset
Binned Wind Climate
"""
# we reset coordinates with drop=True to make sure there is no user added coordinates that cause the dataarrays to be different
if (
not source["sector"]
.reset_coords(drop=True)
.equals(target["sector"].reset_coords(drop=True))
):
raise ValueError(
"Input and target histogram must have the same sector structure!"
)
# source and target wsbin structure are the same, so no resampling needed
same_wsbins = (
source["wsbin"]
.reset_coords(drop=True)
.equals(target["wsbin"].reset_coords(drop=True))
)
if same_wsbins and fit_weibull:
return source
source = _spatial_stack(source).transpose(..., "wsbin", "sector")
target = _spatial_stack(target).transpose(..., "wsbin", "sector")
if not fit_weibull:
if not all(v in target.data_vars for v in ["A", "k"]):
raise ValueError(
"Variables A and k are needed for resampling when fit_weibull is set to False"
)
if fit_weibull:
wsfreq, ierror = xr.apply_ufunc(
_RESAMPLE_HIST_WSBINS,
source.wsfreq,
source.wsceil,
target.wsceil,
target.wsfloor,
input_core_dims=[["wsbin", "sector"], *[["wsbin"]] * 3],
output_core_dims=[["wsbin", "sector"], ["sector"]],
exclude_dims={"wsbin"},
dask="allowed",
vectorize=True,
keep_attrs=True,
)
else:
wsfreq, ierror = xr.apply_ufunc(
_RESAMPLE_HIST_WSBINS_AK,
source.wsfreq,
source.wsceil,
target.A,
target.k,
target.wsceil,
target.wsfloor,
input_core_dims=[
["wsbin", "sector"],
["wsbin"],
*[["sector"]] * 2,
*[["wsbin"]] * 2,
],
output_core_dims=[["wsbin", "sector"], ["sector"]],
exclude_dims={"wsbin"},
dask="allowed",
vectorize=True,
keep_attrs=True,
)
if ierror.sum() > 0 and errors == "raise":
i_points = [i for i, v in enumerate(ierror) if v > 0]
str_points = ",".join(map(str, i_points))
raise ValueError(
f"Resampling of the wsbins in {len(i_points)} out of {source.sizes['point']} points (i={str_points}) exceeded the number of target bins!"
)
result, _ = xr.broadcast(source[["point", "sector"]], target[["wsbin"]])
result = result.drop_vars("point", errors="ignore")
result["wsfreq"] = wsfreq
result["wdfreq"] = source["wdfreq"]
result.coords["wsceil"] = target["wsceil"].reset_coords(drop=True)
result.coords["wsfloor"] = target["wsfloor"].reset_coords(drop=True)
res = _spatial_unstack(result)
return _update_history(res)
def _get_dgdz(z, sp, levs, lats):
"""
Calculate the gradients of the zonal and meridional geostrophic wind components with respect to height.
Parameters
----------
z : array_like
Geopotential height array.
sp : array_like
Surface pressure array.
levs : array_like
Pressure levels array.
lats : array_like
Latitudes array.
Returns
-------
tuple of array_like
Gradients of the zonal (dugdz) and meridional (dvgdz) geostrophic wind components with respect to height.
"""
dugdz, dvgdz = _DGDZ(z, sp, levs, lats)
return (dugdz, dvgdz)
def _get_dgdz_t(z, sp, levs, lats):
"""
Calculate the time-dependent gradients of the zonal and meridional geostrophic wind components with respect to height.
Parameters
----------
z : array_like
Geopotential height array.
sp : array_like
Surface pressure array.
levs : array_like
Pressure levels array.
lats : array_like
Latitudes array.
Returns
-------
tuple of array_like
Time-dependent gradients of the zonal (dugdz) and meridional (dvgdz) geostrophic wind components with respect to height.
"""
dugdz, dvgdz = _DGDZ_TIME(z, sp, levs, lats)
return (dugdz, dvgdz)
def _calc_dgdz(sfc):
"""
Calculate the vertical gradients of the zonal and meridional geostrophic wind components with respect to height for a given surface dataset.
Parameters
----------
sfc : xarray.Dataset
Surface dataset containing geopotential height, surface pressure, wind speed, and wind direction.
Returns
-------
xarray.Dataset
Dataset with added gradients of the zonal (dugdz) and meridional (dvgdz) geostrophic wind components with respect to height.
"""
dugdz, dvgdz = xr.apply_ufunc(
_get_dgdz,
sfc["z"]
.transpose("west_east", "south_north", "isobaricInhPa", "time")
.astype("float32"),
(sfc["sp"].transpose("west_east", "south_north", "time") / 100).astype(
"float32"
),
sfc["isobaricInhPa"],
sfc["south_north"],
input_core_dims=[
["west_east", "south_north", "isobaricInhPa"],
["west_east", "south_north"],
["isobaricInhPa"],
["south_north"],
],
output_core_dims=[["west_east", "south_north"], ["west_east", "south_north"]],
dask="allowed",
vectorize=True,
keep_attrs=True,
)
out = xr.Dataset(
{"dugdz": dugdz, "dvgdz": dvgdz},
)
out = set_crs(out, 4326)
return out.transpose("time", "south_north", "west_east")
def _calc_dgdz_t(sfc):
"""
Calculate vertical gradients of the zonal and meridional geostrophic wind components with respect to height for a given dataset.
Parameters
----------
sfc : xarray.Dataset
Dataset containing geopotential height, surface pressure, wind speed, and wind direction.
Returns
-------
xarray.Dataset
Dataset with added time-dependent gradients of the zonal (dugdz) and meridional (dvgdz) wind components with respect to height.
"""
dugdz, dvgdz = _get_dgdz_t(
sfc["z"].transpose("west_east", "south_north", "isobaricInhPa", "time"),
sfc["sp"].transpose("west_east", "south_north", "time") / 100,
sfc["isobaricInhPa"],
sfc["south_north"],
)
out = xr.Dataset(
{"dugdz": dugdz, "dvgdz": dvgdz},
)
out = set_crs(out, 4326)
return out.transpose("time", "south_north", "west_east")