From d7b07663b626faa9b75b38c2b8000fc085ae1f98 Mon Sep 17 00:00:00 2001 From: akeeste Date: Tue, 21 Nov 2023 15:19:10 -0600 Subject: [PATCH 01/13] loads/general --- mhkit/loads/general.py | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/mhkit/loads/general.py b/mhkit/loads/general.py index 28558eb58..03c5f2d14 100644 --- a/mhkit/loads/general.py +++ b/mhkit/loads/general.py @@ -1,5 +1,6 @@ from scipy.stats import binned_statistic -import pandas as pd +import pandas as pd +import xarray as xr import numpy as np import fatpack @@ -10,7 +11,7 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]): Parameters ----------- - data : pandas DataFrame + data : pandas DataFrame or xarray Dataset Time-series statistics of data signal(s) bin_against : array Data signal to bin data against (e.g. wind speed) @@ -27,9 +28,9 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]): Standard deviation of each bim """ - if not isinstance(data, pd.DataFrame): + if not isinstance(data, (pd.DataFrame, xr.Dataset)): raise TypeError( - f'data must be of type pd.DataFrame. Got: {type(data)}') + f'data must be of type pd.DataFrame or xr.Dataset. Got: {type(data)}') try: bin_against = np.asarray(bin_against) except: @@ -40,10 +41,14 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]): except: raise TypeError( f'bin_edges must be of type np.ndarray. Got: {type(bin_edges)}') + + # If input is pandas, convert to xarray + if isinstance(data,pd.DataFrame): + data = data.to_xarray() # Determine variables to analyze if len(data_signal)==0: # if not specified, bin all variables - data_signal=data.columns.values + data_signal = list(data.keys()) else: if not isinstance(data_signal, list): raise TypeError( @@ -60,12 +65,13 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]): statistic='mean',bins=bin_edges) # Calculate std of bins std = [] - stdev = pd.DataFrame(data[signal_name]) - stdev.set_index(bin_stat.binnumber,inplace=True) + stdev = xr.DataArray(data = data[signal_name].values, + dims = 'binnumber', + coords = {'binnumber': bin_stat.binnumber}) for i in range(1,len(bin_stat.bin_edges)): try: - temp = stdev.loc[i].std(ddof=0) - std.append(temp[0]) + temp = stdev.sel(binnumber=i).std().values + std.append(np.array(temp,ndmin=1)[0]) except: std.append(np.nan) bin_stat_list.append(bin_stat.statistic) From 564010d463d18dfecf4140f65754d1a189753425 Mon Sep 17 00:00:00 2001 From: akeeste Date: Mon, 27 Nov 2023 16:02:17 -0600 Subject: [PATCH 02/13] loads/extreme conversion to xarray --- mhkit/loads/extreme.py | 127 ++++++++++++++++++++++++++++------------- mhkit/loads/general.py | 40 ++++++++----- 2 files changed, 113 insertions(+), 54 deletions(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index 4cff02ca0..d62a31a15 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -1,5 +1,6 @@ import numpy as np import pandas as pd +import xarray as xr from scipy import stats, optimize, signal from mhkit.wave.resource import frequency_moment from mhkit.utils import upcrossing, custom @@ -38,6 +39,7 @@ def _peaks_over_threshold(peaks, threshold, sampling_rate): return independent_storm_peaks + def global_peaks(t, data): """ Find the global peaks of a zero-centered response time-series. @@ -609,7 +611,7 @@ def _cdf(self, x): return _LongTermExtreme(name="long_term_extreme", weights=weights, ste=ste) -def mler_coefficients(rao, wave_spectrum, response_desired): +def mler_coefficients(rao, wave_spectrum, response_desired, to_pandas=True): """ Calculate MLER (most likely extreme response) coefficients from a sea state spectrum and a response RAO. @@ -618,11 +620,13 @@ def mler_coefficients(rao, wave_spectrum, response_desired): ---------- rao: numpy ndarray Response amplitude operator. - wave_spectrum: pd.DataFrame + wave_spectrum: pd.DataFrame or xarray Dataset Wave spectral density [m^2/Hz] indexed by frequency [Hz]. response_desired: int or float Desired response, units should correspond to a motion RAO or units of force for a force RAO. + to_pandas: bool (optional) + Flag to output pandas instead of xarray. Default = True. Returns ------- @@ -638,58 +642,66 @@ def mler_coefficients(rao, wave_spectrum, response_desired): if not isinstance(rao, np.ndarray): raise TypeError( f'rao must be of type np.ndarray. Got: {type(rao)}') - if not isinstance(wave_spectrum, pd.DataFrame): + if not isinstance(wave_spectrum, (pd.DataFrame, xr.Dataset)): raise TypeError( - f'wave_spectrum must be of type pd.DataFrame. Got: {type(wave_spectrum)}') + f'wave_spectrum must be of type pd.DataFrame or xr.Dataset. Got: {type(wave_spectrum)}') if not isinstance(response_desired, (int, float)): raise TypeError( f'response_desired must be of type int or float. Got: {type(response_desired)}') - - freq_hz = wave_spectrum.index.values + if not isinstance(to_pandas, bool): + raise TypeError( + f'to_pandas must be of type bool. Got: {type(to_pandas)}') + + # If input is pandas, convert to xarray + if isinstance(wave_spectrum,pd.DataFrame): + wave_spectrum = wave_spectrum.to_xarray() + # convert from Hz to rad/s + freq_hz = wave_spectrum.index.values freq = freq_hz * (2*np.pi) - # change from Hz to rad/s - wave_spectrum = wave_spectrum.iloc[:, 0].values / (2*np.pi) + key_name = list(wave_spectrum.keys())[0] + wave_spectrum = wave_spectrum[key_name].values / (2*np.pi) + # get delta dw = (2*np.pi - 0.) / (len(freq)-1) - + spectrum_r = np.zeros(len(freq)) # [(response units)^2-s/rad] _s = np.zeros(len(freq)) # [m^2-s/rad] _a = np.zeros(len(freq)) # [m^2-s/rad] _coeff_a_rn = np.zeros(len(freq)) # [1/(response units)] _phase = np.zeros(len(freq)) - + # Note: waves.A is "S" in Quon2016; 'waves' naming convention # matches WEC-Sim conventions (EWQ) # Response spectrum [(response units)^2-s/rad] -- Quon2016 Eqn. 3 spectrum_r[:] = np.abs(rao)**2 * (2*wave_spectrum) - + # calculate spectral moments and other important spectral values. m0 = (frequency_moment(pd.Series(spectrum_r, index=freq), 0)).iloc[0, 0] m1 = (frequency_moment(pd.Series(spectrum_r, index=freq), 1)).iloc[0, 0] m2 = (frequency_moment(pd.Series(spectrum_r, index=freq), 2)).iloc[0, 0] wBar = m1 / m0 - + # calculate coefficient A_{R,n} [(response units)^-1] -- Quon2016 Eqn. 8 # Drummen version. Dietz has negative of this. _coeff_a_rn[:] = np.abs(rao) * np.sqrt(2*wave_spectrum*dw) * \ ((m2 - freq*m1) + wBar*(freq*m0 - m1)) / (m0*m2 - m1**2) - + # save the new spectral info to pass out # Phase delay should be a positive number in this convention (AP) _phase[:] = -np.unwrap(np.angle(rao)) - + # for negative values of Amp, shift phase by pi and flip sign # for negative amplitudes, add a pi phase shift, then flip sign on # negative Amplitudes _phase[_coeff_a_rn < 0] -= np.pi _coeff_a_rn[_coeff_a_rn < 0] *= -1 - + # calculate the conditioned spectrum [m^2-s/rad] _s[:] = wave_spectrum * _coeff_a_rn[:]**2 * response_desired**2 _a[:] = 2*wave_spectrum * _coeff_a_rn[:]**2 * \ response_desired**2 - + # if the response amplitude we ask for is negative, we will add # a pi phase shift to the phase information. This is because # the sign of self.desiredRespAmp is lost in the squaring above. @@ -698,10 +710,21 @@ def mler_coefficients(rao, wave_spectrum, response_desired): # new spectral information, S. (AP) if response_desired < 0: _phase += np.pi - + mler = pd.DataFrame( data={'WaveSpectrum': _s, 'Phase': _phase}, index=freq_hz) mler = mler.fillna(0) + + + mler = xr.Dataset(data_vars = {'WaveSpectrum': (['frequency'], _s), + 'Phase': (['frequency'], _phase)}, + coords = {'frequency': freq_hz} + ) + mler.fillna(0) + + if to_pandas: + mler = mler.to_pandas() + return mler @@ -715,7 +738,7 @@ def mler_simulation(parameters=None): Parameters ---------- - parameters: dict (optional) + parameters: dict (fto_) Simulation parameters. Keys: ----- @@ -764,7 +787,7 @@ def mler_simulation(parameters=None): return sim -def mler_wave_amp_normalize(wave_amp, mler, sim, k): +def mler_wave_amp_normalize(wave_amp, mler, sim, k, to_pandas=True): """ Function that renormalizes the incoming amplitude of the MLER wave to the desired peak height (peak to MSL). @@ -773,13 +796,15 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k): ---------- wave_amp: float Desired wave amplitude (peak to MSL). - mler: pd.DataFrame + mler: pandas DataFrame or xarray Dataset MLER coefficients generated by 'mler_coefficients' function. sim: dict Simulation parameters formatted by output from 'mler_simulation'. k: numpy ndarray - Wave number. + Wave number + to_pandas: bool (optional) + Flag to output pandas instead of xarray. Default = True. Returns ------- @@ -790,9 +815,9 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k): k = np.array(k) except: pass - if not isinstance(mler, pd.DataFrame): + if not isinstance(mler, (pd.DataFrame, xr.Dataset)): raise TypeError( - f'mler must be of type pd.DataFrame. Got: {type(mler)}') + f'mler must be of type pd.DataFrame or xr.Dataset. Got: {type(mler)}') if not isinstance(wave_amp, (int, float)): raise TypeError( f'wave_amp must be of type int or float. Got: {type(wave_amp)}') @@ -802,10 +827,17 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k): if not isinstance(k, np.ndarray): raise TypeError( f'k must be of type ndarray. Got: {type(k)}') + if not isinstance(to_pandas, bool): + raise TypeError( + f'to_pandas must be of type bool. Got: {type(to_pandas)}') - freq = mler.index.values * 2*np.pi + # If input is pandas, convert to xarray + if isinstance(mler,pd.DataFrame): + mler = mler.to_xarray() + + freq = mler.coords['frequency'].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta - + wave_amp_time = np.zeros((sim['maxIX'], sim['maxIT'])) for ix, x in enumerate(sim['X']): for it, t in enumerate(sim['T']): @@ -814,22 +846,24 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k): np.sqrt(2*mler['WaveSpectrum']*dw) * np.cos(freq*(t-sim['T0']) - k*(x-sim['X0']) + mler['Phase']) ) - + tmp_max_amp = np.max(np.abs(wave_amp_time)) - + # renormalization of wave amplitudes rescale_fact = np.abs(wave_amp) / np.abs(tmp_max_amp) + # rescale the wave spectral amplitude coefficients - spectrum = mler['WaveSpectrum'] * rescale_fact**2 - - mler_norm = pd.DataFrame(index=mler.index) - mler_norm['WaveSpectrum'] = spectrum - mler_norm['Phase'] = mler['Phase'] - + mler_norm = mler['WaveSpectrum'] * rescale_fact**2 + mler_norm = mler_norm.to_dataset() + mler_norm = mler_norm.assign({'Phase': ('frequency', mler['Phase'].data)}) + + if to_pandas: + mler_norm = mler_norm.to_pandas() + return mler_norm -def mler_export_time_series(rao, mler, sim, k): +def mler_export_time_series(rao, mler, sim, k, to_pandas=True): """ Generate the wave amplitude time series at X0 from the calculated MLER coefficients @@ -838,13 +872,15 @@ def mler_export_time_series(rao, mler, sim, k): ---------- rao: numpy ndarray Response amplitude operator. - mler: pd.DataFrame + mler: pandas DataFrame or xarray Dataset MLER coefficients dataframe generated from an MLER function. sim: dict Simulation parameters formatted by output from 'mler_simulation'. k: numpy ndarray Wave number. + to_pandas: bool (optional) + Flag to output pandas instead of xarray. Default = True. Returns ------- @@ -864,17 +900,24 @@ def mler_export_time_series(rao, mler, sim, k): if not isinstance(rao, np.ndarray): raise TypeError( f'rao must be of type ndarray. Got: {type(rao)}') - if not isinstance(mler, pd.DataFrame): + if not isinstance(mler, (pd.DataFrame, xr.Dataset)): raise TypeError( - f'mler must be of type pd.DataFrame. Got: {type(mler)}') + f'mler must be of type pd.DataFrame or xr.Dataset. Got: {type(mler)}') if not isinstance(sim, dict): raise TypeError( f'sim must be of type dict. Got: {type(sim)}') if not isinstance(k, np.ndarray): raise TypeError( f'k must be of type ndarray. Got: {type(k)}') + if not isinstance(to_pandas, bool): + raise TypeError( + f'to_pandas must be of type bool. Got: {type(to_pandas)}') - freq = mler.index.values * 2*np.pi # convert Hz to rad/s + # If input is pandas, convert to xarray + if isinstance(mler,pd.DataFrame): + mler = mler.to_xarray() + + freq = mler.coords['frequency'].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta # calculate the series @@ -892,8 +935,12 @@ def mler_export_time_series(rao, mler, sim, k): np.cos(freq*(ti-sim['T0']) - k*(xi-sim['X0'])) ) - mler_ts = pd.DataFrame(wave_amp_time, index=sim['T']) - mler_ts = mler_ts.rename(columns={0: 'WaveHeight', 1: 'LinearResponse'}) + mler_ts = xr.Dataset(data_vars = {'WaveHeight': (['time'], wave_amp_time[:,0]), + 'LinearResponse': (['time'], wave_amp_time[:,1])}, + coords = {'time': sim['T']}) + + if to_pandas: + mler_ts = mler_ts.to_pandas() return mler_ts diff --git a/mhkit/loads/general.py b/mhkit/loads/general.py index 03c5f2d14..514c09ab4 100644 --- a/mhkit/loads/general.py +++ b/mhkit/loads/general.py @@ -4,7 +4,7 @@ import numpy as np import fatpack -def bin_statistics(data,bin_against,bin_edges,data_signal=[]): +def bin_statistics(data,bin_against,bin_edges,data_signal=[],to_pandas=True): """ Bins calculated statistics against data signal (or channel) according to IEC TS 62600-3:2020 ED1. @@ -19,6 +19,8 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]): Bin edges with consistent step size data_signal : list, optional List of data signal(s) to bin, default = all data signals + to_pandas: bool (optional) + Flag to output pandas instead of xarray. Default = True. Returns -------- @@ -41,6 +43,9 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]): except: raise TypeError( f'bin_edges must be of type np.ndarray. Got: {type(bin_edges)}') + if not isinstance(to_pandas, bool): + raise TypeError( + f'to_pandas must be of type bool. Got: {type(to_pandas)}') # If input is pandas, convert to xarray if isinstance(data,pd.DataFrame): @@ -53,11 +58,11 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]): if not isinstance(data_signal, list): raise TypeError( f'data_signal must be of type list. Got: {type(data_signal)}') - - # Pre-allocate list variables - bin_stat_list = [] - bin_std_list = [] - + + # Pre-allocate variable dictionaries + bin_stat_list = {} + bin_std_list = {} + # loop through data_signal and get binned means for signal_name in data_signal: # Bin data @@ -74,17 +79,24 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]): std.append(np.array(temp,ndmin=1)[0]) except: std.append(np.nan) - bin_stat_list.append(bin_stat.statistic) - bin_std_list.append(std) - - # Convert to DataFrames - bin_mean = pd.DataFrame(np.transpose(bin_stat_list),columns=data_signal) - bin_std = pd.DataFrame(np.transpose(bin_std_list),columns=data_signal) - + + bin_stat_list[signal_name] = ('index', bin_stat.statistic) + bin_std_list[signal_name] = ('index', bin_stat.statistic) + + # Convert to Datasets + bin_mean = xr.Dataset(data_vars = bin_stat_list, + coords = {'index':np.arange(0,len(bin_stat.statistic))}) + bin_std = xr.Dataset(data_vars = bin_std_list, + coords = {'index':np.arange(0,len(bin_stat.statistic))}) + # Check for nans if bin_mean.isna().any().any(): print('Warning: some bins may be empty!') - + + if to_pandas: + bin_mean = bin_mean.to_pandas() + bin_std = bin_std.to_pandas() + return bin_mean, bin_std From ad2dcd1ddd1d38db2d1d9c05c24bbbd999828556 Mon Sep 17 00:00:00 2001 From: akeeste Date: Wed, 29 Nov 2023 15:43:51 -0600 Subject: [PATCH 03/13] test and bugfix for loads/general/bin_statistics --- mhkit/loads/general.py | 6 ++++-- mhkit/tests/loads/test_loads.py | 16 ++++++++++++++-- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/mhkit/loads/general.py b/mhkit/loads/general.py index 514c09ab4..db670a855 100644 --- a/mhkit/loads/general.py +++ b/mhkit/loads/general.py @@ -90,8 +90,10 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[],to_pandas=True): coords = {'index':np.arange(0,len(bin_stat.statistic))}) # Check for nans - if bin_mean.isna().any().any(): - print('Warning: some bins may be empty!') + for variable in list(bin_mean.variables): + if bin_mean[variable].isnull().any(): + print('Warning: bins for some variables may be empty!') + break if to_pandas: bin_mean = bin_mean.to_pandas() diff --git a/mhkit/tests/loads/test_loads.py b/mhkit/tests/loads/test_loads.py index 5befa99d9..61356ccb9 100644 --- a/mhkit/tests/loads/test_loads.py +++ b/mhkit/tests/loads/test_loads.py @@ -2,7 +2,6 @@ from numpy.testing import assert_array_almost_equal, assert_allclose from pandas._testing.asserters import assert_series_equal from pandas.testing import assert_frame_equal -from mhkit import utils from mhkit.wave import resource import mhkit.loads as loads import pandas as pd @@ -46,7 +45,20 @@ def test_bin_statistics(self): bin_edges = np.arange(3,26,1) # Apply function to calculate means - load_means =self.data['means'] + load_means = self.data['means'] + bin_against = load_means['uWind_80m'] + [b_means, b_means_std] = loads.general.bin_statistics(load_means, bin_against, bin_edges) + + assert_frame_equal(self.data['bin_means'],b_means) + assert_frame_equal(self.data['bin_means_std'],b_means_std) + + def test_bin_statistics_xarray(self): + # create array containg wind speeds to use as bin edges + bin_edges = np.arange(3,26,1) + + # Apply function to calculate means + load_means = self.data['means'] + load_means = load_means.to_xarray() bin_against = load_means['uWind_80m'] [b_means, b_means_std] = loads.general.bin_statistics(load_means, bin_against, bin_edges) From 03537bdec796fe6945ec4e56df894184d4611952 Mon Sep 17 00:00:00 2001 From: akeeste Date: Fri, 1 Dec 2023 12:19:31 -0600 Subject: [PATCH 04/13] fix bug in bin_statistics where std=0 --- mhkit/loads/extreme.py | 6 ++++-- mhkit/loads/general.py | 23 +++++++---------------- 2 files changed, 11 insertions(+), 18 deletions(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index d62a31a15..2441fe3cb 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -835,7 +835,8 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, to_pandas=True): if isinstance(mler,pd.DataFrame): mler = mler.to_xarray() - freq = mler.coords['frequency'].values * 2*np.pi + coord_names = list(mler.coords) + freq = mler.coords[coord_names[0]].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta wave_amp_time = np.zeros((sim['maxIX'], sim['maxIT'])) @@ -917,7 +918,8 @@ def mler_export_time_series(rao, mler, sim, k, to_pandas=True): if isinstance(mler,pd.DataFrame): mler = mler.to_xarray() - freq = mler.coords['frequency'].values * 2*np.pi + coord_names = list(mler.coords) + freq = mler.coords[coord_names[0]].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta # calculate the series diff --git a/mhkit/loads/general.py b/mhkit/loads/general.py index db670a855..1f70c7428 100644 --- a/mhkit/loads/general.py +++ b/mhkit/loads/general.py @@ -66,28 +66,19 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[],to_pandas=True): # loop through data_signal and get binned means for signal_name in data_signal: # Bin data - bin_stat = binned_statistic(bin_against,data[signal_name], + bin_stat_mean = binned_statistic(bin_against,data[signal_name], statistic='mean',bins=bin_edges) - # Calculate std of bins - std = [] - stdev = xr.DataArray(data = data[signal_name].values, - dims = 'binnumber', - coords = {'binnumber': bin_stat.binnumber}) - for i in range(1,len(bin_stat.bin_edges)): - try: - temp = stdev.sel(binnumber=i).std().values - std.append(np.array(temp,ndmin=1)[0]) - except: - std.append(np.nan) + bin_stat_std = binned_statistic(bin_against,data[signal_name], + statistic='std',bins=bin_edges) - bin_stat_list[signal_name] = ('index', bin_stat.statistic) - bin_std_list[signal_name] = ('index', bin_stat.statistic) + bin_stat_list[signal_name] = ('index', bin_stat_mean.statistic) + bin_std_list[signal_name] = ('index', bin_stat_std.statistic) # Convert to Datasets bin_mean = xr.Dataset(data_vars = bin_stat_list, - coords = {'index':np.arange(0,len(bin_stat.statistic))}) + coords = {'index':np.arange(0,len(bin_stat_mean.statistic))}) bin_std = xr.Dataset(data_vars = bin_std_list, - coords = {'index':np.arange(0,len(bin_stat.statistic))}) + coords = {'index':np.arange(0,len(bin_stat_std.statistic))}) # Check for nans for variable in list(bin_mean.variables): From da27418a4a7090ca7e0e2860c2cd6421b47b0ef7 Mon Sep 17 00:00:00 2001 From: akeeste Date: Fri, 1 Dec 2023 14:00:42 -0600 Subject: [PATCH 05/13] correct bin_statistics test data --- examples/data/loads/loads_data_dict.json | 324 +++++++++++------------ 1 file changed, 162 insertions(+), 162 deletions(-) diff --git a/examples/data/loads/loads_data_dict.json b/examples/data/loads/loads_data_dict.json index 3351ddbb5..9054afe7d 100644 --- a/examples/data/loads/loads_data_dict.json +++ b/examples/data/loads/loads_data_dict.json @@ -763,24 +763,24 @@ "yawoffset": 0.36065239549512096 }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 }, { "ActivePower": NaN, @@ -823,24 +823,24 @@ "yawoffset": NaN }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 }, { "ActivePower": NaN, @@ -863,24 +863,24 @@ "yawoffset": NaN }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 } ], "bin_means": [ @@ -1647,24 +1647,24 @@ "yawoffset": 0.32465542650598184 }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 }, { "ActivePower": NaN, @@ -1707,24 +1707,24 @@ "yawoffset": NaN }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 }, { "ActivePower": NaN, @@ -1747,24 +1747,24 @@ "yawoffset": NaN }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 } ], "bin_mins": [ @@ -2531,24 +2531,24 @@ "yawoffset": 11.605683455992253 }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 }, { "ActivePower": NaN, @@ -2591,24 +2591,24 @@ "yawoffset": NaN }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 }, { "ActivePower": NaN, @@ -2631,24 +2631,24 @@ "yawoffset": NaN }, { - "ActivePower": NaN, - "BL1_EdgeMom": NaN, - "BL1_FlapMom": NaN, - "BL3_EdgeMom": NaN, - "BL3_FlapMom": NaN, - "LSSDW_My": NaN, - "LSSDW_Mz": NaN, - "LSSDW_Tq": NaN, - "TB_ForeAft": NaN, - "TB_SideSide": NaN, - "TTTq": NaN, - "TT_ForeAft": NaN, - "TT_SideSide": NaN, - "WD_ModActive": NaN, - "WD_Nacelle": NaN, - "WD_NacelleMod": NaN, - "uWind_80m": NaN, - "yawoffset": NaN + "ActivePower": 0.0, + "BL1_EdgeMom": 0.0, + "BL1_FlapMom": 0.0, + "BL3_EdgeMom": 0.0, + "BL3_FlapMom": 0.0, + "LSSDW_My": 0.0, + "LSSDW_Mz": 0.0, + "LSSDW_Tq": 0.0, + "TB_ForeAft": 0.0, + "TB_SideSide": 0.0, + "TTTq": 0.0, + "TT_ForeAft": 0.0, + "TT_SideSide": 0.0, + "WD_ModActive": 0.0, + "WD_Nacelle": 0.0, + "WD_NacelleMod": 0.0, + "uWind_80m": 0.0, + "yawoffset": 0.0 } ], "loads": [ From 0d55ed8207fdd7c2ff5c354ec8ed30a0b392f80a Mon Sep 17 00:00:00 2001 From: akeeste Date: Fri, 1 Dec 2023 14:11:04 -0600 Subject: [PATCH 06/13] update bin_statistics test --- mhkit/tests/loads/test_loads.py | 36 ++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/mhkit/tests/loads/test_loads.py b/mhkit/tests/loads/test_loads.py index 61356ccb9..568491a2b 100644 --- a/mhkit/tests/loads/test_loads.py +++ b/mhkit/tests/loads/test_loads.py @@ -48,6 +48,8 @@ def test_bin_statistics(self): load_means = self.data['means'] bin_against = load_means['uWind_80m'] [b_means, b_means_std] = loads.general.bin_statistics(load_means, bin_against, bin_edges) + b_means.index.name = None # compatibility with old test data + b_means_std.index.name = None # compatibility with old test data assert_frame_equal(self.data['bin_means'],b_means) assert_frame_equal(self.data['bin_means_std'],b_means_std) @@ -61,6 +63,8 @@ def test_bin_statistics_xarray(self): load_means = load_means.to_xarray() bin_against = load_means['uWind_80m'] [b_means, b_means_std] = loads.general.bin_statistics(load_means, bin_against, bin_edges) + b_means.index.name = None # compatibility with old test data + b_means_std.index.name = None # compatibility with old test data assert_frame_equal(self.data['bin_means'],b_means) assert_frame_equal(self.data['bin_means_std'],b_means_std) @@ -96,13 +100,13 @@ def test_plot_statistics(self): # Generate plot loads.graphics.plot_statistics( self.data['means']['uWind_80m'], - self.data['means']['TB_ForeAft'], - self.data['maxs']['TB_ForeAft'], - self.data['mins']['TB_ForeAft'], - y_stdev=self.data['std']['TB_ForeAft'], - x_label='Wind Speed [m/s]', - y_label='Tower Base Mom [kNm]', - save_path=savepath) + self.data['means']['TB_ForeAft'], + self.data['maxs']['TB_ForeAft'], + self.data['mins']['TB_ForeAft'], + y_stdev=self.data['std']['TB_ForeAft'], + x_label='Wind Speed [m/s]', + y_label='Tower Base Mom [kNm]', + save_path=savepath) self.assertTrue(isfile(savepath)) @@ -158,6 +162,19 @@ def test_mler_coefficients(self): assert_series_equal(mler_data['Phase'], self.mler['phase'], check_exact=False, check_names=False, rtol=0.001) + def test_mler_coefficients_xarray(self): + Hs = 9.0 # significant wave height + Tp = 15.1 # time period of waves + pm = resource.pierson_moskowitz_spectrum(self.wave_freq, Tp, Hs) + mler_data = loads.extreme.mler_coefficients( + self.mler['RAO'].astype(complex).to_xarray(), pm, 1) + mler_data.reset_index(drop=True, inplace=True) + + assert_series_equal(mler_data['WaveSpectrum'], self.mler['Res_Spec'], + check_exact=False, check_names=False, atol=0.001) + assert_series_equal(mler_data['Phase'], self.mler['phase'], + check_exact=False, check_names=False, rtol=0.001) + def test_mler_simulation(self): T = np.linspace(-150, 150, 301) X = np.linspace(-300, 300, 601) @@ -189,6 +206,7 @@ def test_mler_export_time_series(self): RAO = self.mler['RAO'].astype(complex) mler_ts = loads.extreme.mler_export_time_series( RAO.values, mler, self.sim, k.k.values) + mler_ts.index.name = None # compatibility with old data assert_frame_equal(self.mler_ts, mler_ts, atol=0.0001) @@ -216,8 +234,8 @@ def test_longterm_extreme(self): def test_shortterm_extreme(self): methods = ['peaks_weibull', 'peaks_weibull_tail_fit', - 'peaks_over_threshold', 'block_maxima_gev', - 'block_maxima_gumbel'] + 'peaks_over_threshold', 'block_maxima_gev', + 'block_maxima_gumbel'] filename = "time_series_for_extremes.txt" data = np.loadtxt(os.path.join(datadir, filename)) t = data[:, 0] From 160807857d892af31f9205550323cfe1f3c4522d Mon Sep 17 00:00:00 2001 From: akeeste Date: Fri, 1 Dec 2023 14:11:37 -0600 Subject: [PATCH 07/13] fix dimension name in mler_wave_amp_normalize --- mhkit/loads/extreme.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index 2441fe3cb..e31624eff 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -856,7 +856,7 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, to_pandas=True): # rescale the wave spectral amplitude coefficients mler_norm = mler['WaveSpectrum'] * rescale_fact**2 mler_norm = mler_norm.to_dataset() - mler_norm = mler_norm.assign({'Phase': ('frequency', mler['Phase'].data)}) + mler_norm = mler_norm.assign({'Phase': (coord_names[0], mler['Phase'].data)}) if to_pandas: mler_norm = mler_norm.to_pandas() From ff7b3e7765276ac3e35ba9419586d5cef36e439f Mon Sep 17 00:00:00 2001 From: akeeste Date: Fri, 1 Dec 2023 15:13:20 -0600 Subject: [PATCH 08/13] formatting fixes --- mhkit/loads/extreme.py | 7 +------ mhkit/loads/general.py | 4 ++-- mhkit/tests/loads/test_loads.py | 20 ++++++++++---------- 3 files changed, 13 insertions(+), 18 deletions(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index e31624eff..d2f2549a5 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -711,11 +711,6 @@ def mler_coefficients(rao, wave_spectrum, response_desired, to_pandas=True): if response_desired < 0: _phase += np.pi - mler = pd.DataFrame( - data={'WaveSpectrum': _s, 'Phase': _phase}, index=freq_hz) - mler = mler.fillna(0) - - mler = xr.Dataset(data_vars = {'WaveSpectrum': (['frequency'], _s), 'Phase': (['frequency'], _phase)}, coords = {'frequency': freq_hz} @@ -738,7 +733,7 @@ def mler_simulation(parameters=None): Parameters ---------- - parameters: dict (fto_) + parameters: dict (optional) Simulation parameters. Keys: ----- diff --git a/mhkit/loads/general.py b/mhkit/loads/general.py index 1f70c7428..0244b2719 100644 --- a/mhkit/loads/general.py +++ b/mhkit/loads/general.py @@ -58,11 +58,11 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[],to_pandas=True): if not isinstance(data_signal, list): raise TypeError( f'data_signal must be of type list. Got: {type(data_signal)}') - + # Pre-allocate variable dictionaries bin_stat_list = {} bin_std_list = {} - + # loop through data_signal and get binned means for signal_name in data_signal: # Bin data diff --git a/mhkit/tests/loads/test_loads.py b/mhkit/tests/loads/test_loads.py index 568491a2b..6441aabba 100644 --- a/mhkit/tests/loads/test_loads.py +++ b/mhkit/tests/loads/test_loads.py @@ -99,14 +99,14 @@ def test_plot_statistics(self): savepath = abspath(join(testdir, 'test_scatplotter.png')) # Generate plot - loads.graphics.plot_statistics( self.data['means']['uWind_80m'], - self.data['means']['TB_ForeAft'], - self.data['maxs']['TB_ForeAft'], - self.data['mins']['TB_ForeAft'], - y_stdev=self.data['std']['TB_ForeAft'], - x_label='Wind Speed [m/s]', - y_label='Tower Base Mom [kNm]', - save_path=savepath) + loads.graphics.plot_statistics(self.data['means']['uWind_80m'], + self.data['means']['TB_ForeAft'], + self.data['maxs']['TB_ForeAft'], + self.data['mins']['TB_ForeAft'], + y_stdev=self.data['std']['TB_ForeAft'], + x_label='Wind Speed [m/s]', + y_label='Tower Base Mom [kNm]', + save_path=savepath) self.assertTrue(isfile(savepath)) @@ -234,8 +234,8 @@ def test_longterm_extreme(self): def test_shortterm_extreme(self): methods = ['peaks_weibull', 'peaks_weibull_tail_fit', - 'peaks_over_threshold', 'block_maxima_gev', - 'block_maxima_gumbel'] + 'peaks_over_threshold', 'block_maxima_gev', + 'block_maxima_gumbel'] filename = "time_series_for_extremes.txt" data = np.loadtxt(os.path.join(datadir, filename)) t = data[:, 0] From 136d7b565887305f4bef380edd5a561399e22bff Mon Sep 17 00:00:00 2001 From: akeeste Date: Mon, 11 Dec 2023 15:58:54 -0600 Subject: [PATCH 09/13] update return types, add optional dimension argument --- mhkit/loads/extreme.py | 121 +++++++++++++++++++++++------------------ mhkit/loads/general.py | 4 +- 2 files changed, 69 insertions(+), 56 deletions(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index d2f2549a5..c6fe255c1 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -611,7 +611,7 @@ def _cdf(self, x): return _LongTermExtreme(name="long_term_extreme", weights=weights, ste=ste) -def mler_coefficients(rao, wave_spectrum, response_desired, to_pandas=True): +def mler_coefficients(rao, wave_spectrum, response_desired, dimension="", to_pandas=True): """ Calculate MLER (most likely extreme response) coefficients from a sea state spectrum and a response RAO. @@ -620,17 +620,21 @@ def mler_coefficients(rao, wave_spectrum, response_desired, to_pandas=True): ---------- rao: numpy ndarray Response amplitude operator. - wave_spectrum: pd.DataFrame or xarray Dataset - Wave spectral density [m^2/Hz] indexed by frequency [Hz]. + wave_spectrum: pandas Series, pandas DataFrame, xarray DataArray, or xarray Dataset + Wave spectral density [m^2/Hz] indexed by frequency [Hz]. + DataFrame and Dataset inputs should only have one data variable response_desired: int or float Desired response, units should correspond to a motion RAO or units of force for a force RAO. + dimension: string (optional) + Name of the xarray dimension corresponding to time. + If not supplied, time is assumed to be the first dimension. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. Returns ------- - mler: pd.DataFrame + mler: pandas DataFrame or xarray Dataset DataFrame containing conditioned wave spectral amplitude coefficient [m^2-s], and Phase [rad] indexed by freq [Hz]. """ @@ -651,57 +655,58 @@ def mler_coefficients(rao, wave_spectrum, response_desired, to_pandas=True): if not isinstance(to_pandas, bool): raise TypeError( f'to_pandas must be of type bool. Got: {type(to_pandas)}') + + # Convert input to xarray DataArray + if isinstance(wave_spectrum, (pd.Series, pd.DataFrame)): + wave_spectrum = wave_spectrum.squeeze().to_xarray() - # If input is pandas, convert to xarray - if isinstance(wave_spectrum,pd.DataFrame): - wave_spectrum = wave_spectrum.to_xarray() - + if isinstance(wave_spectrum, xr.Dataset): + if len(wave_spectrum.data_vars) > 1: + raise ValueError( + f'wave_spectrum can only contain one variable. Got {list(wave_spectrum.data_vars)}.') + wave_spectrum = wave_spectrum.to_array() + + if dimension == "": + dimension = list(wave_spectrum.coords)[0] + # convert from Hz to rad/s - freq_hz = wave_spectrum.index.values + freq_hz = wave_spectrum.coords[dimension].values freq = freq_hz * (2*np.pi) - key_name = list(wave_spectrum.keys())[0] - wave_spectrum = wave_spectrum[key_name].values / (2*np.pi) - - # get delta - dw = (2*np.pi - 0.) / (len(freq)-1) - - spectrum_r = np.zeros(len(freq)) # [(response units)^2-s/rad] - _s = np.zeros(len(freq)) # [m^2-s/rad] - _a = np.zeros(len(freq)) # [m^2-s/rad] - _coeff_a_rn = np.zeros(len(freq)) # [1/(response units)] - _phase = np.zeros(len(freq)) - + wave_spectrum = wave_spectrum.to_numpy() / (2*np.pi) + + # get frequency step + dw = 2.0*np.pi / (len(freq)-1) + # Note: waves.A is "S" in Quon2016; 'waves' naming convention # matches WEC-Sim conventions (EWQ) # Response spectrum [(response units)^2-s/rad] -- Quon2016 Eqn. 3 - spectrum_r[:] = np.abs(rao)**2 * (2*wave_spectrum) - + spectrum_r = np.abs(rao)**2 * (2*wave_spectrum) + # calculate spectral moments and other important spectral values. m0 = (frequency_moment(pd.Series(spectrum_r, index=freq), 0)).iloc[0, 0] m1 = (frequency_moment(pd.Series(spectrum_r, index=freq), 1)).iloc[0, 0] m2 = (frequency_moment(pd.Series(spectrum_r, index=freq), 2)).iloc[0, 0] wBar = m1 / m0 - + # calculate coefficient A_{R,n} [(response units)^-1] -- Quon2016 Eqn. 8 # Drummen version. Dietz has negative of this. - _coeff_a_rn[:] = np.abs(rao) * np.sqrt(2*wave_spectrum*dw) * \ + _coeff_a_rn = np.abs(rao) * np.sqrt(2*wave_spectrum*dw) * \ ((m2 - freq*m1) + wBar*(freq*m0 - m1)) / (m0*m2 - m1**2) - + # save the new spectral info to pass out # Phase delay should be a positive number in this convention (AP) - _phase[:] = -np.unwrap(np.angle(rao)) - + _phase = -np.unwrap(np.angle(rao)) + # for negative values of Amp, shift phase by pi and flip sign # for negative amplitudes, add a pi phase shift, then flip sign on # negative Amplitudes _phase[_coeff_a_rn < 0] -= np.pi _coeff_a_rn[_coeff_a_rn < 0] *= -1 - + # calculate the conditioned spectrum [m^2-s/rad] - _s[:] = wave_spectrum * _coeff_a_rn[:]**2 * response_desired**2 - _a[:] = 2*wave_spectrum * _coeff_a_rn[:]**2 * \ - response_desired**2 - + _s = wave_spectrum * _coeff_a_rn**2 * response_desired**2 + _a = 2*wave_spectrum * _coeff_a_rn**2 * response_desired**2 + # if the response amplitude we ask for is negative, we will add # a pi phase shift to the phase information. This is because # the sign of self.desiredRespAmp is lost in the squaring above. @@ -710,16 +715,16 @@ def mler_coefficients(rao, wave_spectrum, response_desired, to_pandas=True): # new spectral information, S. (AP) if response_desired < 0: _phase += np.pi - + mler = xr.Dataset(data_vars = {'WaveSpectrum': (['frequency'], _s), 'Phase': (['frequency'], _phase)}, coords = {'frequency': freq_hz} ) mler.fillna(0) - + if to_pandas: mler = mler.to_pandas() - + return mler @@ -782,7 +787,7 @@ def mler_simulation(parameters=None): return sim -def mler_wave_amp_normalize(wave_amp, mler, sim, k, to_pandas=True): +def mler_wave_amp_normalize(wave_amp, mler, sim, k, dimension="", to_pandas=True): """ Function that renormalizes the incoming amplitude of the MLER wave to the desired peak height (peak to MSL). @@ -798,12 +803,15 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, to_pandas=True): 'mler_simulation'. k: numpy ndarray Wave number + dimension: string (optional) + Name of the xarray dimension corresponding to time. + If not supplied, time is assumed to be the first dimension. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. Returns ------- - mler_norm : pd.DataFrame + mler_norm : pandas DataFrame or xarray Dataset MLER coefficients """ try: @@ -829,11 +837,12 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, to_pandas=True): # If input is pandas, convert to xarray if isinstance(mler,pd.DataFrame): mler = mler.to_xarray() - - coord_names = list(mler.coords) - freq = mler.coords[coord_names[0]].values * 2*np.pi + + if dimension == "": + dimension = list(mler.coords)[0] + freq = mler.coords[dimension].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta - + wave_amp_time = np.zeros((sim['maxIX'], sim['maxIT'])) for ix, x in enumerate(sim['X']): for it, t in enumerate(sim['T']): @@ -842,24 +851,24 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, to_pandas=True): np.sqrt(2*mler['WaveSpectrum']*dw) * np.cos(freq*(t-sim['T0']) - k*(x-sim['X0']) + mler['Phase']) ) - + tmp_max_amp = np.max(np.abs(wave_amp_time)) - + # renormalization of wave amplitudes rescale_fact = np.abs(wave_amp) / np.abs(tmp_max_amp) - + # rescale the wave spectral amplitude coefficients mler_norm = mler['WaveSpectrum'] * rescale_fact**2 mler_norm = mler_norm.to_dataset() - mler_norm = mler_norm.assign({'Phase': (coord_names[0], mler['Phase'].data)}) - + mler_norm = mler_norm.assign({'Phase': (dimension, mler['Phase'].data)}) + if to_pandas: mler_norm = mler_norm.to_pandas() - + return mler_norm -def mler_export_time_series(rao, mler, sim, k, to_pandas=True): +def mler_export_time_series(rao, mler, sim, k, dimension="", to_pandas=True): """ Generate the wave amplitude time series at X0 from the calculated MLER coefficients @@ -875,12 +884,15 @@ def mler_export_time_series(rao, mler, sim, k, to_pandas=True): 'mler_simulation'. k: numpy ndarray Wave number. + dimension: string (optional) + Name of the xarray dimension corresponding to time. + If not supplied, time is assumed to be the first dimension. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. Returns ------- - mler_ts: pd.DataFrame + mler_ts: pandas DataFrame or xarray Dataset Time series of wave height [m] and linear response [*] indexed by time [s]. @@ -912,9 +924,10 @@ def mler_export_time_series(rao, mler, sim, k, to_pandas=True): # If input is pandas, convert to xarray if isinstance(mler,pd.DataFrame): mler = mler.to_xarray() - - coord_names = list(mler.coords) - freq = mler.coords[coord_names[0]].values * 2*np.pi + + if dimension == "": + dimension = list(mler.coords)[0] + freq = mler.coords[dimension].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta # calculate the series @@ -935,7 +948,7 @@ def mler_export_time_series(rao, mler, sim, k, to_pandas=True): mler_ts = xr.Dataset(data_vars = {'WaveHeight': (['time'], wave_amp_time[:,0]), 'LinearResponse': (['time'], wave_amp_time[:,1])}, coords = {'time': sim['T']}) - + if to_pandas: mler_ts = mler_ts.to_pandas() diff --git a/mhkit/loads/general.py b/mhkit/loads/general.py index 0244b2719..48b466b15 100644 --- a/mhkit/loads/general.py +++ b/mhkit/loads/general.py @@ -24,9 +24,9 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[],to_pandas=True): Returns -------- - bin_mean : pandas DataFrame + bin_mean : pandas DataFrame or xarray Dataset Mean of each bin - bin_std : pandas DataFrame + bin_std : pandas DataFrame or xarray Dataset Standard deviation of each bim """ From 2fac9db46a1700b118465234ffbdb4064f8d17a7 Mon Sep 17 00:00:00 2001 From: akeeste Date: Tue, 19 Dec 2023 14:10:51 -0600 Subject: [PATCH 10/13] add Series and DataArray to type check for mler_coefficients --- mhkit/loads/extreme.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index c6fe255c1..ceeba6dfd 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -646,9 +646,9 @@ def mler_coefficients(rao, wave_spectrum, response_desired, dimension="", to_pan if not isinstance(rao, np.ndarray): raise TypeError( f'rao must be of type np.ndarray. Got: {type(rao)}') - if not isinstance(wave_spectrum, (pd.DataFrame, xr.Dataset)): + if not isinstance(wave_spectrum, (pd.Series, pd.DataFrame, xr.DataArray, xr.Dataset)): raise TypeError( - f'wave_spectrum must be of type pd.DataFrame or xr.Dataset. Got: {type(wave_spectrum)}') + f'wave_spectrum must be of type pd.Series, pd.DataFrame, xr.DataArray, or xr.Dataset. Got: {type(wave_spectrum)}') if not isinstance(response_desired, (int, float)): raise TypeError( f'response_desired must be of type int or float. Got: {type(response_desired)}') From a3ee356f84393864ade7edb58f7ebf323f90e55c Mon Sep 17 00:00:00 2001 From: akeeste Date: Tue, 19 Dec 2023 14:25:01 -0600 Subject: [PATCH 11/13] update argument and docstring for time_dimension --- mhkit/loads/extreme.py | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index ceeba6dfd..3071212af 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -611,7 +611,7 @@ def _cdf(self, x): return _LongTermExtreme(name="long_term_extreme", weights=weights, ste=ste) -def mler_coefficients(rao, wave_spectrum, response_desired, dimension="", to_pandas=True): +def mler_coefficients(rao, wave_spectrum, response_desired, time_dimension="", to_pandas=True): """ Calculate MLER (most likely extreme response) coefficients from a sea state spectrum and a response RAO. @@ -626,9 +626,11 @@ def mler_coefficients(rao, wave_spectrum, response_desired, dimension="", to_pan response_desired: int or float Desired response, units should correspond to a motion RAO or units of force for a force RAO. - dimension: string (optional) - Name of the xarray dimension corresponding to time. - If not supplied, time is assumed to be the first dimension. + time_dimension: string (optional) + Name of the xarray dimension corresponding to time. Pandas only defines + 1 dimension, so for Pandas input the index is assumed to be time. + If not supplied for xarray input, time is assumed to be the first + dimension. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. @@ -787,7 +789,7 @@ def mler_simulation(parameters=None): return sim -def mler_wave_amp_normalize(wave_amp, mler, sim, k, dimension="", to_pandas=True): +def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas=True): """ Function that renormalizes the incoming amplitude of the MLER wave to the desired peak height (peak to MSL). @@ -803,9 +805,11 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, dimension="", to_pandas=True 'mler_simulation'. k: numpy ndarray Wave number - dimension: string (optional) - Name of the xarray dimension corresponding to time. - If not supplied, time is assumed to be the first dimension. + time_dimension: string (optional) + Name of the xarray dimension corresponding to time. Pandas only defines + 1 dimension, so for Pandas input the index is assumed to be time. + If not supplied for xarray input, time is assumed to be the first + dimension. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. @@ -868,7 +872,7 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, dimension="", to_pandas=True return mler_norm -def mler_export_time_series(rao, mler, sim, k, dimension="", to_pandas=True): +def mler_export_time_series(rao, mler, sim, k, time_dimension="", to_pandas=True): """ Generate the wave amplitude time series at X0 from the calculated MLER coefficients @@ -884,9 +888,11 @@ def mler_export_time_series(rao, mler, sim, k, dimension="", to_pandas=True): 'mler_simulation'. k: numpy ndarray Wave number. - dimension: string (optional) - Name of the xarray dimension corresponding to time. - If not supplied, time is assumed to be the first dimension. + time_dimension: string (optional) + Name of the xarray dimension corresponding to time. Pandas only defines + 1 dimension, so for Pandas input the index is assumed to be time. + If not supplied for xarray input, time is assumed to be the first + dimension. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. From a5182f167584bd141a5ae937a9462ecc164ffe95 Mon Sep 17 00:00:00 2001 From: akeeste Date: Tue, 19 Dec 2023 14:40:36 -0600 Subject: [PATCH 12/13] update dimension variables to time_dimension --- mhkit/loads/extreme.py | 38 ++++++++++++++++---------------------- 1 file changed, 16 insertions(+), 22 deletions(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index 3071212af..e0aaa393d 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -627,10 +627,8 @@ def mler_coefficients(rao, wave_spectrum, response_desired, time_dimension="", t Desired response, units should correspond to a motion RAO or units of force for a force RAO. time_dimension: string (optional) - Name of the xarray dimension corresponding to time. Pandas only defines - 1 dimension, so for Pandas input the index is assumed to be time. - If not supplied for xarray input, time is assumed to be the first - dimension. + Name of the xarray dimension corresponding to time. If not supplied, + defaults to the first dimension. Does not affect pandas input. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. @@ -668,11 +666,11 @@ def mler_coefficients(rao, wave_spectrum, response_desired, time_dimension="", t f'wave_spectrum can only contain one variable. Got {list(wave_spectrum.data_vars)}.') wave_spectrum = wave_spectrum.to_array() - if dimension == "": - dimension = list(wave_spectrum.coords)[0] + if time_dimension == "": + time_dimension = list(wave_spectrum.coords)[0] # convert from Hz to rad/s - freq_hz = wave_spectrum.coords[dimension].values + freq_hz = wave_spectrum.coords[time_dimension].values freq = freq_hz * (2*np.pi) wave_spectrum = wave_spectrum.to_numpy() / (2*np.pi) @@ -806,10 +804,8 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas k: numpy ndarray Wave number time_dimension: string (optional) - Name of the xarray dimension corresponding to time. Pandas only defines - 1 dimension, so for Pandas input the index is assumed to be time. - If not supplied for xarray input, time is assumed to be the first - dimension. + Name of the xarray dimension corresponding to time. If not supplied, + defaults to the first dimension. Does not affect pandas input. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. @@ -842,9 +838,9 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas if isinstance(mler,pd.DataFrame): mler = mler.to_xarray() - if dimension == "": - dimension = list(mler.coords)[0] - freq = mler.coords[dimension].values * 2*np.pi + if time_dimension == "": + time_dimension = list(mler.coords)[0] + freq = mler.coords[time_dimension].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta wave_amp_time = np.zeros((sim['maxIX'], sim['maxIT'])) @@ -864,7 +860,7 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas # rescale the wave spectral amplitude coefficients mler_norm = mler['WaveSpectrum'] * rescale_fact**2 mler_norm = mler_norm.to_dataset() - mler_norm = mler_norm.assign({'Phase': (dimension, mler['Phase'].data)}) + mler_norm = mler_norm.assign({'Phase': (time_dimension, mler['Phase'].data)}) if to_pandas: mler_norm = mler_norm.to_pandas() @@ -889,10 +885,8 @@ def mler_export_time_series(rao, mler, sim, k, time_dimension="", to_pandas=True k: numpy ndarray Wave number. time_dimension: string (optional) - Name of the xarray dimension corresponding to time. Pandas only defines - 1 dimension, so for Pandas input the index is assumed to be time. - If not supplied for xarray input, time is assumed to be the first - dimension. + Name of the xarray dimension corresponding to time. If not supplied, + defaults to the first dimension. Does not affect pandas input. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. @@ -931,9 +925,9 @@ def mler_export_time_series(rao, mler, sim, k, time_dimension="", to_pandas=True if isinstance(mler,pd.DataFrame): mler = mler.to_xarray() - if dimension == "": - dimension = list(mler.coords)[0] - freq = mler.coords[dimension].values * 2*np.pi + if time_dimension == "": + time_dimension = list(mler.coords)[0] + freq = mler.coords[time_dimension].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta # calculate the series From 6c10342e143e7596c1aee52eeed8b26dbdf0f729 Mon Sep 17 00:00:00 2001 From: akeeste Date: Tue, 19 Dec 2023 15:10:57 -0600 Subject: [PATCH 13/13] rename time_dimension to frequency_dimension where required --- mhkit/loads/extreme.py | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py index e0aaa393d..bbf4b813e 100644 --- a/mhkit/loads/extreme.py +++ b/mhkit/loads/extreme.py @@ -611,7 +611,7 @@ def _cdf(self, x): return _LongTermExtreme(name="long_term_extreme", weights=weights, ste=ste) -def mler_coefficients(rao, wave_spectrum, response_desired, time_dimension="", to_pandas=True): +def mler_coefficients(rao, wave_spectrum, response_desired, frequency_dimension="", to_pandas=True): """ Calculate MLER (most likely extreme response) coefficients from a sea state spectrum and a response RAO. @@ -626,8 +626,8 @@ def mler_coefficients(rao, wave_spectrum, response_desired, time_dimension="", t response_desired: int or float Desired response, units should correspond to a motion RAO or units of force for a force RAO. - time_dimension: string (optional) - Name of the xarray dimension corresponding to time. If not supplied, + frequency_dimension: string (optional) + Name of the xarray dimension corresponding to frequency. If not supplied, defaults to the first dimension. Does not affect pandas input. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. @@ -666,11 +666,11 @@ def mler_coefficients(rao, wave_spectrum, response_desired, time_dimension="", t f'wave_spectrum can only contain one variable. Got {list(wave_spectrum.data_vars)}.') wave_spectrum = wave_spectrum.to_array() - if time_dimension == "": - time_dimension = list(wave_spectrum.coords)[0] + if frequency_dimension == "": + frequency_dimension = list(wave_spectrum.coords)[0] # convert from Hz to rad/s - freq_hz = wave_spectrum.coords[time_dimension].values + freq_hz = wave_spectrum.coords[frequency_dimension].values freq = freq_hz * (2*np.pi) wave_spectrum = wave_spectrum.to_numpy() / (2*np.pi) @@ -787,7 +787,7 @@ def mler_simulation(parameters=None): return sim -def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas=True): +def mler_wave_amp_normalize(wave_amp, mler, sim, k, frequency_dimension="", to_pandas=True): """ Function that renormalizes the incoming amplitude of the MLER wave to the desired peak height (peak to MSL). @@ -803,8 +803,8 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas 'mler_simulation'. k: numpy ndarray Wave number - time_dimension: string (optional) - Name of the xarray dimension corresponding to time. If not supplied, + frequency_dimension: string (optional) + Name of the xarray dimension corresponding to frequency. If not supplied, defaults to the first dimension. Does not affect pandas input. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. @@ -838,9 +838,9 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas if isinstance(mler,pd.DataFrame): mler = mler.to_xarray() - if time_dimension == "": - time_dimension = list(mler.coords)[0] - freq = mler.coords[time_dimension].values * 2*np.pi + if frequency_dimension == "": + frequency_dimension = list(mler.coords)[0] + freq = mler.coords[frequency_dimension].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta wave_amp_time = np.zeros((sim['maxIX'], sim['maxIT'])) @@ -860,7 +860,7 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas # rescale the wave spectral amplitude coefficients mler_norm = mler['WaveSpectrum'] * rescale_fact**2 mler_norm = mler_norm.to_dataset() - mler_norm = mler_norm.assign({'Phase': (time_dimension, mler['Phase'].data)}) + mler_norm = mler_norm.assign({'Phase': (frequency_dimension, mler['Phase'].data)}) if to_pandas: mler_norm = mler_norm.to_pandas() @@ -868,7 +868,7 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k, time_dimension="", to_pandas return mler_norm -def mler_export_time_series(rao, mler, sim, k, time_dimension="", to_pandas=True): +def mler_export_time_series(rao, mler, sim, k, frequency_dimension="", to_pandas=True): """ Generate the wave amplitude time series at X0 from the calculated MLER coefficients @@ -884,8 +884,8 @@ def mler_export_time_series(rao, mler, sim, k, time_dimension="", to_pandas=True 'mler_simulation'. k: numpy ndarray Wave number. - time_dimension: string (optional) - Name of the xarray dimension corresponding to time. If not supplied, + frequency_dimension: string (optional) + Name of the xarray dimension corresponding to frequency. If not supplied, defaults to the first dimension. Does not affect pandas input. to_pandas: bool (optional) Flag to output pandas instead of xarray. Default = True. @@ -925,9 +925,9 @@ def mler_export_time_series(rao, mler, sim, k, time_dimension="", to_pandas=True if isinstance(mler,pd.DataFrame): mler = mler.to_xarray() - if time_dimension == "": - time_dimension = list(mler.coords)[0] - freq = mler.coords[time_dimension].values * 2*np.pi + if frequency_dimension == "": + frequency_dimension = list(mler.coords)[0] + freq = mler.coords[frequency_dimension].values * 2*np.pi dw = (max(freq) - min(freq)) / (len(freq)-1) # get delta # calculate the series