Source code for pywasp.lincom.spectral_correction_method

"""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)