Speed up wave.resource module

mhkit/loads/extreme/mler.py

@@ -61,9 +61,9 @@ def _calculate_spectral_values(
  spectrum_r = np.abs(rao_array) ** 2 * (2 * wave_spectrum)
 
  # Calculate spectral moments
- m_0 = frequency_moment(pd.Series(spectrum_r, index=freq_hz), 0).iloc[0, 0]
- m_1 = frequency_moment(pd.Series(spectrum_r, index=freq_hz), 1).iloc[0, 0]
- m_2 = frequency_moment(pd.Series(spectrum_r, index=freq_hz), 2).iloc[0, 0]
+ m_0 = frequency_moment(pd.Series(spectrum_r, index=freq_hz), 0).item()
+ m_1 = frequency_moment(pd.Series(spectrum_r, index=freq_hz), 1).item()
+ m_2 = frequency_moment(pd.Series(spectrum_r, index=freq_hz), 2).item()
 
  # Calculate coefficient A_{R,n}
  coeff_a_rn = (

mhkit/tests/loads/test_loads.py

@@ -422,9 +422,9 @@ def test_mler_wave_amp_normalize(self):
  mler["WaveSpectrum"] = self.mler["Res_Spec"].values
  mler["Phase"] = self.mler["phase"].values
  k = resource.wave_number(wave_freq, 70)
- k = k.fillna(0)
+ np.nan_to_num(k, 0)
  mler_norm = loads.extreme.mler_wave_amp_normalize(
- 4.5 * 1.9, mler, self.sim, k.k.values
+ 4.5 * 1.9, mler, self.sim, k
  )
  mler_norm.reset_index(drop=True, inplace=True)
 
@@ -442,10 +442,10 @@ def test_mler_export_time_series(self):
  mler["WaveSpectrum"] = self.mler["Norm_Spec"].values
  mler["Phase"] = self.mler["phase"].values
  k = resource.wave_number(wave_freq, 70)
- k = k.fillna(0)
+ np.nan_to_num(k, 0)
  RAO = self.mler["RAO"].astype(complex)
  mler_ts = loads.extreme.mler_export_time_series(
- RAO.values, mler, self.sim, k.k.values
+ RAO.values, mler, self.sim, k
  )
  mler_ts.index.name = None # compatibility with old data
 

mhkit/tests/wave/test_performance.py

@@ -126,6 +126,7 @@ def test_powerperformance_workflow(self):
  if isfile(filename):
  os.remove(filename)
  P = pd.Series(np.random.normal(200, 40, 743), index=self.S.columns)
+ P.index.name = "variable"
  statistic = ["mean"]
  savepath = plotdir
  show_values = True

mhkit/tests/wave/test_resource_metrics.py

@@ -95,7 +95,8 @@ def test_kfromw(self):
 
  expected = self.valdata1[i]["k"]
  k = wave.resource.wave_number(f, h, rho)
- calculated = k.loc[:, "k"].values
+ # calculated = k.loc[:, "k"].values
+ calculated = k
  error = ((expected - calculated) ** 2).sum() # SSE
 
  self.assertLess(error, 1e-6)
@@ -106,7 +107,7 @@ def test_kfromw_one_freq(self):
  h = 1e9
  w = np.pi * 2 * f # deep water dispersion
  expected = w**2 / g
- calculated = wave.resource.wave_number(f=f, h=h, g=g).values[0][0]
+ calculated = wave.resource.wave_number(f=f, h=h, g=g).item()
  error = np.abs(expected - calculated)
  self.assertLess(error, 1e-6)
 
@@ -218,7 +219,7 @@ def test_moments(self):
 
  calculated = wave.resource.frequency_moment(
  S, int(m), frequency_bins=f_bins
- ).iloc[0, 0]
+ ).item()
  error = np.abs(expected - calculated) / expected
 
  self.assertLess(error, 0.01)
@@ -238,7 +239,7 @@ def test_energy_period_to_peak_period(self):
  f = np.linspace(1 / (10 * Tp), 3 / Tp, 100)
  S = wave.resource.jonswap_spectrum(f, Tp, Hs, g)
 
- Te_calc = wave.resource.energy_period(S).values[0][0]
+ Te_calc = wave.resource.energy_period(S).item()
 
  error = np.abs(T - Te_calc) / Te_calc
  self.assertLess(error, 0.01)
@@ -259,16 +260,16 @@ def test_metrics(self):
  expected = data["metrics"]["Hm0"]
  calculated = wave.resource.significant_wave_height(
  S, frequency_bins=f_bins
- ).iloc[0, 0]
+ )
  error = np.abs(expected - calculated) / expected
  # print('Hm0', expected, calculated, error)
  self.assertLess(error, 0.01)
 
  # Te
  expected = data["metrics"]["Te"]
