"""Routines to calculate the spectral correction factor for extreme winds"""
import numpy as np
from windkit.spatial import spatial_stack, spatial_unstack
import xarray as xr
from windkit.spatial import create_dataset, get_crs
from windkit.metadata import _SPECTRUM_ATTRS, update_history
from ..core import Rvea0292
get_spectrum = Rvea0292.tstospectrumpoints
get_corr_fac = Rvea0292.getspcffromextrapolation
get_corr_fac_auto = Rvea0292.getspcffromextrapolationauto
def _calc_spectrum_len(len_ts, num_pd=25):
"""
Calculates spectral length
Parameters
----------
len_ts : int
Length of the input timeseries
num_pd : int
Number of periods to smooth
Returns
-------
len_smooth : int
Spectral length
"""
len_trim_ts = int(2 ** np.floor(np.log(len_ts) / np.log(2)))
len_spectrum = int(len_trim_ts / 2 + 1)
smooth_max = int(num_pd * np.log10(len_spectrum) + 1)
len_smooth = 0
nlo = 0
nhi = 0
for i in np.arange(smooth_max) + 1:
nlo = nhi + 1
nhi = int(10 ** (i / num_pd))
if nlo <= nhi:
len_smooth += 1
return len_smooth
[docs]
def get_spec_corr_fac(wind_ts, fc=0.8, fh=72, f01=0.6, f02=0.9, auto=False):
"""Calculate spectral correction factor for model time-series
.. warning::
This function is experimental and its signature may change.
Parameters
----------
wind_ts : xarray.Dataset
Dataset with wind_speed variable that has (time, height, south_north, and
west_east) dimensions.
fc : float
Frequency where to shift from the model spectrum to the extrapolated spectrum
fh : float
Highest sample frequency of the extrapolated spectrum #/day (Default: 72 = 20minutes)
f01 : float
Frequency used to fit a linear relationship to determine slope of the model spectrum
f02 : float
Frequency used to fit a linear relationship to determine slope of the model spectrum
auto : bool
Run the automatic extrapolation identification script? If true, this ignores the
fc, f01, and f02 parameters, instead using pre-set values.
Returns
-------
out_ds : xarray.Dataset
PyWAsP xr.Dataset containing spectral correction factors.
Notes
-----
This routine fits a spectrum to model time-series, and then uses the 5/3
extrapolation to create a hybrid spectrum that is used to determine the
adjustment of the V50 winds due to the missing spectral information.
"""
wind_ds_pt = spatial_stack(wind_ts)
wind_ts_pt = wind_ds_pt.wind_speed
# Xarray returns the time difference in nano-seconds, we want # per day
model_sample_freq = (24 * 3600 * 1e9) / float(
wind_ts_pt.time[1] - wind_ts_pt.time[0]
)
# Fit the spectrum
len_spec = _calc_spectrum_len(wind_ts.sizes["time"])
spe_freq, spe_pow = get_spectrum(wind_ts_pt, model_sample_freq, len_spec)
# Calculate the correction factors
wind_means = wind_ts_pt.groupby("time.year").mean("time")
wind_sdevs = wind_ts_pt.groupby("time.year").std("time")
if auto:
result = get_corr_fac_auto(
spe_freq, spe_pow, model_sample_freq, fh, wind_means, wind_sdevs
)
else:
result = get_corr_fac(
spe_freq,
spe_pow,
model_sample_freq,
fc,
fh,
f01,
f02,
wind_means,
wind_sdevs,
)
result = (result,)
##########################
# Build the output dataset
##########################
result_names = ("spec_corr_fac",)
# Combine the results into a dataset
out_ds = create_dataset(
wind_ds_pt["west_east"],
wind_ds_pt["south_north"],
wind_ds_pt["height"],
get_crs(wind_ds_pt),
).drop_vars("output")
for i, name in enumerate(result_names):
out_ds[name] = xr.DataArray(result[i], dims=["point"])
out_ds[name].attrs = _SPECTRUM_ATTRS[name]
# Add attributes from gewc
out_ds.attrs = wind_ds_pt.attrs
out_ds.attrs["title"] = "WEng Spectral correction factor"
# Convert to original projection and raster if necessary
out_ds = spatial_unstack(out_ds)
return update_history(out_ds)