Convert wave module to xarray (#302)

* fix type error message in power.quality

* convert wave.performance to xarray

* remove extraneous variables from test_performance

* wave/resource 90 percent converted to xarray

* speed up conversion to dataset

* wave.resource conversion complete

* allow xarray input to wave.contours functions

* allow xarray in wave.graphics

* minor cleanup

* update type handling in wave_number and depth_regime

* update error messages in convert_to_dataarray

* test_wave_metrics passing

* correct typo in graphics

* wave tests passing except frequency binned surface_elevation calls

* fix bug/memory error with surface_elevation frequency bins

* black formatting

* allow xarray output in wave.io.cdip

* allow xarray output from wave.io.wecsim

* make utils function to convert nesed dicts of dataframes to xarray

* allow xarray input/output in wave.io.swan

* allow xarray input/output in wave.io.ndbc

* add xarray test for test_wecsim.py

* add xarray test in test_cdip

* update dataframe to dataarray case in utils.type_handling

* add xarray test to test_ndbc

* add xarray tests to test_swan

* black formatting

* fix single value dataframe to dataarray conversion

* black formatting

* update docstring in type_handling

* update description of to_pandas flag when returning a dict of pd or xr

* update ndbc docstrings

* move TODO from comment to an issue

* update error message

* update output types for wavelength

* add frequency_dimension and time_dimension kwargs

* add check for shape of phases and S in surface_elevation

* black formatting

* pylint suggestion in power.characteristics

* fix cdip example

* allow xarray output in wave.io.hindcast functions

mhkit/power/characteristics.py

@@ -229,11 +229,11 @@ def ac_power_three_phase(
     current = convert_to_dataset(current, "current")
 
     # Check that sizes are the same
-    if not len(voltage.data_vars) == 3:
+    if len(voltage.data_vars) != 3:
         raise ValueError("voltage must have three columns")
-    if not len(current.data_vars) == 3:
+    if len(current.data_vars) != 3:
         raise ValueError("current must have three columns")
-    if not current.sizes == voltage.sizes:
+    if current.sizes != voltage.sizes:
         raise ValueError("current and voltage must be of the same size")
 
     power = dc_power(voltage, current, to_pandas=False)["Gross"]

mhkit/power/quality.py

@@ -167,7 +167,7 @@ def harmonic_subgroups(
 
     if not isinstance(frequency_dimension, str):
         raise TypeError(
-            f"frequency_dimension must be of type bool. Got: {type(frequency_dimension)}"
+            f"frequency_dimension must be of type str. Got: {type(frequency_dimension)}"
         )
 
     # Convert input to xr.Dataset

mhkit/tests/wave/io/test_cdip.py

@@ -84,7 +84,9 @@ def test_dates_to_timestamp(self):
         self.assertEqual(end_dt, end_date)
 
     def test_get_netcdf_variables_all2Dvars(self):
-        data = wave.io.cdip.get_netcdf_variables(self.test_nc, all_2D_variables=True)
+        data = wave.io.cdip.get_netcdf_variables(
+            self.test_nc, all_2D_variables=True, to_pandas=False
+        )
         returned_keys = [key for key in data["data"]["wave2D"].keys()]
         self.assertTrue(set(returned_keys) == set(self.vars2D))
 

mhkit/tests/wave/io/test_ndbc.py

@@ -5,6 +5,7 @@
 import mhkit.wave as wave
 from io import StringIO
 import pandas as pd
+import xarray as xr
 import numpy as np
 import contextlib
 import unittest
@@ -117,8 +118,9 @@ def test_ndbnc_read_historical_met(self):
 
     # Spectral data
     def test_ndbc_read_spectral(self):
-        data, units = wave.io.ndbc.read_file(join(datadir, "data.txt"))
-        self.assertEqual(data.shape, (743, 47))
+        data, units = wave.io.ndbc.read_file(join(datadir, "data.txt"), to_pandas=False)
+        self.assertEqual(len(data.data_vars), 47)
+        self.assertEqual(len(data["dim_0"]), 743)
         self.assertEqual(units, None)
 
     # Continuous wind data
@@ -157,8 +159,8 @@ def test__ndbc_parse_filenames(self):
 
     def test_ndbc_request_data(self):
         filenames = pd.Series(self.filenames[0])
-        ndbc_data = wave.io.ndbc.request_data("swden", filenames)
-        self.assertTrue(self.swden.equals(ndbc_data["1996"]))
+        ndbc_data = wave.io.ndbc.request_data("swden", filenames, to_pandas=False)
+        self.assertTrue(xr.Dataset(self.swden).equals(ndbc_data["1996"]))
 
     def test_ndbc_request_data_from_dataframe(self):
         filenames = pd.DataFrame(pd.Series(data=self.filenames[0]))

mhkit/tests/wave/io/test_swan.py

@@ -9,6 +9,7 @@
 import mhkit.wave as wave
 from io import StringIO
 import pandas as pd
+import xarray as xr
 import numpy as np
 import contextlib
 import unittest
@@ -62,6 +63,14 @@ def test_read_block_txt(self):
         sumSum = swanBlockTxt["Significant wave height"].sum().sum()
         self.assertAlmostEqual(self.expected_table["Hsig"].sum(), sumSum, places=-2)
 
+    def test_read_block_txt_xarray(self):
+        swanBlockTxt, metaData = wave.io.swan.read_block(
+            self.swan_block_txt_file, to_pandas=False
+        )
+        self.assertEqual(len(swanBlockTxt), 4)
+        sumSum = swanBlockTxt["Significant wave height"].sum().sum()
+        self.assertAlmostEqual(self.expected_table["Hsig"].sum(), sumSum, places=-2)
+
     def test_block_to_table(self):
         x = np.arange(5)
         y = np.arange(5, 10)

mhkit/tests/wave/io/test_wecsim.py

@@ -51,12 +51,13 @@ def test_read_wecSim_no_mooring(self):
     ### WEC-Sim data, with cable
     def test_read_wecSim_cable(self):
         ws_output = wave.io.wecsim.read_output(
-            join(datadir, "Cable_matlabWorkspace_structure.mat")
+            join(datadir, "Cable_matlabWorkspace_structure.mat"),
+            to_pandas=False,
+        )
+        self.assertEqual(ws_output["wave"]["elevation"].name, "elevation")
+        self.assertEqual(
+            ws_output["bodies"]["body1"]["position_dof1"].name, "position_dof1"
         )
-        self.assertEqual(ws_output["wave"].elevation.name, "elevation")
-        self.assertEqual(ws_output["bodies"]["body1"].name, "BuoyDraft5cm")
-        self.assertEqual(ws_output["cables"].name, "Cable")
-        self.assertEqual(ws_output["constraints"]["constraint1"].name, "Mooring")
         self.assertEqual(len(ws_output["mooring"]), 0)
         self.assertEqual(len(ws_output["moorDyn"]), 0)
         self.assertEqual(len(ws_output["ptosim"]), 0)

mhkit/tests/wave/test_contours.py

@@ -68,7 +68,7 @@ def test_environmental_contour(self):
     def test_environmental_contours_invalid_inputs(self):
         # Invalid x1 tests
         x1_non_numeric = "not an array"
-        with self.assertRaises(ValueError):
+        with self.assertRaises(TypeError):
             wave.contours.environmental_contours(
                 x1_non_numeric, self.wdrt_Te, 3600, 50, "PCA"
             )
@@ -81,7 +81,7 @@ def test_environmental_contours_invalid_inputs(self):
 
         # Invalid x2 tests
         x2_non_numeric = "not an array"
-        with self.assertRaises(ValueError):
+        with self.assertRaises(TypeError):
             wave.contours.environmental_contours(
                 self.wdrt_Hm0, x2_non_numeric, 3600, 50, "PCA"
             )
@@ -264,7 +264,7 @@ def test_PCA_contours_invalid_inputs(self):
 
         # Invalid x1 tests
         x1_non_numeric = "not an array"
-        with self.assertRaises(ValueError):
+        with self.assertRaises(TypeError):
             wave.contours.PCA_contour(
                 x1_non_numeric, self.wdrt_Te, copula["PCA_fit"], PCA_args
             )
@@ -277,7 +277,7 @@ def test_PCA_contours_invalid_inputs(self):
 
         # Invalid x2 tests
         x2_non_numeric = "not an array"
-        with self.assertRaises(ValueError):
+        with self.assertRaises(TypeError):
             wave.contours.PCA_contour(
                 self.wdrt_Hm0, x2_non_numeric, copula["PCA_fit"], PCA_args
             )

mhkit/tests/wave/test_performance.py

@@ -1,23 +1,10 @@
 from os.path import abspath, dirname, join, isfile, normpath, relpath
-from pandas.testing import assert_frame_equal
-from numpy.testing import assert_allclose
-from scipy.interpolate import interp1d
-from random import seed, randint
 import matplotlib.pylab as plt
-from datetime import datetime
 import xarray.testing as xrt
 import mhkit.wave as wave
-from io import StringIO
 import pandas as pd
 import numpy as np
-import contextlib
 import unittest
-import netCDF4
-import inspect
-import pickle
-import time
-import json
-import sys
 import os
 
 
@@ -144,7 +131,6 @@ def test_powerperformance_workflow(self):
         show_values = True
         h = 60
         expected = 401239.4822345051
-        x = self.S.T
         CM, MAEP = wave.performance.power_performance_workflow(
             self.S, h, P, statistic, savepath=savepath, show_values=show_values
         )

mhkit/tests/wave/test_resource_metrics.py

@@ -9,6 +9,7 @@
 import mhkit.wave as wave
 from io import StringIO
 import pandas as pd
+import xarray as xr
 import numpy as np
 import contextlib
 import unittest
@@ -110,37 +111,44 @@ def test_kfromw_one_freq(self):
         self.assertLess(error, 1e-6)
 
     def test_wave_length(self):
-        k_list = [1, 2, 10, 3]
-        l_expected = (2.0 * np.pi / np.array(k_list)).tolist()
+        k_array = np.asarray([1.0, 2.0, 10.0, 3.0])
 
-        k_df = pd.DataFrame(k_list, index=[1, 2, 3, 4])
+        k_int = int(k_array[0])
+        k_float = k_array[0]
+        k_df = pd.DataFrame(k_array, index=[1, 2, 3, 4])
         k_series = k_df[0]
-        k_array = np.array(k_list)
 
-        for l in [k_list, k_df, k_series, k_array]:
+        for l in [k_array, k_int, k_float, k_df, k_series]:
             l_calculated = wave.resource.wave_length(l)
-            self.assertListEqual(l_expected, l_calculated.tolist())
-
-        idx = 0
-        k_int = k_list[idx]
-        l_calculated = wave.resource.wave_length(k_int)
-        self.assertEqual(l_expected[idx], l_calculated)
+            self.assertTrue(np.all(2.0 * np.pi / l == l_calculated))
 
     def test_depth_regime(self):
-        expected = [True, True, False, True]
-        l_list = [1, 2, 10, 3]
-        l_df = pd.DataFrame(l_list, index=[1, 2, 3, 4])
-        l_series = l_df[0]
-        l_array = np.array(l_list)
         h = 10
-        for l in [l_list, l_df, l_series, l_array]:
+
+        # non-array like formats
+        l_int = 1
+        l_float = 1.0
+        expected = True
+        for l in [l_int, l_float]:
+            calculated = wave.resource.depth_regime(l, h)
+            self.assertTrue(np.all(expected == calculated))
+
+        # array-like formats
+        l_array = np.array([1, 2, 10, 3])
+        l_df = pd.DataFrame(l_array, index=[1, 2, 3, 4])
+        l_series = l_df[0]
+        l_da = xr.DataArray(l_series)
+        l_da.name = "data"
+        l_ds = l_da.to_dataset()
+        expected = [True, True, False, True]
+        for l in [l_array, l_series, l_da, l_ds]:
             calculated = wave.resource.depth_regime(l, h)
-            self.assertListEqual(expected, calculated.tolist())
+            self.assertTrue(np.all(expected == calculated))
 
-        idx = 0
-        l_int = l_list[idx]
-        calculated = wave.resource.depth_regime(l_int, h)
-        self.assertEqual(expected[idx], calculated)
+        # special formatting for pd.DataFrame
+        for l in [l_df]:
+            calculated = wave.resource.depth_regime(l, h)
+            self.assertTrue(np.all(expected == calculated[0]))
 
     def test_wave_celerity(self):
         # Depth regime ratio
@@ -166,10 +174,10 @@ def test_wave_celerity(self):
         self.assertTrue(all(np.pi * f / k.squeeze().values == cg.squeeze().values))
 
     def test_energy_flux_deep(self):
-        # Dependent on mhkit.resource.BS spectrum
         S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
         Te = wave.resource.energy_period(S)
         Hm0 = wave.resource.significant_wave_height(S)
+
         rho = 1025
         g = 9.80665
         coeff = rho * (g**2) / (64 * np.pi)
@@ -180,6 +188,21 @@ def test_energy_flux_deep(self):
 
         self.assertTrue(J_calc.squeeze() == J)
 
+    def test_energy_flux_shallow(self):
+        S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
+        Te = wave.resource.energy_period(S)
+        Hm0 = wave.resource.significant_wave_height(S)
+
+        rho = 1025
+        g = 9.80665
+        coeff = rho * (g**2) / (64 * np.pi)
+        J = coeff * (Hm0.squeeze() ** 2) * Te.squeeze()
+
+        h = 1000  # effectively deep but without assumptions
+        J_calc = wave.resource.energy_flux(S, h, deep=False)
+        err = np.abs(J_calc.squeeze() - J)
+        self.assertLess(err, 1e-6)
+
     def test_moments(self):
         for file_i in self.valdata2.keys():  # for each file MC, AH, CDiP
             datasets = self.valdata2[file_i]

mhkit/tests/wave/test_resource_spectrum.py

@@ -2,22 +2,12 @@
 from pandas.testing import assert_frame_equal
 from numpy.testing import assert_allclose
 from scipy.interpolate import interp1d
-from random import seed, randint
 import matplotlib.pylab as plt
-from datetime import datetime
-import xarray.testing as xrt
+import xarray as xr
 import mhkit.wave as wave
-from io import StringIO
 import pandas as pd
 import numpy as np
-import contextlib
 import unittest
-import netCDF4
-import inspect
-import pickle
-import time
-import json
-import sys
 import os
 
 
@@ -89,18 +79,19 @@ def test_jonswap_spectrum_zero_freq(self):
         self.assertEqual(S_zero.values.squeeze()[0], 0.0)
         self.assertGreater(S_nonzero.values.squeeze()[0], 0.0)
 
-    def test_surface_elevation_phases_np_and_pd(self):
+    def test_surface_elevation_phases_xr_and_pd(self):
         S0 = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
         S1 = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs * 1.1)
         S = pd.concat([S0, S1], axis=1)
 
         phases_np = np.random.rand(S.shape[0], S.shape[1]) * 2 * np.pi
         phases_pd = pd.DataFrame(phases_np, index=S.index, columns=S.columns)
+        phases_xr = xr.Dataset(phases_pd)
 
-        eta_np = wave.resource.surface_elevation(S, self.t, phases=phases_np, seed=1)
+        eta_xr = wave.resource.surface_elevation(S, self.t, phases=phases_xr, seed=1)
         eta_pd = wave.resource.surface_elevation(S, self.t, phases=phases_pd, seed=1)
 
-        assert_frame_equal(eta_np, eta_pd)
+        assert_frame_equal(eta_xr, eta_pd)
 
     def test_surface_elevation_frequency_bins_np_and_pd(self):
         S0 = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
@@ -151,7 +142,8 @@ def test_surface_elevation_rmse(self):
         )
 
         fSn = interp1d(Sn.index.values, Sn.values, axis=0)
-        rmse = (S.values - fSn(S.index.values)) ** 2
+        Sn_interp = fSn(S.index.values).squeeze()
+        rmse = (S.values.squeeze() - Sn_interp) ** 2
         rmse_sum = (np.sum(rmse) / len(rmse)) ** 0.5
 
         self.assertLess(rmse_sum, 0.02)

mhkit/utils/__init__.py

@@ -8,6 +8,11 @@
 )
 from .cache import handle_caching, clear_cache
 from .upcrossing import upcrossing, peaks, troughs, heights, periods, custom
-from .type_handling import to_numeric_array, convert_to_dataset, convert_to_dataarray
+from .type_handling import (
+    to_numeric_array,
+    convert_to_dataset,
+    convert_to_dataarray,
+    convert_nested_dict_and_pandas,
+)
 
 _matlab = False  # Private variable indicating if mhkit is run through matlab