- calculated = wave.resource.energy_period(S, frequency_bins=f_bins).iloc[
- 0, 0
- ]
+ calculated = wave.resource.energy_period(
+ S, frequency_bins=f_bins
+ ).item()
  error = np.abs(expected - calculated) / expected
  # print('Te', expected, calculated, error)
  self.assertLess(error, 0.01)
@@ -277,7 +278,7 @@ def test_metrics(self):
  expected = data["metrics"]["T0"]
  calculated = wave.resource.average_zero_crossing_period(
  S, frequency_bins=f_bins
- ).iloc[0, 0]
+ ).item()
  error = np.abs(expected - calculated) / expected
  # print('T0', expected, calculated, error)
  self.assertLess(error, 0.01)
@@ -289,7 +290,7 @@ def test_metrics(self):
  S,
  # Tc = Tavg**2
  frequency_bins=f_bins,
- ).iloc[0, 0]
+ ).item()
  ** 2
  )
  error = np.abs(expected - calculated) / expected
@@ -300,14 +301,14 @@ def test_metrics(self):
  expected = np.sqrt(data["metrics"]["Tm"])
  calculated = wave.resource.average_wave_period(
  S, frequency_bins=f_bins
- ).iloc[0, 0]
+ ).item()
  error = np.abs(expected - calculated) / expected
  # print('Tm', expected, calculated, error)
  self.assertLess(error, 0.01)
 
  # Tp
  expected = data["metrics"]["Tp"]
- calculated = wave.resource.peak_period(S).iloc[0, 0]
+ calculated = wave.resource.peak_period(S)
  error = np.abs(expected - calculated) / expected
  # print('Tp', expected, calculated, error)
  self.assertLess(error, 0.001)
@@ -316,7 +317,7 @@ def test_metrics(self):
  expected = data["metrics"]["e"]
  calculated = wave.resource.spectral_bandwidth(
  S, frequency_bins=f_bins
- ).iloc[0, 0]
+ ).item()
  error = np.abs(expected - calculated) / expected
  # print('e', expected, calculated, error)
  self.assertLess(error, 0.001)
@@ -325,8 +326,8 @@ def test_metrics(self):
  if file_i != "CDiP":
  for i, j in zip(data["h"], data["J"]):
  expected = data["J"][j]
- calculated = wave.resource.energy_flux(S, i)
- error = np.abs(expected - calculated.values) / expected
+ calculated = wave.resource.energy_flux(S, i).item()
+ error = np.abs(expected - calculated) / expected
  self.assertLess(error, 0.1)
 
  # v
@@ -335,14 +336,14 @@ def test_metrics(self):
  expected = data["metrics"]["v"]
  calculated = wave.resource.spectral_width(
  S, frequency_bins=f_bins
- ).iloc[0, 0]
+ ).item()
  error = np.abs(expected - calculated) / expected
  self.assertLess(error, 0.01)
 
  if file_i == "MC":
  expected = data["metrics"]["v"]
  # testing that default uniform frequency bin widths works
- calculated = wave.resource.spectral_width(S).iloc[0, 0]
+ calculated = wave.resource.spectral_width(S).item()
  error = np.abs(expected - calculated) / expected
  self.assertLess(error, 0.01)
 

mhkit/tests/wave/test_resource_spectrum.py

@@ -36,8 +36,8 @@ def tearDownClass(self):
 
  def test_pierson_moskowitz_spectrum(self):
  S = wave.resource.pierson_moskowitz_spectrum(self.f, self.Tp, self.Hs)
- Hm0 = wave.resource.significant_wave_height(S).iloc[0, 0]
- Tp0 = wave.resource.peak_period(S).iloc[0, 0]
+ Hm0 = wave.resource.significant_wave_height(S).item()
+ Tp0 = wave.resource.peak_period(S).item()
 
  errorHm0 = np.abs(self.Tp - Tp0) / self.Tp
  errorTp0 = np.abs(self.Hs - Hm0) / self.Hs
@@ -60,8 +60,8 @@ def test_pierson_moskowitz_spectrum_zero_freq(self):
 
  def test_jonswap_spectrum(self):
  S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
- Hm0 = wave.resource.significant_wave_height(S).iloc[0, 0]
- Tp0 = wave.resource.peak_period(S).iloc[0, 0]
+ Hm0 = wave.resource.significant_wave_height(S).item()
+ Tp0 = wave.resource.peak_period(S).item()
 
  errorHm0 = np.abs(self.Tp - Tp0) / self.Tp
  errorTp0 = np.abs(self.Hs - Hm0) / self.Hs
@@ -122,14 +122,14 @@ def test_surface_elevation_moments(self):
  eta, 1 / dt, len(eta.values), detrend=False, window="boxcar", noverlap=0
  )
 
- m0 = wave.resource.frequency_moment(S, 0).m0.values[0]
- m0n = wave.resource.frequency_moment(Sn, 0).m0.values[0]
+ m0 = wave.resource.frequency_moment(S, 0).item()
+ m0n = wave.resource.frequency_moment(Sn, 0).item()
  errorm0 = np.abs((m0 - m0n) / m0)
 
  self.assertLess(errorm0, 0.01)
 
- m1 = wave.resource.frequency_moment(S, 1).m1.values[0]
- m1n = wave.resource.frequency_moment(Sn, 1).m1.values[0]
+ m1 = wave.resource.frequency_moment(S, 1).item()
+ m1n = wave.resource.frequency_moment(Sn, 1).item()
  errorm1 = np.abs((m1 - m1n) / m1)
 
  self.assertLess(errorm1, 0.01)
@@ -168,7 +168,7 @@ def test_user_spectrum_without_frequency_index_name_defined(self):
 
  expected_magnitude = [-0.983917, 1.274248, -2.129812]
 
- assert_allclose(result["magnitude"], expected_magnitude, atol=1e-6)
+ assert_allclose(result, expected_magnitude, atol=1e-6)
 
  def test_ifft_sum_of_sines(self):
  S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)

mhkit/utils/type_handling.py

@@ -160,27 +160,25 @@ def convert_to_dataarray(data, name="data"):
  # iloc returns a Series with one value as expected.
  data = data.iloc[:, 0]
  else:
- index = data.index.values
- columns = data.columns.values
- data = xr.DataArray(
- data=data.T,
- dims=("variable", "index"),
- coords={"variable": columns, "index": index},
- )
+ # With this conversion, dataframe columns always become "dim_1".
+ # Rename to "variable" to match how multiple Dataset variables get converted into a DataArray dimension
+ data = xr.DataArray(data)
+ if data.dims[1] == "dim_1":
+ # Slight chance their is already a name for the columns
+ data = data.rename({"dim_1": "variable"})
 
  # Checks xr.Dataset input and converts to xr.DataArray if possible
  if isinstance(data, xr.Dataset):
  keys = list(data.keys())
  if len(keys) == 1:
- # if only one variable, remove the "variable" dimension and rename the DataArray to simplify
- data = data.to_array()
- data = data.sel(variable=keys[0])
- data.name = keys[0]
- data.drop_vars("variable")
+ # if only one variable, select that variable so reduce the Dataset to a DataArray
+ data = data[keys[0]]
  else:
  # Allow multiple variables if they have the same dimensions
  if all([data[keys[0]].dims == data[key].dims for key in keys]):
- data = data.to_array()
+ data = (
+ data.to_array().T
+ ) # transpose so that the new "variable dimension" is the last dimension (matches DataFrame to DataArray behavior)
  else:
  raise ValueError(
  "Multivariate Datasets can only be input if all variables have the same dimensions."

mhkit/wave/performance.py

@@ -6,7 +6,7 @@
 from mhkit import wave
 import matplotlib.pylab as plt
 from os.path import join
-from mhkit.utils import convert_to_dataarray, convert_to_dataset
+from mhkit.utils import convert_to_dataarray
 
 
 def capture_length(P, J, to_pandas=True):
@@ -246,22 +246,22 @@ def power_matrix(LM, JM):
 
  Parameters
  ------------
- LM: pandas DataFrame or xarray Dataset
+ LM: pandas DataFrame, xarray DatArray, or xarray Dataset
  Capture length matrix
- JM: pandas DataFrame or xarray Dataset
+ JM: pandas DataFrame, xarray DatArray, or xarray Dataset
  Wave energy flux matrix
 
  Returns
  ---------
- PM: pandas DataFrame or xarray Dataset
+ PM: pandas DataFrame, xarray DatArray, or xarray Dataset
  Power matrix
 
  """
- if not isinstance(LM, (pd.DataFrame, xr.Dataset)):
+ if not isinstance(LM, (pd.DataFrame, xr.DataArray, xr.Dataset)):
  raise TypeError(
  f"LM must be of type pd.DataFrame or xr.Dataset. Got: {type(LM)}"
  )
- if not isinstance(JM, (pd.DataFrame, xr.Dataset)):
+ if not isinstance(JM, (pd.DataFrame, xr.DataArray, xr.Dataset)):
  raise TypeError(
  f"JM must be of type pd.DataFrame or xr.Dataset. Got: {type(JM)}"
  )
@@ -306,11 +306,11 @@ def mean_annual_energy_production_matrix(LM, JM, frequency):
 
  Parameters
  ------------
- LM: pandas DataFrame or xarray Dataset
+ LM: pandas DataFrame, xarray DatArray, or xarray Dataset
  Capture length
- JM: pandas DataFrame or xarray Dataset
+ JM: pandas DataFrame, xarray DatArray, or xarray Dataset
  Wave energy flux
- frequency: pandas DataFrame or xarray Dataset
+ frequency: pandas DataFrame, xarray DatArray, or xarray Dataset
  Data frequency for each bin
 
  Returns
@@ -393,7 +393,7 @@ def power_performance_workflow(
  maep_matrix: float
  Mean annual energy production
  """
- S = convert_to_dataset(S)
+ S = convert_to_dataarray(S)
  if not isinstance(h, (int, float)):
  raise TypeError(f"h must be of type int or float. Got: {type(h)}")
  P = convert_to_dataarray(P)
@@ -408,19 +408,16 @@ def power_performance_workflow(
 
  # Compute the enegy periods from the spectra data
  Te = wave.resource.energy_period(S, frequency_bins=frequency_bins, to_pandas=False)
- Te = Te["Te"]
 
  # Compute the significant wave height from the NDBC spectra data
  Hm0 = wave.resource.significant_wave_height(
  S, frequency_bins=frequency_bins, to_pandas=False
  )
- Hm0 = Hm0["Hm0"]
 
  # Compute the energy flux from spectra data and water depth
  J = wave.resource.energy_flux(
  S, h, deep=deep, rho=rho, g=g, ratio=ratio, to_pandas=False
  )
- J = J["J"]
 
  # Calculate capture length from power and energy flux
  L = wave.performance.capture_length(P, J, to_pandas=False)