Repository: MHKiT-Software/MHKiT-Python (language: python)

Issue:
Assert Statements are for Debugging
Currently MHKiT uses assert staements in its production code (primarilty for type checking functional inputs). As stated in the docs (https://docs.python.org/3/reference/simple_stmts.html#the-assert-statement) assert stements are for debugging. 

We should raise an error if the user does not pass the correct information. 

So going forward we would convert:
```python
assert isinstance(station_number, (str, type(None))), (f'station_number must be  of type string')
```
to this

```python
if not isinstance(station_number, (str, type(None))):
    raise ValueError(f'station_number must be of type string. Got: {station_number}')
```
D3D
Added limits to `variable_interpolation`  and added 3 array input capability to `create_points`


Discussion / hints:
Here is an article covering this idea in more detail:
https://snyk.io/blog/the-dangers-of-assert-in-python/

![image](https://github.com/MHKiT-Software/MHKiT-Python/assets/13438942/8affe511-ab42-4a41-8316-65e1a5dc36f9)

Hello! I'm trying to get my feet wet with some python open source. This seems straight forward so I'm happy to take care of it. Does this seem like a good first issue? Thanks!
> Hello! I'm trying to get my feet wet with some python open source. This seems straight forward so I'm happy to take care of it. Does this seem like a good first issue? Thanks!

@onedeeper  I am so sorry I missed this. Yes this would be a great first issue. Please lmk if you need any help contributing!


Fix patch:
diff --git a/examples/directional_waves.ipynb b/examples/directional_waves.ipynb
index 73ed2c6ff..4d05a3822 100644
--- a/examples/directional_waves.ipynb
+++ b/examples/directional_waves.ipynb
@@ -1165,7 +1165,7 @@
     }
    ],
    "source": [
-    "wave.graphics.plot_directional_spectrum(spectrum, min=0.3)"
+    "wave.graphics.plot_directional_spectrum(spectrum, color_level_min=0.3)"
    ]
   },
   {
@@ -1195,7 +1195,7 @@
     }
    ],
    "source": [
-    "wave.graphics.plot_directional_spectrum(spectrum, min=0.3, fill=False, nlevels=4)"
+    "wave.graphics.plot_directional_spectrum(spectrum, color_level_min=0.3, fill=False, nlevels=4)"
    ]
   },
   {
@@ -1318,7 +1318,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.9.13 ('.venv': venv)",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -1332,9 +1332,8 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.9.17"
   },
-  "orig_nbformat": 4,
   "vscode": {
    "interpreter": {
     "hash": "15fd306e44580d8cf431083454c399b84f9cc4f7f2c761501397671836835f49"
diff --git a/mhkit/dolfyn/velocity.py b/mhkit/dolfyn/velocity.py
index 47d3a6528..4ab0194bc 100644
--- a/mhkit/dolfyn/velocity.py
+++ b/mhkit/dolfyn/velocity.py
@@ -861,8 +861,9 @@ def power_spectral_density(self, veldat,
 
         # Spectra, if input is full velocity or a single array
         if len(vel.shape) == 2:
-            assert vel.shape[0] == 3, "Function can only handle 1D or 3D arrays." \
-                " If ADCP data, please select a specific depth bin."
+            if not vel.shape[0] == 3:
+                raise Exception("Function can only handle 1D or 3D arrays." \
+                " If ADCP data, please select a specific depth bin.")
             if (noise is not None) and (np.shape(noise)[0] != 3):
                 raise Exception(
                     'Noise should have same first dimension as velocity')
diff --git a/mhkit/loads/extreme.py b/mhkit/loads/extreme.py
index c364b14a5..4cff02ca0 100644
--- a/mhkit/loads/extreme.py
+++ b/mhkit/loads/extreme.py
@@ -59,8 +59,10 @@ def global_peaks(t, data):
     peaks: np.array
         Peak values of the response time-series
     """
-    assert isinstance(t, np.ndarray), 't must be of type np.ndarray'
-    assert isinstance(data, np.ndarray), 'data must be of type np.ndarray'
+    if not isinstance(t, np.ndarray):
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
 
     # Find zero up-crossings
     inds = upcrossing(t, data)
@@ -99,9 +101,12 @@ def number_of_short_term_peaks(n, t, t_st):
     n_st : float
         Number of peaks in short term period.
     """
-    assert isinstance(n, int), 'n must be of type int'
-    assert isinstance(t, float), 't must be of type float'
-    assert isinstance(t_st, float), 't_st must be of type float'
+    if not isinstance(n, int):
+        raise TypeError(f'n must be of type int. Got: {type(n)}')
+    if not isinstance(t, float):
+        raise TypeError(f't must be of type float. Got: {type(t)}')
+    if not isinstance(t_st, float):
+        raise TypeError(f't_st must be of type float. Got: {type(t_st)}')
 
     return n * t_st / t
 
@@ -124,7 +129,8 @@ def peaks_distribution_weibull(x):
     peaks: scipy.stats.rv_frozen
         Probability distribution of the peaks.
     """
-    assert isinstance(x, np.ndarray), 'x must be of type np.ndarray'
+    if not isinstance(x, np.ndarray):
+        raise TypeError(f'x must be of type np.ndarray. Got: {type(x)}')
 
     # peaks distribution
     peaks_params = stats.exponweib.fit(x, f0=1, floc=0)
@@ -154,7 +160,8 @@ def peaks_distribution_weibull_tail_fit(x):
     peaks: scipy.stats.rv_frozen
         Probability distribution of the peaks.
     """
-    assert isinstance(x, np.ndarray), 'x must be of type np.ndarray'
+    if not isinstance(x, np.ndarray):
+        raise TypeError(f'x must be of type np.ndarray. Got: {type(x)}')
 
     # Initial guess for Weibull parameters
     p0 = stats.exponweib.fit(x, f0=1, floc=0)
@@ -231,16 +238,17 @@ def automatic_hs_threshold(
         Threshold that results in the best correlation.
     """
     if not isinstance(sampling_rate, (float, int)):
-        raise TypeError('sampling_rate must be of type float or int')
-    
+        raise TypeError(
+            f'sampling_rate must be of type float or int. Got: {type(sampling_rate)}')
     if not isinstance(peaks, np.ndarray):
-        raise TypeError('peaks must be of type np.ndarray')
-    
+        raise TypeError(
+            f'peaks must be of type np.ndarray. Got: {type(peaks)}')
     if not len(initial_threshold_range) == 3:
-        raise ValueError('initial_threshold_range must be length 3')
-    
+        raise ValueError(
+            f'initial_threshold_range must be length 3. Got: {len(initial_threshold_range)}')
     if not isinstance(max_refinement, int):
-        raise TypeError('max_refinement must be of type int')
+        raise TypeError(
+            f'max_refinement must be of type int. Got: {type(max_refinement)}')
 
     range_min, range_max, range_step = initial_threshold_range
     best_threshold = -1
@@ -309,11 +317,13 @@ def peaks_distribution_peaks_over_threshold(x, threshold=None):
     peaks: scipy.stats.rv_frozen
         Probability distribution of the peaks.
     """
-    assert isinstance(x, np.ndarray), 'x must be of type np.ndarray'
+    if not isinstance(x, np.ndarray):
+        raise TypeError(f'x must be of type np.ndarray. Got: {type(x)}')
     if threshold is None:
         threshold = np.mean(x) + 1.4 * np.std(x)
-    assert isinstance(threshold, float
-                      ), 'threshold must be of type float'
+    if not isinstance(threshold, float):
+        raise TypeError(
+            f'If specified, threshold must be of type float. Got: {type(threshold)}')
 
     # peaks over threshold
     x = np.sort(x)
@@ -370,9 +380,10 @@ def ste_peaks(peaks_distribution, npeaks):
     ste: scipy.stats.rv_frozen
             Short-term extreme distribution.
     """
-    assert callable(peaks_distribution.cdf
-                    ), 'peaks_distribution must be a scipy.stat distribution.'
-    assert isinstance(npeaks, float), 'npeaks must be of type float'
+    if not callable(peaks_distribution.cdf):
+        raise TypeError('peaks_distribution must be a scipy.stat distribution.')
+    if not isinstance(npeaks, float):
+        raise TypeError(f'npeaks must be of type float. Got: {type(npeaks)}')
 
     class _ShortTermExtreme(stats.rv_continuous):
 
@@ -415,9 +426,12 @@ def block_maxima(t, x, t_st):
     block_maxima: np.array
         Block maxima (i.e. largest peak in each block).
     """
-    assert isinstance(t, np.ndarray), 't must be of type np.ndarray'
-    assert isinstance(x, np.ndarray), 'x must be of type np.ndarray'
-    assert isinstance(t_st, float), 't_st must be of type float'
+    if not isinstance(t, np.ndarray):
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
+    if not isinstance(x, np.ndarray):
+        raise TypeError(f'x must be of type np.ndarray. Got: {type(x)}')
+    if not isinstance(t_st, float):
+        raise TypeError(f't_st must be of type float. Got: {type(t_st)}')
 
     nblock = int(t[-1] / t_st)
     block_maxima = np.zeros(int(nblock))
@@ -442,8 +456,9 @@ def ste_block_maxima_gev(block_maxima):
     ste: scipy.stats.rv_frozen
             Short-term extreme distribution.
     """
-    assert isinstance(
-        block_maxima, np.ndarray), 'block_maxima must be of type np.ndarray'
+    if not isinstance(block_maxima, np.ndarray):
+        raise TypeError(
+            f'block_maxima must be of type np.ndarray. Got: {type(block_maxima)}')
 
     ste_params = stats.genextreme.fit(block_maxima)
     param_names = ['c', 'loc', 'scale']
@@ -468,8 +483,9 @@ def ste_block_maxima_gumbel(block_maxima):
     ste: scipy.stats.rv_frozen
             Short-term extreme distribution.
     """
-    assert isinstance(
-        block_maxima, np.ndarray), 'block_maxima must be of type np.ndarray'
+    if not isinstance(block_maxima, np.ndarray):
+        raise TypeError(
+            f'block_maxima must be of type np.ndarray. Got: {type(block_maxima)}')
 
     ste_params = stats.gumbel_r.fit(block_maxima)
     param_names = ['loc', 'scale']
@@ -514,10 +530,14 @@ def short_term_extreme(t, data, t_st, method):
     ste: scipy.stats.rv_frozen
             Short-term extreme distribution.
     """
-    assert isinstance(t, np.ndarray), 't must be of type np.ndarray'
-    assert isinstance(data, np.ndarray), 'x must be of type np.ndarray'
-    assert isinstance(t_st, float), 't_st must be of type float'
-    assert isinstance(method, str), 'method must be of type string'
+    if not isinstance(t, np.ndarray):
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
+    if not isinstance(t_st, float):
+        raise TypeError(f't_st must be of type float. Got: {type(t_st)}')
+    if not isinstance(method, str):
+        raise TypeError(f'method must be of type string. Got: {type(method)}')
 
     peaks_methods = {
         'peaks_weibull': peaks_distribution_weibull,
@@ -555,7 +575,7 @@ def full_seastate_long_term_extreme(ste, weights):
     ste: list[scipy.stats.rv_frozen]
         Short-term extreme distribution of the quantity of interest for
         each sample sea state.
-    weights: list[floats]
+    weights: list, np.ndarray
         The weights from the full sea state sampling
 
     Returns
@@ -563,10 +583,12 @@ def full_seastate_long_term_extreme(ste, weights):
     ste: scipy.stats.rv_frozen
         Short-term extreme distribution.
     """
-    assert isinstance(
-        ste, list), 'ste must be of type list[scipy.stats.rv_frozen]'
-    assert isinstance(weights, (list, np.ndarray)
-                      ), 'weights must be of type list[floats]'
+    if not isinstance(ste, list):
+        raise TypeError(
+            f'ste must be of type list[scipy.stats.rv_frozen]. Got: {type(ste)}')
+    if not isinstance(weights, (list, np.ndarray)):
+        raise TypeError(
+            f'weights must be of type list or np.ndarray. Got: {type(weights)}')
 
     class _LongTermExtreme(stats.rv_continuous):
 
@@ -612,11 +634,16 @@ def mler_coefficients(rao, wave_spectrum, response_desired):
         rao = np.array(rao)
     except:
         pass
-    assert isinstance(rao, np.ndarray), 'rao must be of type np.ndarray'
-    assert isinstance(wave_spectrum, pd.DataFrame
-                      ), 'wave_spectrum must be of type pd.DataFrame'
-    assert isinstance(response_desired, (int, float)
-                      ), 'response_desired must be of type int or float'
+    
+    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):
+        raise TypeError(
+            f'wave_spectrum must be of type pd.DataFrame. 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
     # convert from Hz to rad/s
@@ -707,8 +734,9 @@ def mler_simulation(parameters=None):
         Simulation parameters including spatial and time calculated
         arrays.
     """
-    if not parameters == None:
-        assert isinstance(parameters, dict), 'parameters must be of type dict'
+    if not isinstance(parameters, (type(None), dict)):
+        raise TypeError(
+            f'If specified, parameters must be of type dict. Got: {type(parameters)}')
 
     sim = {}
 
@@ -762,11 +790,18 @@ def mler_wave_amp_normalize(wave_amp, mler, sim, k):
         k = np.array(k)
     except:
         pass
-    assert isinstance(mler, pd.DataFrame), 'mler must be of type pd.DataFrame'
-    assert isinstance(wave_amp, (int, float)
-                      ), 'wave_amp must be of type int or float'
-    assert isinstance(sim, dict), 'sim must be of type dict'
-    assert isinstance(k, np.ndarray), 'k must be of type ndarray'
+    if not isinstance(mler, pd.DataFrame):
+        raise TypeError(
+            f'mler must be of type pd.DataFrame. 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)}')
+    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)}')
 
     freq = mler.index.values * 2*np.pi
     dw = (max(freq) - min(freq)) / (len(freq)-1)  # get delta
@@ -826,10 +861,18 @@ def mler_export_time_series(rao, mler, sim, k):
         k = np.array(k)
     except:
         pass
-    assert isinstance(rao, np.ndarray), 'rao must be of type ndarray'
-    assert isinstance(mler, pd.DataFrame), 'mler must be of type pd.DataFrame'
-    assert isinstance(sim, dict), 'sim must be of type dict'
-    assert isinstance(k, np.ndarray), 'k must be of type ndarray'
+    if not isinstance(rao, np.ndarray):
+        raise TypeError(
+            f'rao must be of type ndarray. Got: {type(rao)}')
+    if not isinstance(mler, pd.DataFrame):
+        raise TypeError(
+            f'mler must be of type pd.DataFrame. 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)}')
 
     freq = mler.index.values * 2*np.pi  # convert Hz to rad/s
     dw = (max(freq) - min(freq)) / (len(freq)-1)  # get delta
@@ -874,9 +917,14 @@ def return_year_value(ppf, return_year, short_term_period_hr):
     value: float
         The value corresponding to the return period from the distribution.
     """
-    assert callable(ppf)
-    assert isinstance(return_year, (float, int))
-    assert isinstance(short_term_period_hr, (float, int))
+    if not callable(ppf):
+        raise TypeError('ppf must be a callable Percentage Point Function')
+    if not isinstance(return_year, (float, int)):
+        raise TypeError(
+            f'return_year must be of type float or int. Got: {type(return_year)}')
+    if not isinstance(short_term_period_hr, (float, int)):
+        raise TypeError(
+            f'short_term_period_hr must be of type float or int. Got: {type(short_term_period_hr)}')
 
     p = 1 / (return_year * 365.25 * 24 / short_term_period_hr)
 
diff --git a/mhkit/loads/general.py b/mhkit/loads/general.py
index 0c38b6bc7..28558eb58 100644
--- a/mhkit/loads/general.py
+++ b/mhkit/loads/general.py
@@ -27,17 +27,27 @@ def bin_statistics(data,bin_against,bin_edges,data_signal=[]):
         Standard deviation of each bim
     """
 
-    assert isinstance(data, pd.DataFrame), 'data must be of type pd.DataFram'   
-    try: bin_against = np.asarray(bin_against) 
-    except: 'bin_against must be of type np.ndarray'
-    try: bin_edges = np.asarray(bin_edges)
-    except: 'bin_edges must be of type np.ndarray'    
+    if not isinstance(data, pd.DataFrame):
+        raise TypeError(
+            f'data must be of type pd.DataFrame. Got: {type(data)}')
+    try:
+        bin_against = np.asarray(bin_against) 
+    except:
+        raise TypeError(
+            f'bin_against must be of type np.ndarray. Got: {type(bin_against)}')
+    try:
+        bin_edges = np.asarray(bin_edges)
+    except:
+        raise TypeError(
+            f'bin_edges must be of type np.ndarray. Got: {type(bin_edges)}')
 
     # Determine variables to analyze
     if len(data_signal)==0: # if not specified, bin all variables
         data_signal=data.columns.values
     else:
-        assert isinstance(data_signal, list), 'must be of type list'
+        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 = []
@@ -97,15 +107,28 @@ def blade_moments(blade_coefficients,flap_offset,flap_raw,edge_offset,edge_raw):
         Blade edgewise moment in SI units
     '''
     
-    try: blade_coefficients = np.asarray(blade_coefficients)
-    except: 'blade_coefficients must be of type np.ndarray'
-    try: flap_raw = np.asarray(flap_raw)
-    except: 'flap_raw must be of type np.ndarray'    
-    try: edge_raw = np.asarray(edge_raw)
-    except:  'edge_raw must be of type np.ndarray'    
+    try:
+        blade_coefficients = np.asarray(blade_coefficients)
+    except:
+        raise TypeError(
+            f'blade_coefficients must be of type np.ndarray. Got: {type(blade_coefficients)}')
+    try:
+        flap_raw = np.asarray(flap_raw)
+    except:
+        raise TypeError(
+            f'flap_raw must be of type np.ndarray. Got: {type(flap_raw)}')
+    try:
+        edge_raw = np.asarray(edge_raw)
+    except:
+        raise TypeError(
+            f'edge_raw must be of type np.ndarray. Got: {type(edge_raw)}')
     
-    assert isinstance(flap_offset, (float,int)), 'flap_offset must be of type int or float'
-    assert isinstance(edge_offset, (float,int)), 'edge_offset must be of type int or float'
+    if not isinstance(flap_offset, (float,int)):
+        raise TypeError(
+            f'flap_offset must be of type int or float. Got: {type(flap_offset)}')
+    if not isinstance(edge_offset, (float,int)):
+        raise TypeError(
+            f'edge_offset must be of type int or float. Got: {type(edge_offset)}')
     
     # remove offset from raw signal
     flap_signal = flap_raw - flap_offset
@@ -150,11 +173,19 @@ def damage_equivalent_load(data_signal, m, bin_num=100, data_length=600):
         Damage equivalent load (DEL) of single data signal
     '''
     
-    try: data_signal = np.array(data_signal)
-    except: 'data_signal must be of type np.ndarray'
-    assert isinstance(m, (float,int)), 'm must be of type float or int'
-    assert isinstance(bin_num, (float,int)), 'bin_num must be of type float or int'
-    assert isinstance(data_length, (float,int)), 'data_length must be of type float or int'
+    try:
+        data_signal = np.array(data_signal)
+    except:
+        raise TypeError(
+            f'data_signal must be of type np.ndarray. Got: {type(data_signal)}')
+    if not isinstance(m, (float,int)):
+        raise TypeError(f'm must be of type float or int. Got: {type(m)}')
+    if not isinstance(bin_num, (float,int)):
+        raise TypeError(
+            f'bin_num must be of type float or int. Got: {type(bin_num)}')
+    if not isinstance(data_length, (float,int)):
+        raise TypeError(
+            f'data_length must be of type float or int. Got: {type(data_length)}')
 
     rainflow_ranges = fatpack.find_rainflow_ranges(data_signal,k=256)
 
diff --git a/mhkit/loads/graphics.py b/mhkit/loads/graphics.py
index 291e028af..93a246196 100644
--- a/mhkit/loads/graphics.py
+++ b/mhkit/loads/graphics.py
@@ -32,24 +32,37 @@ def plot_statistics(x,y_mean,y_max,y_min,y_stdev=[],**kwargs):
     ax : matplotlib pyplot axes
     '''
     
-    try: x = np.array(x)
-    except: 'x must be of type np.ndarray'       
-    try: y_mean = np.array(y_mean)
-    except: 'y_mean must be of type np.ndarray'           
-    try:y_max = np.array(y_max)
-    except: 'y_max must be of type np.ndarray'
-    try: y_min = np.array(y_min)
-    except: 'y_min must be of type np.ndarray'
+    try:
+        x = np.array(x)
+    except:
+        raise TypeError(f'x must be of type np.ndarray. Got: {type(x)}')
+    try:
+        y_mean = np.array(y_mean)
+    except:
+        raise TypeError(f'y_mean must be of type np.ndarray. Got: {type(y_mean)}')       
+    try:
+        y_max = np.array(y_max)
+    except:
+        raise TypeError(f'y_max must be of type np.ndarray. Got: {type(y_max)}')
+    try:
+        y_min = np.array(y_min)
+    except:
+        raise TypeError(f'y_min must be of type np.ndarray. Got: {type(y_min)}')
     
     x_label   = kwargs.get("x_label", None)
     y_label   = kwargs.get("y_label", None)
     title     = kwargs.get("title", None)
     save_path = kwargs.get("save_path", None)
     
-    assert isinstance(x_label, (str, type(None))), 'x_label must be of type str'
-    assert isinstance(y_label, (str, type(None))), 'y_label must be of type str'
-    assert isinstance(title, (str, type(None))), 'title must be of type str'
-    assert isinstance(save_path, (str, type(None))), 'save_path must be of type str'
+    if not isinstance(x_label, (str, type(None))):
+        raise TypeError(f'x_label must be of type str. Got: {type(x_label)}')
+    if not isinstance(y_label, (str, type(None))):
+        raise TypeError(f'y_label must be of type str. Got: {type(y_label)}')
+    if not isinstance(title, (str, type(None))):
+        raise TypeError(f'title must be of type str. Got: {type(title)}')
+    if not isinstance(save_path, (str, type(None))):
+        raise TypeError(
+            f'save_path must be of type str. Got: {type(save_path)}')
 
     fig, ax = plt.subplots(figsize=(6,4))
     ax.plot(x,y_max,'^',label='max',mfc='none')
@@ -132,10 +145,15 @@ def plot_bin_statistics(bin_centers, bin_mean,bin_max, bin_min,
     title     = kwargs.get("title", None)
     save_path = kwargs.get("save_path", None)
     
-    assert isinstance(x_label, (str, type(None))), 'x_label must be of type str'
-    assert isinstance(y_label, (str, type(None))), 'y_label must be of type str'
-    assert isinstance(title, (str, type(None))), 'title must be of type str'
-    assert isinstance(save_path, (str, type(None))), 'save_path must be of type str'
+    if not isinstance(x_label, (str, type(None))):
+        raise TypeError(f'x_label must be of type str. Got: {type(x_label)}')
+    if not isinstance(y_label, (str, type(None))):
+        raise TypeError(f'y_label must be of type str. Got: {type(y_label)}')
+    if not isinstance(title, (str, type(None))):
+        raise TypeError(f'title must be of type str. Got: {type(title)}')
+    if not isinstance(save_path, (str, type(None))):
+        raise TypeError(
+            f'save_path must be of type str. Got: {type(save_path)}')
     
     fig, ax = plt.subplots(figsize=(7,5))
     ax.errorbar(bin_centers,bin_max,marker='^',mfc='none',
diff --git a/mhkit/power/characteristics.py b/mhkit/power/characteristics.py
index 08578f984..0569db572 100644
--- a/mhkit/power/characteristics.py
+++ b/mhkit/power/characteristics.py
@@ -21,7 +21,8 @@ def instantaneous_frequency(um):
         Frequency of the measured voltage (Hz) indexed by time  
         with signal name columns
     """  
-    assert isinstance(um, (pd.Series, pd.DataFrame)), 'um must be of type pd.Series or pd.DataFrame'
+    if not isinstance(um, (pd.Series, pd.DataFrame)):
+        raise TypeError(f'um must be of type pd.Series or pd.DataFrame. Got: {type(um)}')
    
     if isinstance(um.index[0], datetime.datetime):
         t = (um.index - datetime.datetime(1970,1,1)).total_seconds()
@@ -59,9 +60,12 @@ def dc_power(voltage, current):
     P: pandas DataFrame
         DC power [W] from each channel and gross power indexed by time
     """
-    assert isinstance(voltage, (pd.Series, pd.DataFrame)), 'voltage must be of type pd.Series or pd.DataFrame'
-    assert isinstance(current, (pd.Series, pd.DataFrame)), 'current must be of type pd.Series or pd.DataFrame'
-    assert voltage.shape == current.shape, 'current and volatge must have the same shape'
+    if not isinstance(voltage, (pd.Series, pd.DataFrame)):
+        raise TypeError(f'voltage must be of type pd.Series or pd.DataFrame. Got: {type(voltage)}')
+    if not isinstance(current, (pd.Series, pd.DataFrame)):
+        raise TypeError(f'current must be of type pd.Series or pd.DataFrame. Got: {type(current)}')
+    if not voltage.shape == current.shape:
+        raise ValueError('current and volatge must have the same shape')
     
     
     P = current.values * voltage.values
@@ -90,11 +94,16 @@ def ac_power_three_phase(voltage, current, power_factor, line_to_line=False):
     P: pandas DataFrame
         Magnitude of active AC power [W] indexed by time with Power column 
     """
-    assert isinstance(voltage, pd.DataFrame), 'voltage must be of type pd.DataFrame'
-    assert isinstance(current, pd.DataFrame), 'current must be of type pd.DataFrame'
-    assert len(voltage.columns) == 3, 'voltage must have three columns'
-    assert len(current.columns) == 3, 'current must have three columns'
-    assert current.shape == voltage.shape, 'current and voltage must be of the same size'
+    if not isinstance(voltage, pd.DataFrame):
+        raise TypeError(f'voltage must be of type pd.DataFrame. Got: {type(voltage)}')
+    if not isinstance(current, pd.DataFrame):
+        raise TypeError(f'current must be of type pd.DataFrame. Got: {type(current)}')
+    if not len(voltage.columns) == 3:
+        raise ValueError('voltage must have three columns')
+    if not len(current.columns) == 3:
+        raise ValueError('current must have three columns')
+    if not current.shape == voltage.shape:
+        raise ValueError('current and voltage must be of the same size')
     
 
     abs_current = np.abs(current.values)
diff --git a/mhkit/river/io/usgs.py b/mhkit/river/io/usgs.py
index 4282f6f56..4a69e2dd2 100644
--- a/mhkit/river/io/usgs.py
+++ b/mhkit/river/io/usgs.py
@@ -92,8 +92,8 @@ def request_usgs_data(
         Data indexed by datetime with columns named according to the parameter's 
         variable description
     """
-    assert data_type in [
-        'Daily', 'Instantaneous'], 'data_type must be Daily or Instantaneous'
+    if not data_type in ['Daily', 'Instantaneous']:
+        raise ValueError(f'data_type must be Daily or Instantaneous. Got: {data_type}')
 
     # Define the path to the cache directory
     cache_dir = os.path.join(os.path.expanduser("~"),
diff --git a/mhkit/river/performance.py b/mhkit/river/performance.py
index ac4529793..ac28393c3 100644
--- a/mhkit/river/performance.py
+++ b/mhkit/river/performance.py
@@ -17,7 +17,8 @@ def circular(diameter):
     projected_capture_area : float
         Projected capture area [m^2]
     """
-    assert isinstance(diameter, (int,float)), 'diameter must be of type int or float'
+    if not isinstance(diameter, (int,float)):
+        raise TypeError(f'diameter must be of type int or float. Got: {type(diameter)}')
     
     equivalent_diameter = diameter
     projected_capture_area = (1/4)*np.pi*(equivalent_diameter**2) 
@@ -41,7 +42,8 @@ def ducted(duct_diameter):
     projected_capture_area : float
         Projected capture area [m^2]
     """
-    assert isinstance(duct_diameter, (int,float)), 'duct_diameter must be of type int or float'
+    if not isinstance(duct_diameter, (int,float)):
+        raise TypeError(f'duct_diameter must be of type int or float. Got: {type(duct_diameter)}')
     
     equivalent_diameter = duct_diameter
     projected_capture_area = (1/4)*np.pi*(equivalent_diameter**2) 
@@ -67,8 +69,10 @@ def rectangular(h, w):
     projected_capture_area : float
         Projected capture area [m^2]
     """
-    assert isinstance(h, (int,float)), 'h must be of type int or float'
-    assert isinstance(w, (int,float)), 'w must be of type int or float'
+    if not isinstance(h, (int,float)):
+        raise TypeError(f'h must be of type int or float. Got: {type(h)}')
+    if not isinstance(w, (int,float)):
+        raise TypeError(f'w must be of type int or float. Got: {type(w)}')
     
     equivalent_diameter = np.sqrt(4.*h*w / np.pi) 
     projected_capture_area = h*w
@@ -92,7 +96,8 @@ def multiple_circular(diameters):
     projected_capture_area : float
         Projected capture area [m^2]
     """
-    assert isinstance(diameters, list), 'diameters must be of type list'
+    if not isinstance(diameters, list):
+        raise TypeError(f'diameters must be of type list. Got: {type(diameters)}')
     
     diameters_squared = [x**2 for x in diameters]
     equivalent_diameter = np.sqrt(sum(diameters_squared))
@@ -124,7 +129,8 @@ def tip_speed_ratio(rotor_speed,rotor_diameter,inflow_speed):
     try: inflow_speed = np.asarray(inflow_speed)
     except: 'inflow_speed must be of type np.ndarray'
     
-    assert isinstance(rotor_diameter, (float,int)), 'rotor diameter must be of type int or float'
+    if not isinstance(rotor_diameter, (float,int)):
+        raise TypeError(f'rotor_diameter must be of type int or float. Got: {type(rotor_diameter)}')
 
 
     rotor_velocity = rotor_speed * np.pi*rotor_diameter
@@ -159,8 +165,10 @@ def power_coefficient(power,inflow_speed,capture_area,rho):
     try: inflow_speed = np.asarray(inflow_speed)
     except: 'inflow_speed must be of type np.ndarray'
     
-    assert isinstance(capture_area, (float,int)), 'capture_area must be of type int or float'
-    assert isinstance(rho, (float,int)), 'rho must be of type int or float'
+    if not isinstance(capture_area, (float,int)):
+        raise TypeError(f'capture_area must be of type int or float. Got: {type(capture_area)}')
+    if not isinstance(rho, (float,int)):
+        raise TypeError(f'rho must be of type int or float. Got: {type(rho)}')
 
     # Predicted power from inflow
     power_in = (0.5 * rho * capture_area * inflow_speed**3)
diff --git a/mhkit/river/resource.py b/mhkit/river/resource.py
index c4a0e760c..fcf2e0d07 100644
--- a/mhkit/river/resource.py
+++ b/mhkit/river/resource.py
@@ -14,7 +14,7 @@ def Froude_number(v, h, g=9.80665):
     v : int/float 
         Average velocity [m/s].
     h : int/float
-        Mean hydrolic depth float [m].
+        Mean hydraulic depth float [m].
     g : int/float
         Gravitational acceleration [m/s2].
 
@@ -24,9 +24,12 @@ def Froude_number(v, h, g=9.80665):
         Froude Number of the river [unitless].
 
     """
-    assert isinstance(v, (int,float)), 'v must be of type int or float'
-    assert isinstance(h, (int,float)), 'h must be of type int or float'
-    assert isinstance(g, (int,float)), 'g must be of type int or float'
+    if not isinstance(v, (int,float)):
+        raise TypeError(f'v must be of type int or float. Got: {type(v)}')
+    if not isinstance(h, (int,float)):
+        raise TypeError(f'h must be of type int or float. Got: {type(h)}')
+    if not isinstance(g, (int,float)):
+        raise TypeError(f'g must be of type int or float. Got: {type(g)}')
     
     Fr = v / np.sqrt( g * h )
     
@@ -46,8 +49,10 @@ def exceedance_probability(D):
     -------
     F : pandas DataFrame    
         Exceedance probability [unitless] indexed by time [datetime or s]
-    """  
-    assert isinstance(D, (pd.DataFrame, pd.Series)), 'D must be of type pd.Series' # dataframe allowed for matlab
+    """
+    # dataframe allowed for matlab
+    if not isinstance(D, (pd.DataFrame, pd.Series)):
+        raise TypeError(f'D must be of type pd.Series or pd.DataFrame. Got: {type(D)}')
     
     if isinstance(D, pd.DataFrame) and len(D.columns) == 1: # for matlab
         D = D.squeeze().copy()
@@ -91,9 +96,12 @@ def polynomial_fit(x, y, n):
         y = np.array(y)
     except:
         pass
-    assert isinstance(x, np.ndarray), 'x must be of type np.ndarray'
-    assert isinstance(y, np.ndarray), 'y must be of type np.ndarray'
-    assert isinstance(n, int), 'n must be of type int'
+    if not isinstance(x, np.ndarray):
+        raise TypeError(f'x must be of type np.ndarray. Got: {type(x)}')
+    if not isinstance(y, np.ndarray):
+        raise TypeError(f'y must be of type np.ndarray. Got: {type(y)}')
+    if not isinstance(n, int):
+        raise TypeError(f'n must be of type int. Got: {type(n)}')
     
     # Get coeffcients of polynomial of order n 
     polynomial_coefficients = np.poly1d(np.polyfit(x, y, n))
@@ -122,9 +130,12 @@ def discharge_to_velocity(D, polynomial_coefficients):
     ------------
     V: pandas DataFrame   
         Velocity [m/s] indexed by time [datetime or s]
-    """  
-    assert isinstance(D, (pd.DataFrame, pd.Series)), 'D must be of type pd.Series' # dataframe allowed for matlab
-    assert isinstance(polynomial_coefficients, np.poly1d), 'polynomial_coefficients must be of type np.poly1d'
+    """
+    # dataframe allowed for matlab
+    if not isinstance(D, (pd.DataFrame, pd.Series)):
+        raise TypeError(f'D must be of type pd.Series. Got: {type(D)}')
+    if not isinstance(polynomial_coefficients, np.poly1d):
+        raise TypeError(f'polynomial_coefficients must be of type np.poly1d. Got: {type(polynomial_coefficients)}')
     
     if isinstance(D, pd.DataFrame) and len(D.columns) == 1: # for matlab
         D = D.squeeze().copy()
@@ -159,11 +170,16 @@ def velocity_to_power(V, polynomial_coefficients, cut_in, cut_out):
     -------
     P : pandas DataFrame
         Power [W] indexed by time [datetime or s]
-    """  
-    assert isinstance(V, (pd.DataFrame, pd.Series)), 'V must be of type pd.Series' # dataframe allowed for matlab
-    assert isinstance(polynomial_coefficients, np.poly1d), 'polynomial_coefficients must be of type np.poly1d'
-    assert isinstance(cut_in, (int,float)), 'cut_in must be of type int or float'
-    assert isinstance(cut_out, (int,float)), 'cut_out must be of type int or float'
+    """
+    # dataframe allowed for matlab
+    if not isinstance(V, (pd.DataFrame, pd.Series)):
+        raise TypeError(f'V must be of type pd.Series or pd.DataFrame. Got: {type(V)}')
+    if not isinstance(polynomial_coefficients, np.poly1d):
+        raise TypeError(f'polynomial_coefficients must be of type np.poly1d. Got: {type(polynomial_coefficients)}')
+    if not isinstance(cut_in, (int,float)):
+        raise TypeError(f'cut_in must be of type int or float. Got: {type(cut_in)}')
+    if not isinstance(cut_out, (int,float)):
+        raise TypeError(f'cut_out must be of type int or float. Got: {type(cut_out)}')
     
     if isinstance(V, pd.DataFrame) and len(V.columns) == 1:
         V = V.squeeze().copy()
@@ -199,8 +215,11 @@ def energy_produced(P, seconds):
     E : float
         Energy [J] produced in the given time frame
     """
-    assert isinstance(P, (pd.DataFrame, pd.Series)), 'D must be of type pd.Series' # dataframe allowed for matlab
-    assert isinstance(seconds, (int, float)), 'seconds must be of type int or float' 
+    # dataframe allowed for matlab
+    if not isinstance(P, (pd.DataFrame, pd.Series)):
+        raise TypeError(f'P must be of type pd.Series or pd.DataFrame. Got: {type(P)}')
+    if not isinstance(seconds, (int, float)):
+        raise TypeError(f'seconds must be of type int or float. Got: {type(seconds)}')
 
     if isinstance(P, pd.DataFrame) and len(P.columns) == 1: # for matlab
         P = P.squeeze().copy()
diff --git a/mhkit/tidal/resource.py b/mhkit/tidal/resource.py
index ef9961780..3206001f4 100644
--- a/mhkit/tidal/resource.py
+++ b/mhkit/tidal/resource.py
@@ -93,8 +93,9 @@ def principal_flow_directions(directions, width_dir):
 
     if isinstance(directions, np.ndarray):
         directions=pd.Series(directions)
-    assert(all(directions>=0) and all(directions<=360),
-           'flood must be between 0 and 360 degrees')
+    if any(directions<0) or any(directions>360):
+        violating_values = [d for d in directions if d < 0 or d > 360]
+        raise ValueError(f'directions must be between 0 and 360 degrees. Values out of range: {violating_values}')
 
     # Number of directional bins 
     N_dir=int(360/width_dir)
diff --git a/mhkit/utils/__init__.py b/mhkit/utils/__init__.py
index 442d4f949..074232541 100644
--- a/mhkit/utils/__init__.py
+++ b/mhkit/utils/__init__.py
@@ -1,4 +1,4 @@
-from .time_utils import matlab_to_datetime, excel_to_datetime, index_to_datetime
+from .time_utils import matlab_to_datetime, excel_to_datetime
 from .stat_utils import get_statistics, vector_statistics, unwrap_vector, magnitude_phase, unorm
 from .cache import handle_caching, clear_cache
 from .upcrossing import upcrossing, peaks, troughs, heights, periods, custom
diff --git a/mhkit/utils/stat_utils.py b/mhkit/utils/stat_utils.py
index 92ed08fa9..639517de7 100644
--- a/mhkit/utils/stat_utils.py
+++ b/mhkit/utils/stat_utils.py
@@ -1,4 +1,3 @@
-import matplotlib.pyplot as plt
 from mhkit import qc
 import pandas as pd
 import numpy as np
@@ -28,15 +27,17 @@ def get_statistics(data, freq, period=600, vector_channels=[]):
         Calculated statistical values from the data, indexed by the first timestamp
     """
     # Check data type
-    assert isinstance(data, pd.DataFrame), 'data must be of type pd.DataFrame'
-    assert isinstance(freq, (float, int)), 'freq must be of type int or float'
-    assert isinstance(period, (float, int)
-                      ), 'freq must be of type int or float'
+    if not isinstance(data, pd.DataFrame):
+        raise TypeError(f'data must be of type pd.DataFrame. Got: {type(data)}')
+    if not isinstance(freq, (float, int)):
+        raise TypeError(f'freq must be of type int or float. Got: {type(freq)}')
+    if not isinstance(period, (float, int)):
+        raise TypeError(f'period must be of type int or float. Got: {type(period)}')
     # catch if vector_channels is not an string array
     if isinstance(vector_channels, str):
         vector_channels = [vector_channels]
-    assert isinstance(
-        vector_channels, list), 'vector_channels must be a list of strings'
+    if not isinstance(vector_channels, list):
+        raise TypeError(f'vector_channels must be a list of strings. Got: {type(vector_channels)}')
 
     # Check timestamp using qc module
     data.index = data.index.round('1ms')
@@ -110,7 +111,8 @@ def vector_statistics(data):
         data = np.array(data)
     except:
         pass
-    assert isinstance(data, np.ndarray), 'data must be of type np.ndarray'
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
 
     # calculate mean
     Ux = sum(np.sin(data*np.pi/180))/len(data)
@@ -152,7 +154,8 @@ def unwrap_vector(data):
         data = np.array(data)
     except:
         pass
-    assert isinstance(data, np.ndarray), 'data must be of type np.ndarray'
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
 
     # Loop through and unwrap points
     for i in range(len(data)):
@@ -196,9 +199,12 @@ def magnitude_phase(x, y, z=None):
         z = np.array(z)
         threeD = True
 
-    assert isinstance(x, (float, int, np.ndarray))
-    assert isinstance(y, (float, int, np.ndarray))
-    assert isinstance(z, (type(None), float, int, np.ndarray))
+    if not isinstance(x, (float, int, np.ndarray)):
+        raise TypeError(f'x must be of type float, int, or np.ndarray. Got: {type(x)}')
+    if not isinstance(y, (float, int, np.ndarray)):
+        raise TypeError(f'y must be of type float, int, or np.ndarray. Got: {type(y)}')
+    if not isinstance(z, (type(None), float, int, np.ndarray)):
+        raise TypeError(f'If specified, z must be of type float, int, or np.ndarray. Got: {type(z)}')
 
     if threeD:
         mag = np.sqrt(x**2 + y**2 + z**2)
@@ -236,14 +242,14 @@ def unorm(x, y, z):
     The resulting output is [ 8.1240384,  9.64365076, 11.22497216].
     '''
 
-    assert isinstance(x, (np.ndarray, np.float64, pd.Series)
-                      ), 'x must be an array'
-    assert isinstance(y, (np.ndarray, np.float64, pd.Series)
-                      ), 'y must be an array'
-    assert isinstance(z, (np.ndarray, np.float64, pd.Series)
-                      ), 'z must be an array'
-    assert all([len(x) == len(y), len(y) == len(z)]), ('lengths of arrays must'
-                                                       + ' match')
+    if not isinstance(x, (np.ndarray, np.float64, pd.Series)):
+        raise TypeError(f'x must be of type np.ndarray, np.float64, or pd.Series. Got: {type(x)}')
+    if not isinstance(y, (np.ndarray, np.float64, pd.Series)):
+        raise TypeError(f'y must be of type np.ndarray, np.float64, or pd.Series. Got: {type(y)}')
+    if not isinstance(z, (np.ndarray, np.float64, pd.Series)):
+        raise TypeError(f'z must be of type np.ndarray, np.float64, or pd.Series. Got: {type(z)}')
+    if not all([len(x) == len(y), len(y) == len(z)]):
+        raise ValueError('lengths of arrays must match')
 
     xyz = np.array([x, y, z])
     unorm = np.linalg.norm(xyz, axis=0)
diff --git a/mhkit/utils/time_utils.py b/mhkit/utils/time_utils.py
index 0a380ebea..b6b7ea494 100644
--- a/mhkit/utils/time_utils.py
+++ b/mhkit/utils/time_utils.py
@@ -1,9 +1,6 @@
-
-
 import datetime as dt
 import pandas as pd
 import numpy as np
-from pecos.utils import index_to_datetime
 
 
 def matlab_to_datetime(matlab_datenum):
@@ -25,8 +22,8 @@ def matlab_to_datetime(matlab_datenum):
         matlab_datenum = np.array(matlab_datenum, ndmin=1)
     except:
         pass
-    assert isinstance(
-        matlab_datenum, np.ndarray), 'data must be of type np.ndarray'
+    if not isinstance(matlab_datenum, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
 
     # Pre-allocate
     time = []
@@ -60,7 +57,8 @@ def excel_to_datetime(excel_num):
         excel_num = np.array(excel_num)
     except:
         pass
-    assert isinstance(excel_num, np.ndarray), 'data must be of type np.ndarray'
+    if not isinstance(excel_num, np.ndarray):
+        raise TypeError(f'excel_num must be of type np.ndarray. Got: {type(excel_num)}')
 
     # Convert to datetime
     time = pd.to_datetime('1899-12-30')+pd.to_timedelta(excel_num, 'D')
diff --git a/mhkit/utils/upcrossing.py b/mhkit/utils/upcrossing.py
index 5a244124f..24762a946 100644
--- a/mhkit/utils/upcrossing.py
+++ b/mhkit/utils/upcrossing.py
@@ -65,12 +65,11 @@ def upcrossing(t, data):
     inds: np.array
         Zero crossing indices
     """
+    # Check data types
     if not isinstance(t, np.ndarray):
-        raise TypeError('t must be of type np.ndarray')
-    
-    if not isinstance(t, np.ndarray):
-        'data must be of type np.ndarray'
-    
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
     if len(data.shape) != 1:
         raise ValueError('only 1D data supported, try calling squeeze()')
 
@@ -108,11 +107,11 @@ def peaks(t, data, inds=None):
         Peak values of the time-series
 
     """
+    # Check data types
     if not isinstance(t, np.ndarray):
-        TypeError('t must be of type np.ndarray')
-        
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
     if not isinstance(data, np.ndarray):
-        TypeError('data must be of type np.ndarray')
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
 
     return _apply(t, data, lambda ind1, ind2: np.max(data[ind1:ind2]), inds)
 
@@ -139,8 +138,11 @@ def troughs(t, data, inds=None):
         Trough values of the time-series
 
     """
-    assert isinstance(t, np.ndarray), 't must be of type np.ndarray'
-    assert isinstance(data, np.ndarray), 'data must be of type np.ndarray'
+    # Check data types
+    if not isinstance(t, np.ndarray):
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
 
     return _apply(t, data, lambda ind1, ind2: np.min(data[ind1:ind2]), inds)
 
@@ -169,8 +171,11 @@ def heights(t, data, inds=None):
     heights: np.array
         Height values of the time-series
     """
-    assert isinstance(t, np.ndarray), 't must be of type np.ndarray'
-    assert isinstance(data, np.ndarray), 'data must be of type np.ndarray'
+    # Check data types
+    if not isinstance(t, np.ndarray):
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
 
     def func(ind1, ind2): 
         return np.max(data[ind1:ind2]) - np.min(data[ind1:ind2])
@@ -199,8 +204,11 @@ def periods(t, data, inds=None):
     periods: np.array
         Period values of the time-series
     """
-    assert isinstance(t, np.ndarray), 't must be of type np.ndarray'
-    assert isinstance(data, np.ndarray), 'data must be of type np.ndarray'
+    # Check data types
+    if not isinstance(t, np.ndarray):
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
 
     return _apply(t, data, lambda ind1, ind2: t[ind2] - t[ind1], inds)
 
@@ -230,9 +238,12 @@ def custom(t, data, func, inds=None):
     values: np.array
         Custom values of the time-series
     """
-
-    assert isinstance(t, np.ndarray), 't must be of type np.ndarray'
-    assert isinstance(data, np.ndarray), 'data must of type np.ndarray'
-    assert callable(func), 'f must be callable'
+    # Check data types
+    if not isinstance(t, np.ndarray):
+        raise TypeError(f't must be of type np.ndarray. Got: {type(t)}')
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
+    if not callable(func):
+        raise ValueError('func must be callable')
 
     return _apply(t, data, func, inds)
diff --git a/mhkit/wave/contours.py b/mhkit/wave/contours.py
index 39a70a640..3007448b8 100644
--- a/mhkit/wave/contours.py
+++ b/mhkit/wave/contours.py
@@ -11,8 +11,6 @@
 mpl_version = tuple(map(int, matplotlib.__version__.split('.')))
 
 # Contours
-
-
 def environmental_contours(x1, x2, sea_state_duration, return_period,
                            method, **kwargs):
     """
@@ -85,12 +83,14 @@ def environmental_contours(x1, x2, sea_state_duration, return_period,
         x2 = np.array(x2)
     except:
         pass
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(sea_state_duration, (int, float)), (
-        'sea_state_duration must be of type int or float')
-    assert isinstance(return_period, (int, float, np.ndarray)), (
-        'return_period must be of type int, float, or array')
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(sea_state_duration, (int, float)):
+        raise TypeError(f'sea_state_duration must be of type int or float. Got: {type(sea_state_duration)}')
+    if not isinstance(return_period, (int, float, np.ndarray)):
+        raise TypeError(f'return_period must be of type int, float, or np.ndarray. Got: {type(return_period)}')
 
     bin_val_size = kwargs.get("bin_val_size", 0.25)
     nb_steps = kwargs.get("nb_steps", 1000)
@@ -104,27 +104,29 @@ def environmental_contours(x1, x2, sea_state_duration, return_period,
     PCA_bin_size = kwargs.get("PCA_bin_size", 250)
     return_fit = kwargs.get("return_fit", False)
 
-    assert isinstance(PCA, (dict, type(None))), (
-        'If specified PCA must be a dict')
-    assert isinstance(PCA_bin_size, int), 'PCA_bin_size must be of type int'
-    assert isinstance(return_fit, bool), 'return_fit must be of type bool'
-    assert isinstance(bin_val_size, (int, float)), (
-        'bin_val_size must be of type int or float')
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
-    assert isinstance(min_bin_count, int), ('min_bin_count must be of '
-                                            + 'type int')
-    assert isinstance(initial_bin_max_val, (int, float)), (
-        'initial_bin_max_val must be of type int or float')
-    if bandwidth == None:
-        assert (not 'bivariate_KDE' in method), (
-            'Must specify keyword bandwidth with bivariate KDE method')
+    if not isinstance(PCA, (dict, type(None))):
+        raise TypeError(f'If specified, PCA must be a dict. Got: {type(PCA)}')
+    if not isinstance(PCA_bin_size, int):
+        raise TypeError(f'PCA_bin_size must be of type int. Got: {type(PCA_bin_size)}')
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'return_fit must be of type bool. Got: {type(return_fit)}')
+    if not isinstance(bin_val_size, (int, float)):
+        raise TypeError(f'bin_val_size must be of type int or float. Got: {type(bin_val_size)}')
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
+    if not isinstance(min_bin_count, int):
+        raise TypeError(f'min_bin_count must be of type int. Got: {type(min_bin_count)}')
+    if not isinstance(initial_bin_max_val, (int, float)):
+        raise TypeError(f'initial_bin_max_val must be of type int or float. Got: {type(initial_bin_max_val)}')
+    if 'bivariate_KDE' in method and bandwidth == None:
+        raise TypeError(f'Must specify keyword bandwidth with bivariate KDE method. Got: {type(bandwidth)}')
 
     if isinstance(method, str):
         method = [method]
-    assert (len(set(method)) == len(method)), (
-        'Can only pass a unique '
-        + 'method once per function call. Consider wrapping this '
-        + 'function in a for loop to investage variations on the same method')
+    if not (len(set(method)) == len(method)):
+        raise ValueError(f'Can only pass a unique '
+                        + 'method once per function call. Consider wrapping this '
+                        + 'function in a for loop to investage variations on the same method')
 
     method_class = {'PCA': 'parametric',
                     'gaussian': 'parametric',
@@ -311,16 +313,21 @@ def PCA_contour(x1, x2, fit, kwargs):
         x2 = np.array(x2)
     except:
         pass
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
 
     bin_size = kwargs.get("bin_size", 250)
     nb_steps = kwargs.get("nb_steps", 1000)
     return_fit = kwargs.get("return_fit", False)
 
-    assert isinstance(bin_size, int), 'bin_size must be of type int'
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
-    assert isinstance(return_fit, bool), 'return_fit must be of type bool'
+    if not isinstance(bin_size, int):
+        raise TypeError(f'bin_size must be of type int. Got: {type(bin_size)}')
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'return_fit must be of type bool. Got: {type(return_fit)}')
 
     if 'x1_fit' not in fit:
         pca_fit = _principal_component_analysis(x1, x2, bin_size=bin_size)
@@ -414,9 +421,13 @@ def _principal_component_analysis(x1, x2, bin_size=250):
        'mu_param'      : fit to _mu_fcn
        'sigma_param'   : fit to _sig_fits
     """
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(bin_size, int), 'bin_size must be of type int'
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(bin_size, int):
+        raise TypeError(f'bin_size must be of type int. Got: {type(bin_size)}')
+    
     # Step 0: Perform Standard PCA
     mean_location = 0
     x1_mean_centered = x1 - x1.mean(axis=0)
@@ -457,7 +468,7 @@ def _principal_component_analysis(x1, x2, bin_size=250):
     minimum_4_bins = np.floor(N*0.25)
     if bin_size > minimum_4_bins:
         bin_size = minimum_4_bins
-        msg = ('To allow for a minimum of 4 bins the bin size has been' +
+        msg = ('To allow for a minimum of 4 bins, the bin size has been' +
                f'set to {minimum_4_bins}')
         print(msg)
 
@@ -545,12 +556,12 @@ def _iso_prob_and_quantile(sea_state_duration, return_period, nb_steps):
         'y_quantile' - CDF of y-component
     """
 
-    assert isinstance(sea_state_duration, (int, float)
-                      ), 'sea_state_duration must be of type int or float'
-    assert isinstance(return_period, (int, float)), (
-        'return_period must be of type int or float')
-
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
+    if not isinstance(sea_state_duration, (int, float)):
+        raise TypeError(f'sea_state_duration must be of type int or float. Got: {type(sea_state_duration)}')
+    if not isinstance(return_period, (int, float)):
+        raise TypeError(f'return_period must be of type int or float. Got: {type(return_period)}')
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
 
     dt_yrs = sea_state_duration / (3600 * 24 * 365)
     exceedance_probability = 1 / (return_period / dt_yrs)
@@ -606,14 +617,16 @@ def _copula_parameters(x1, x2, min_bin_count, initial_bin_max_val,
     std_cond: array
         Estimate coefficients of the standard deviation of Ln(x2|x1)
     """
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(min_bin_count, int), ('min_bin_count must be of'
-                                            + 'type int')
-    assert isinstance(bin_val_size, (int, float)), (
-        'bin_val_size must be of type int or float')
-    assert isinstance(initial_bin_max_val, (int, float)), (
-        'initial_bin_max_val must be of type int or float')
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(min_bin_count, int):
+        raise TypeError(f'min_bin_count must be of type int. Got: {type(min_bin_count)}')
+    if not isinstance(bin_val_size, (int, float)):
+        raise TypeError(f'bin_val_size must be of type int or float. Got: {type(bin_val_size)}')
+    if not isinstance(initial_bin_max_val, (int, float)):
+        raise TypeError(f'initial_bin_max_val must be of type int or float. Got: {type(initial_bin_max_val)}')
 
     # Binning
     x1_sorted_index = x1.argsort()
@@ -757,14 +770,15 @@ def _gaussian_copula(x1, x2, fit, component_1, kwargs):
         x2 = np.array(x2)
     except:
         pass
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(component_1, np.ndarray), (
-        'x2 must be of type np.ndarray')
-
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(component_1, np.ndarray):
+        raise TypeError(f'component_1 must be of type np.ndarray. Got: {type(component_1)}')
     return_fit = kwargs.get("return_fit", False)
-    assert isinstance(return_fit, bool), (
-        'If specified return_fit must be a bool')
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'If specified, return_fit must be of type bool. Got: {type(return_fit)}')
 
     x_component_iso_prob = fit['x_component_iso_prob']
     y_component_iso_prob = fit['y_component_iso_prob']
@@ -873,13 +887,15 @@ def _gumbel_copula(x1, x2, fit, component_1, nb_steps, kwargs):
         x2 = np.array(x2)
     except:
         pass
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(component_1, np.ndarray), 'x2 must be of type np.ndarray'
-
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(component_1, np.ndarray):
+        raise TypeError(f'component_1 must be of type np.ndarray. Got: {type(component_1)}')
     return_fit = kwargs.get("return_fit", False)
-    assert isinstance(
-        return_fit, bool), 'If specified return_fit must be a bool'
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'If specified, return_fit must be of type bool. Got: {type(return_fit)}')
 
     x_quantile = fit['x_quantile']
     y_quantile = fit['y_quantile']
@@ -971,12 +987,15 @@ def _clayton_copula(x1, x2, fit, component_1, kwargs):
         If return_fit=True. Dictionary with iso-probabilities passed
         with additional fit metrics from the copula method.
     """
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(component_1, np.ndarray), 'x2 must be of type np.ndarray'
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(component_1, np.ndarray):
+        raise TypeError(f'component_1 must be of type np.ndarray. Got: {type(component_1)}')
     return_fit = kwargs.get("return_fit", False)
-    assert isinstance(
-        return_fit, bool), 'If specified return_fit must be a bool'
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'If specified, return_fit must be of type bool. Got: {type(return_fit)}')
 
     x_quantile = fit['x_quantile']
     y_quantile = fit['y_quantile']
@@ -1051,12 +1070,15 @@ def _rosenblatt_copula(x1, x2, fit, component_1, kwargs):
         x2 = np.array(x2)
     except:
         pass
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(component_1, np.ndarray), 'x2 must be of type np.ndarray'
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(component_1, np.ndarray):
+        raise TypeError(f'component_1 must be of type np.ndarray. Got: {type(component_1)}')
     return_fit = kwargs.get("return_fit", False)
-    assert isinstance(
-        return_fit, bool), 'If specified return_fit must be a bool'
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'If specified, return_fit must be of type bool. Got: {type(return_fit)}')
 
     y_quantile = fit['y_quantile']
     mean_cond = fit['mean_cond']
@@ -1106,15 +1128,20 @@ def _nonparametric_copula_parameters(x1, x2, max_x1=None, max_x2=None,
     nonpara_pdf_2:
         x2 points in KDE space and Nonparametric PDF for x2
     """
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
     if not max_x1:
         max_x1 = x1.max()*2
     if not max_x2:
         max_x2 = x2.max()*2
-    assert isinstance(max_x1, float), 'max_x1 must be of type float'
-    assert isinstance(max_x2, float), 'max_x2 must be of type float'
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
+    if not isinstance(max_x1, float):
+        raise TypeError(f'max_x1 must be of type float. Got: {type(max_x1)}')
+    if not isinstance(max_x2, float):
+        raise TypeError(f'max_x2 must be of type float. Got: {type(max_x2)}')
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
 
     # Binning
     x1_sorted_index = x1.argsort()
@@ -1180,7 +1207,8 @@ def _nonparametric_component(z, nonpara_dist, nb_steps):
     component: array
         nonparametic component values
     """
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
 
     component = np.zeros(nb_steps)
     for k in range(0, nb_steps):
@@ -1227,13 +1255,15 @@ def _nonparametric_gaussian_copula(x1, x2, fit, nb_steps, kwargs):
         If return_fit=True. Dictionary with iso-probabilities passed
         with additional fit metrics from the copula method.
     """
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
-
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
     return_fit = kwargs.get("return_fit", False)
-    assert isinstance(
-        return_fit, bool), 'If specified return_fit must be a bool'
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'If specified, return_fit must be of type bool. Got: {type(return_fit)}')
 
     x_component_iso_prob = fit['x_component_iso_prob']
     y_component_iso_prob = fit['y_component_iso_prob']
@@ -1306,13 +1336,15 @@ def _nonparametric_clayton_copula(x1, x2, fit, nb_steps, kwargs):
         If return_fit=True. Dictionary with iso-probabilities passed
         with additional fit metrics from the copula method.
     """
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
-
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
     return_fit = kwargs.get("return_fit", False)
-    assert isinstance(return_fit, bool), ('If specified return_fit '
-                                          + 'must be a bool')
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'If specified, return_fit must be of type bool. Got: {type(return_fit)}')
 
     x_component_iso_prob = fit['x_component_iso_prob']
     x_quantile = fit['x_quantile']
@@ -1388,13 +1420,15 @@ def _nonparametric_gumbel_copula(x1, x2, fit, nb_steps, kwargs):
         If return_fit=True. Dictionary with iso-probabilities passed
         with additional fit metrics from the copula method.
     """
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
-
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
     return_fit = kwargs.get("return_fit", False)
-    assert isinstance(return_fit, bool), ('If specified return_fit '
-                                          + 'must be a bool')
+    if not isinstance(return_fit, bool):
+        raise TypeError(f'If specified, return_fit must be a bool. Got: {type(return_fit)}')
 
     Ndata = 1000
 
@@ -1470,7 +1504,7 @@ def _bivariate_KDE(x1, x2, bw, fit, nb_steps, Ndata_bivariate_KDE, kwargs):
     fit: Dictionay
         Dictionary of the iso-probability results
     nb_steps: int
-        number of points used to discritize KDE space
+        number of points used to discretize KDE space
     max_x1: float
         Defines the max value of x1 to discretize the KDE space
     max_x2: float
@@ -1491,9 +1525,12 @@ def _bivariate_KDE(x1, x2, bw, fit, nb_steps, Ndata_bivariate_KDE, kwargs):
         If return_fit=True. Dictionary with iso-probabilities passed
         with additional fit metrics from the copula method.
     """
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(nb_steps, int), 'nb_steps must be of type int'
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(nb_steps, int):
+        raise TypeError(f'nb_steps must be of type int. Got: {type(nb_steps)}')
 
     max_x1 = kwargs.get("max_x1", None)
     max_x2 = kwargs.get("max_x2", None)
@@ -1504,12 +1541,14 @@ def _bivariate_KDE(x1, x2, bw, fit, nb_steps, Ndata_bivariate_KDE, kwargs):
         max_x1 = x1.max()*2
     if isinstance(max_x2, type(None)):
         max_x2 = x2.max()*2
-    assert isinstance(max_x1, float), 'max_x1 must be of type float'
-    assert isinstance(max_x2, float), 'max_x2 must be of type float'
-    assert isinstance(log_transform, bool), ('If specified log_transform'
-                                             + 'must be a bool')
-    assert isinstance(return_fit, bool), ('If specified return_fit must '
-                                          + 'be a bool')
+    if not isinstance(max_x1, float):
+        raise TypeError(f'max_x1 must be of type float. Got: {type(max_x1)}')
+    if not isinstance(max_x2, float):
+       raise TypeError(f'max_x2 must be of type float. Got: {type(max_x2)}')
+    if not isinstance(log_transform, bool):
+       raise TypeError(f'If specified, log_transform must be of type bool. Got: {type(log_transform)}')
+    if not isinstance(return_fit, bool):
+       raise TypeError(f'If specified, return_fit must be of type bool. Got: {type(return_fit)}')
 
     p_f = fit['exceedance_probability']
 
@@ -1625,17 +1664,20 @@ def samples_full_seastate(x1, x2, points_per_interval, return_periods,
         raise NotImplementedError(
             "Full sea state sampling is currently only implemented using " +
             "the 'PCA' method.")
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(points_per_interval,
-                      int), 'points_per_interval must be of int'
-    assert isinstance(return_periods, np.ndarray
-                      ), 'return_periods must be of type np.ndarray'
-    assert isinstance(sea_state_duration, (int, float)
-                      ), 'sea_state_duration must be of int or float'
-    assert isinstance(method, (str, list)
-                      ), 'method must be of type string or list'
-    assert isinstance(bin_size, int), 'bin_size must be of int'
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(points_per_interval,int):
+        raise TypeError(f'points_per_interval must be of int. Got: {type(points_per_interval)}')
+    if not isinstance(return_periods, np.ndarray):
+        raise TypeError(f'return_periods must be of type np.ndarray. Got: {type(return_periods)}')
+    if not isinstance(sea_state_duration, (int, float)):
+        raise TypeError(f'sea_state_duration must be of int or float. Got: {type(sea_state_duration)}')
+    if not isinstance(method, (str, list)):
+        raise TypeError(f'method must be of type string or list. Got: {type(method)}')
+    if not isinstance(bin_size, int):
+        raise TypeError(f'bin_size must be of int. Got: {type(bin_size)}')
 
     pca_fit = _principal_component_analysis(x1, x2, bin_size)
 
@@ -1744,12 +1786,12 @@ def samples_contour(t_samples, t_contour, hs_contour):
     hs_samples : nparray
         points sampled along return contour
     """
-    assert isinstance(
-        t_samples, np.ndarray), 't_samples must be of type np.ndarray'
-    assert isinstance(
-        t_contour, np.ndarray), 't_contour must be of type np.ndarray'
-    assert isinstance(
-        hs_contour, np.ndarray), 'hs_contour must be of type np.ndarray'
+    if not isinstance(t_samples, np.ndarray):
+        raise TypeError(f't_samples must be of type np.ndarray. Got: {type(t_samples)}')
+    if not isinstance(t_contour, np.ndarray):
+        raise TypeError(f't_contour must be of type np.ndarray. Got: {type(t_contour)}')
+    if not isinstance(hs_contour, np.ndarray):
+        raise TypeError(f'hs_contour must be of type np.ndarray. Got: {type(hs_contour)}')
 
     # finds minimum and maximum energy period values
     amin = np.argmin(t_contour)
@@ -1801,16 +1843,16 @@ def _generate_sample_data(beta_lines, rho_zeroline, theta_zeroline,
     weight_points: np.array
         Array of weights for each point.
     """
-    assert isinstance(
-        beta_lines, np.ndarray), 'beta_lines must be of type np.ndarray'
-    assert isinstance(
-        rho_zeroline, np.ndarray), 'rho_zeroline must be of type np.ndarray'
-    assert isinstance(theta_zeroline, np.ndarray
-                      ), 'theta_zeroline must be of type np.ndarray'
-    assert isinstance(points_per_interval, int
-                      ), 'points_per_interval must be of type int'
-    assert isinstance(
-        contour_probs, np.ndarray), 'contour_probs must be of type np.ndarray'
+    if not isinstance(beta_lines, np.ndarray):
+        raise TypeError(f'beta_lines must be of type np.ndarray. Got: {type(beta_lines)}')
+    if not isinstance(rho_zeroline, np.ndarray):
+        raise TypeError(f'rho_zeroline must be of type np.ndarray. Got: {type(rho_zeroline)}')
+    if not isinstance(theta_zeroline, np.ndarray):
+        raise TypeError(f'theta_zeroline must be of type np.ndarray. Got: {type(theta_zeroline)}')
+    if not isinstance(points_per_interval, int):
+        raise TypeError(f'points_per_interval must be of type int. Got: {type(points_per_interval)}')
+    if not isinstance(contour_probs, np.ndarray):
+        raise TypeError(f'contour_probs must be of type np.ndarray. Got: {type(contour_probs)}')
 
     num_samples = (len(beta_lines) - 1) * points_per_interval
     alpha_bounds = np.zeros((len(beta_lines) - 1, 2))
@@ -1889,12 +1931,14 @@ def _princomp_inv(princip_data1, princip_data2, coeff, shift):
     original2: np.array
         T values following rotation from principal component space.
     """
-    assert isinstance(
-        princip_data1, np.ndarray), 'princip_data1 must be of type np.ndarray'
-    assert isinstance(
-        princip_data2, np.ndarray), 'princip_data2 must be of type np.ndarray'
-    assert isinstance(coeff, np.ndarray), 'coeff must be of type np.ndarray'
-    assert isinstance(shift, float), 'float must be of type float'
+    if not isinstance(princip_data1, np.ndarray):
+        raise TypeError(f'princip_data1 must be of type np.ndarray. Got: {type(princip_data1)}')
+    if not isinstance(princip_data2, np.ndarray):
+        raise TypeError(f'princip_data2 must be of type np.ndarray. Got: {type(princip_data2)}')
+    if not isinstance(coeff, np.ndarray):
+        raise TypeError(f'coeff must be of type np.ndarray. Got: {type(coeff)}')
+    if not isinstance(shift, float):
+        raise TypeError(f'shift must be of type float. Got: {type(shift)}')
 
     original1 = np.zeros(len(princip_data1))
     original2 = np.zeros(len(princip_data1))
diff --git a/mhkit/wave/graphics.py b/mhkit/wave/graphics.py
index afb24016a..79d44f3b6 100644
--- a/mhkit/wave/graphics.py
+++ b/mhkit/wave/graphics.py
@@ -26,7 +26,8 @@ def plot_spectrum(S, ax=None):
     ---------
     ax : matplotlib pyplot axes
     """
-    assert isinstance(S, pd.DataFrame), 'S must be of type pd.DataFrame'
+    if not isinstance(S, pd.DataFrame):
+        raise TypeError(f'S must be of type pd.DataFrame. Got: {type(S)}')
 
     f = S.index
     for key in S.keys():
@@ -52,7 +53,8 @@ def plot_elevation_timeseries(eta, ax=None):
     ax : matplotlib pyplot axes
     """
 
-    assert isinstance(eta, pd.DataFrame), 'eta must be of type pd.DataFrame'
+    if not isinstance(eta, pd.DataFrame):
+        raise TypeError(f'eta must be of type pd.DataFrame. Got: {type(eta)}')
 
     for key in eta.keys():
         ax = _xy_plot(eta.index, eta[key], fmt='-', xlabel='Time',
@@ -93,7 +95,8 @@ def plot_matrix(
     ax : matplotlib pyplot axes
     """
 
-    assert isinstance(M, pd.DataFrame), 'M must be of type pd.DataFrame'
+    if not isinstance(M, pd.DataFrame):
+        raise TypeError(f'M must be of type pd.DataFrame. Got: {type(M)}')
 
     if ax is None:
         plt.figure()
@@ -175,22 +178,22 @@ def plot_chakrabarti(H, lambda_w, D, ax=None):
     -------
     ax : matplotlib pyplot axes
     """
-    assert isinstance(H, (np.ndarray, float, int, np.int64,pd.Series)), \
-           'H must be a real numeric type'
-    assert isinstance(lambda_w, (np.ndarray, float, int, np.int64,pd.Series)), \
-           'lambda_w must be a real numeric type'
-    assert isinstance(D, (np.ndarray, float, int, np.int64,pd.Series)), \
-           'D must be a real numeric type'
-
-    if any([(isinstance(H, np.ndarray) or isinstance(H, pd.Series)),        \
-            (isinstance(lambda_w, np.ndarray) or isinstance(H, pd.Series)), \
-            (isinstance(D, np.ndarray) or isinstance(H, pd.Series))\
+    if not isinstance(H, (np.ndarray, float, int, np.int64, pd.Series)):
+        raise TypeError(f'H must be of type float, int, np.int64, np.ndarray, or pd.Series. Got: {type(H)}')
+    if not isinstance(lambda_w, (np.ndarray, float, int, np.int64, pd.Series)):
+        raise TypeError(f'lambda_w must be of type float, int, np.int64, np.ndarray, or pd.Series. Got: {type(lambda_w)}')
+    if not isinstance(D, (np.ndarray, float, int, np.int64, pd.Series)):
+        raise TypeError(f'D must be of type float, int, np.int64, np.ndarray, or pd.Series. Got: {type(D)}')
+
+    if any([isinstance(H, (np.ndarray, pd.Series)),
+            isinstance(lambda_w, (np.ndarray, pd.Series)),
+            isinstance(D, (np.ndarray, pd.Series))
            ]):
-        errMsg = 'D, H, and lambda_w must be same shape'
         n_H = H.squeeze().shape
         n_lambda_w = lambda_w.squeeze().shape
         n_D = D.squeeze().shape
-        assert n_H == n_lambda_w and n_H == n_D, errMsg
+        if not (n_H == n_lambda_w and n_H == n_D):
+            raise ValueError('D, H, and lambda_w must be same shape')
 
         if isinstance(H, np.ndarray):
             mvals = pd.DataFrame(H.reshape(len(H),1), columns=['H'])
@@ -345,55 +348,57 @@ def plot_environmental_contour(x1, x2, x1_contour, x2_contour, **kwargs):
     except: pass
     try: x2 = x2.values
     except: pass
-    assert isinstance(x1, np.ndarray), 'x1 must be of type np.ndarray'
-    assert isinstance(x2, np.ndarray), 'x2 must be of type np.ndarray'
-    assert isinstance(x1_contour, (np.ndarray,list)), ('x1_contour must be of '
-                                                'type np.ndarray or list')
-    assert isinstance(x2_contour, (np.ndarray,list)), ('x2_contour must be of '
-                                               'type np.ndarray or list')
+    if not isinstance(x1, np.ndarray):
+        raise TypeError(f'x1 must be of type np.ndarray. Got: {type(x1)}')
+    if not isinstance(x2, np.ndarray):
+        raise TypeError(f'x2 must be of type np.ndarray. Got: {type(x2)}')
+    if not isinstance(x1_contour, (np.ndarray,list)):
+        raise TypeError(f'x1_contour must be of type np.ndarray or list. Got: {type(x1_contour)}')
+    if not isinstance(x2_contour, (np.ndarray,list)):
+        raise TypeError(f'x2_contour must be of type np.ndarray or list. Got: {type(x2_contour)}')
+    
     x_label = kwargs.get("x_label", None)
     y_label = kwargs.get("y_label", None)
-    data_label=kwargs.get("data_label", None)
-    contour_label=kwargs.get("contour_label", None)
-    ax=kwargs.get("ax", None)
-    markers=kwargs.get("markers", '-')
-    assert isinstance(data_label, (str,type(None))), 'data_label must be of type str'
-    assert isinstance(contour_label, (str,list, type(None))), ('contour_label be of '
-                                                  'type str')
-
-    if isinstance(markers, list):
-        assert all( [isinstance(marker, (str)) for marker in markers] )
-    elif isinstance(markers, str):
-        markers=[markers]
-        assert all( [isinstance(marker, (str)) for marker in markers] )
-    else:
-        assert isinstance(markers, (str,list)), ('markers must be of type str or list of strings')
-
-    assert len(x2_contour) == len(x1_contour),  ('contour must be of'
-            f'equal dimesion got {len(x2_contour)} and {len(x1_contour)}')
-
+    data_label = kwargs.get("data_label", None)
+    contour_label = kwargs.get("contour_label", None)
+    ax = kwargs.get("ax", None)
+    markers = kwargs.get("markers", '-')
+    if not isinstance(data_label, (str, type(None))):
+        raise TypeError(f'If specified, data_label must be of type str. Got: {type(data_label)}')
+    if not isinstance(contour_label, (str, list, type(None))):
+        raise TypeError(f'If specified, contour_label be of type str. Got: {type(contour_label)}')
+
+    if isinstance(markers, str):
+        markers = [markers]
+    if not isinstance(markers, list) or not all( [isinstance(marker, (str)) for marker in markers] ):
+        raise TypeError(f'markers must be of type str or list of strings. Got: {markers}')
+
+    if not len(x2_contour) == len(x1_contour):
+        raise ValueError(f'contour must be of equal dimension got {len(x2_contour)} and {len(x1_contour)}')
 
     if isinstance(x1_contour, np.ndarray):
-        N_contours=1
-        x2_contour  = [x2_contour]
+        N_contours = 1
+        x2_contour = [x2_contour]
         x1_contour = [x1_contour]
     elif isinstance(x1_contour, list):
-        N_contours=len(x1_contour)
+        N_contours = len(x1_contour)
 
     if contour_label != None:
         if isinstance(contour_label, str):
             contour_label = [contour_label]
         N_c_labels = len(contour_label)
-        assert  N_c_labels == N_contours, ('If specified, the '
-             'number of contour lables must be equal to number the '
-            f'number of contour years. Got {N_c_labels} and {N_contours}')
+        if not N_c_labels == N_contours:
+            raise ValueError('If specified, the number of contour labels must'
+                             ' be equal to number the number of contour years.'
+                             f' Got: {N_c_labels} and {N_contours}')
     else:
         contour_label = [None] * N_contours
 
     if len(markers)==1:
-        markers=markers*N_contours
-    assert len(markers) == N_contours,  ('Markers must be same length'
-            f'as N contours specified. Got: {len(markers)} and {len(x1_contour)}')
+        markers = markers*N_contours
+    if not len(markers) == N_contours:
+        raise ValueError('Markers must be same length as N contours specified.'
+                         f'Got: {len(markers)} and {len(x1_contour)}')
 
     for i in range(N_contours):
         contour1 = np.array(x1_contour[i]).T
@@ -543,9 +548,10 @@ def monthly_cumulative_distribution(J):
     ax: axes
         Figure of monthly cumulative distribution
     '''
-    assert isinstance(J, pd.Series), 'J must be of type pd.Series'
-    cumSum={}
-    months=J.index.month.unique()
+    if not isinstance(J, pd.Series):
+        raise TypeError(f'J must be of type pd.Series. Got: {type(J)}')
+    cumSum = {}
+    months = J.index.month.unique()
     for month in months:
         F = exceedance_probability(J[J.index.month==month])
         cumSum[month] = 1-F/100
@@ -592,10 +598,14 @@ def plot_compendium(Hs, Tp, Dp, buoy_title=None, ax=None):
     ax : matplotlib pyplot axes
 
     """
-    assert isinstance(Hs, pd.Series), 'Hs must be of type pd.Series'
-    assert isinstance(Tp, pd.Series), 'Tp must be of type pd.Series'
-    assert isinstance(Dp, pd.Series), 'Dp must be of type pd.Series'
-    assert isinstance(buoy_title, (str, type(None))), 'buoy_title must be of type string'
+    if not isinstance(Hs, pd.Series):
+        raise TypeError(f'Hs must be of type pd.Series. Got: {type(Hs)}')
+    if not isinstance(Tp, pd.Series):
+        raise TypeError(f'Tp must be of type pd.Series. Got: {type(Tp)}')
+    if not isinstance(Dp, pd.Series):
+        raise TypeError(f'Dp must be of type pd.Series. Got: {type(Dp)}')
+    if not isinstance(buoy_title, (str, type(None))):
+        raise TypeError(f'If specified, buoy_title must be of type string. Got: {type(buoy_title)}')
 
     f, (pHs, pTp, pDp) = plt.subplots(3, 1, sharex=True, figsize=(15,10))
 
@@ -668,8 +678,10 @@ def plot_boxplot(Hs, buoy_title=None):
     ---------
     ax : matplotlib pyplot axes
     """
-    assert isinstance(Hs, pd.Series), 'Hs must be of type pd.Series'
-    assert isinstance(buoy_title, (str, type(None))), 'buoy_title must be of type string'
+    if not isinstance(Hs, pd.Series):
+        raise TypeError(f'Hs must be of type pd.Series. Got: {type(Hs)}')
+    if not isinstance(buoy_title, (str, type(None))):
+        raise TypeError(f'If specified, buoy_title must be of type string. Got: {type(buoy_title)}')
 
     months = Hs.index.month
     means = Hs.groupby(months).mean()
@@ -755,7 +767,7 @@ def plot_boxplot(Hs, buoy_title=None):
 
 def plot_directional_spectrum(
         spectrum, 
-        min=None, 
+        color_level_min=None, 
         fill=True, 
         nlevels=11,
         name="Elevation Variance", 
@@ -768,8 +780,8 @@ def plot_directional_spectrum(
     ------------
     spectrum: xarray.DataArray
         Spectral data indexed frequency [Hz] and wave direction [deg].
-    min: float (optional)
-        Minimum value to plot.
+    color_level_min: float (optional)
+        Minimum color bar level.
     fill: bool
         Whether to use `contourf` (filled) instead of `contour` (lines).
     nlevels: int
@@ -783,20 +795,25 @@ def plot_directional_spectrum(
     ---------
     ax : matplotlib pyplot axes
     """
-    assert isinstance(spectrum, xr.DataArray), 'spectrum must be a DataArray'
-    if min is not None:
-        assert isinstance(min, float), 'min must be a float'
-    assert isinstance(fill, bool), 'fill must be a bool'
-    assert isinstance(nlevels, int), 'nlevels must be an int'
-    assert isinstance(name, str), 'name must be a string'
-    assert isinstance(units, str), 'units must be a string'
+    if not isinstance(spectrum, xr.DataArray):
+        raise TypeError(f'spectrum must be of type xr.DataArray. Got: {type(spectrum)}')
+    if not isinstance(color_level_min, (type(None), float)):
+        raise TypeError(f'If specified, color_level_min must be of type float. Got: {type(color_level_min)}')
+    if not isinstance(fill, bool):
+        raise TypeError(f'If specified, fill must be of type bool. Got: {type(fill)}')
+    if not isinstance(nlevels, int):
+        raise TypeError(f'If specified, nlevels must be of type int. Got: {type(nlevels)}')
+    if not isinstance(name, str):
+        raise TypeError(f'If specified, name must be of type string. Got: {type(name)}')
+    if not isinstance(units, str):
+        raise TypeError(f'If specified, units must be of type string. Got: {type(units)}')
 
     a,f = np.meshgrid(np.deg2rad(spectrum.direction), spectrum.frequency)
     _, ax = plt.subplots(subplot_kw=dict(projection='polar'))
     tmp = np.floor(np.min(spectrum.data)*10)/10
-    min = tmp if (min is None) else min
-    max = np.ceil(np.max(spectrum.data)*10)/10
-    levels = np.linspace(min, max, nlevels)
+    color_level_min = tmp if (color_level_min is None) else color_level_min
+    color_level_max = np.ceil(np.max(spectrum.data)*10)/10
+    levels = np.linspace(color_level_min, color_level_max, nlevels)
     if fill:
         c = ax.contourf(a, f, spectrum, levels=levels)
     else:
diff --git a/mhkit/wave/io/hindcast/hindcast.py b/mhkit/wave/io/hindcast/hindcast.py
index f17c14371..61e5b4d20 100644
--- a/mhkit/wave/io/hindcast/hindcast.py
+++ b/mhkit/wave/io/hindcast/hindcast.py
@@ -55,10 +55,10 @@ def region_selection(lat_lon):
         Name of predefined region for given coordinates
     '''
     if not isinstance(lat_lon, (list, tuple)):
-        raise TypeError('lat_lon must be of type list or tuple')
+        raise TypeError(f'lat_lon must be of type list or tuple. Got: {type(lat_lon)}')
 
     if not all(isinstance(coord, (float, int)) for coord in lat_lon):
-        raise TypeError('lat_lon values must be of type float or int')
+        raise TypeError(f'lat_lon values must be of type float or int. Got: {type(lat_lon[0])}')
 
     regions = {
         'Hawaii': {
@@ -85,7 +85,7 @@ def region_search(lat_lon, region, regions):
         lat_lon, region, regions)]
 
     if not region:
-        raise ValueError('ERROR: coordinates out of bounds')
+        raise ValueError('ERROR: coordinates out of bounds.')
 
     return region[0]
 
@@ -165,18 +165,26 @@ def request_wpto_point_data(
     meta: DataFrame 
         Location metadata for the requested data location   
     """
-    assert isinstance(parameter, (str, list)
-                      ), 'parameter must be of type string or list'
-    assert isinstance(lat_lon, (list, tuple)
-                      ), 'lat_lon must be of type list or tuple'
-    assert isinstance(data_type, str), 'data_type must be a string'
-    assert isinstance(years, list), 'years must be a list'
-    assert isinstance(tree, (str, type(None))), 'tree must be a string'
-    assert isinstance(unscale, bool), 'unscale must be bool type'
-    assert isinstance(str_decode, bool), 'str_decode must be bool type'
-    assert isinstance(hsds, bool), 'hsds must be bool type'
-    assert isinstance(path, (str, type(None))), 'path must be a string'
-    assert isinstance(as_xarray, bool), 'as_xarray must be bool type'
+    if not isinstance(parameter, (str, list)):
+        raise TypeError(f'parameter must be of type string or list. Got: {type(parameter)}')
+    if not isinstance(lat_lon, (list, tuple)):
+        raise TypeError(f'lat_lon must be of type list or tuple. Got: {type(lat_lon)}')
+    if not isinstance(data_type, str):
+        raise TypeError(f'data_type must be a string. Got: {type(data_type)}')
+    if not isinstance(years, list):
+        raise TypeError(f'years must be a list. Got: {type(years)}')
+    if not isinstance(tree, (str, type(None))):
+        raise TypeError(f'If specified, tree must be a string. Got: {type(tree)}')
+    if not isinstance(unscale, bool):
+        raise TypeError(f'If specified, unscale must be bool type. Got: {type(unscale)}')
+    if not isinstance(str_decode, bool):
+        raise TypeError(f'If specified, str_decode must be bool type. Got: {type(str_decode)}')
+    if not isinstance(hsds, bool):
+        raise TypeError(f'If specified, hsds must be bool type. Got: {type(hsds)}')
+    if not isinstance(path, (str, type(None))):
+        raise TypeError(f'If specified, path must be a string. Got: {type(path)}')
+    if not isinstance(as_xarray, bool):
+        raise TypeError(f'If specified, as_xarray must be bool type. Got: {type(as_xarray)}')
 
     # Attempt to load data from cache
     # Construct a string representation of the function parameters
@@ -330,13 +338,20 @@ def request_wpto_directional_spectrum(
     meta: DataFrame
         Location metadata for the requested data location
     """
-    assert isinstance(lat_lon, (list, tuple)
-                      ), 'lat_lon must be of type list or tuple'
-    assert isinstance(year, str), 'years must be a string'
-    assert isinstance(tree, (str, type(None))), 'tree must be a sring'
-    assert isinstance(unscale, bool), 'unscale must be bool type'
-    assert isinstance(str_decode, bool), 'str_decode must be bool type'
-    assert isinstance(hsds, bool), 'hsds must be bool type'
+    if not isinstance(lat_lon, (list, tuple)):
+        raise TypeError(f'lat_lon must be of type list or tuple. Got: {type(lat_lon)}')
+    if not isinstance(year, str):
+        raise TypeError(f'year must be a string. Got: {type(year)}')
+    if not isinstance(tree, (str, type(None))):
+        raise TypeError(f'If specified, tree must be a string. Got: {type(tree)}')
+    if not isinstance(unscale, bool):
+        raise TypeError(f'If specified, unscale must be bool type. Got: {type(unscale)}')
+    if not isinstance(str_decode, bool):
+        raise TypeError(f'If specified, str_decode must be bool type. Got: {type(str_decode)}')
+    if not isinstance(hsds, bool):
+        raise TypeError(f'If specified, hsds must be bool type. Got: {type(hsds)}')
+    if not isinstance(path, (str, type(None))):
+        raise TypeError(f'If specified, path must be a string. Got: {type(path)}')
 
     # check for multiple region selection
     if isinstance(lat_lon[0], float):
diff --git a/mhkit/wave/io/ndbc.py b/mhkit/wave/io/ndbc.py
index 752ea0910..e5613553d 100644
--- a/mhkit/wave/io/ndbc.py
+++ b/mhkit/wave/io/ndbc.py
@@ -47,9 +47,10 @@ def read_file(file_name, missing_values=['MM', 9999, 999, 99]):
         Dictionary with {column name: units} key value pairs when the NDBC file
         contains unit information, otherwise None is returned
     """
-    assert isinstance(file_name, str), 'file_name must be of type str'
-    assert isinstance(
-        missing_values, list), 'missing_values must be of type list'
+    if not isinstance(file_name, str):
+        raise TypeError(f'file_name must be of type str. Got: {type(file_name)}')
+    if not isinstance(missing_values, list):
+        raise TypeError(f'If specified, missing_values must be of type list. Got: {type(missing_values)}')
 
     # Open file and get header rows
     f = open(file_name, "r")
@@ -149,19 +150,22 @@ def available_data(parameter, buoy_number=None, proxy=None, clear_cache=False):
     available_data: DataFrame
         DataFrame with station ID, years, and NDBC file names.
     '''
-    assert isinstance(parameter, str), 'parameter must be a string'
-    assert isinstance(buoy_number, (str, type(None), list)), ('If '
-                                                              'specified the buoy number must be a string or list of strings')
-    assert isinstance(proxy, (dict, type(None))
-                      ), 'If specified proxy must be a dict'
+    if not isinstance(parameter, str):
+        raise TypeError(f'parameter must be a string. Got: {type(parameter)}')
+    if not isinstance(buoy_number, (str, type(None), list)):
+        raise TypeError(f'If specified, buoy_number must be a string or list of strings. Got: {type(buoy_number)}')
+    if not isinstance(proxy, (dict, type(None))):
+        raise TypeError(f'If specified, proxy must be a dict. Got: {type(proxy)}')
     _supported_params(parameter)
     if isinstance(buoy_number, str):
-        assert len(buoy_number) == 5, ('Buoy must be 5-character'
-                                       f'alpha-numeric station identifier got: {buoy_number}')
+        if not len(buoy_number) == 5:
+            raise ValueError('buoy_number must be 5-character'
+                             f'alpha-numeric station identifier. Got: {buoy_number}')
     elif isinstance(buoy_number, list):
         for buoy in buoy_number:
-            assert len(buoy) == 5, ('Each buoy must be a 5-character'
-                                    f'alpha-numeric station identifier got: {buoy}')
+            if not len(buoy) == 5:
+                raise ValueError('Each value in the buoy_number list must be a 5-character'
+                                 f'alpha-numeric station identifier. Got: {buoy_number}')
 
     # Generate a unique hash_params based on the function parameters
     hash_params = f"parameter:{parameter}_buoy_number:{buoy_number}_proxy:{proxy}"
@@ -240,9 +244,10 @@ def _parse_filenames(parameter, filenames):
     buoys: DataFrame
         DataFrame with keys=['id','year','file_name']
     '''
-    assert isinstance(
-        filenames, pd.Series), 'filenames must be of type pd.Series'
-    assert isinstance(parameter, str), 'parameter must be a string'
+    if not isinstance(filenames, pd.Series):
+        raise TypeError(f'filenames must be of type pd.Series. Got: {type(filenames)}')
+    if not isinstance(parameter, str):
+        raise TypeError(f'parameter must be a string. Got: {type(parameter)}')
     supported = _supported_params(parameter)
 
     file_seps = {
@@ -299,16 +304,18 @@ def request_data(parameter, filenames, proxy=None, clear_cache=False):
     ndbc_data: dict
         Dictionary of DataFrames indexed by buoy and year.
     '''
-    assert isinstance(filenames, (pd.Series, pd.DataFrame)), (
-        'filenames must be of type pd.Series')
-    assert isinstance(parameter, str), 'parameter must be a string'
-    assert isinstance(proxy, (dict, type(None))), (
-        'If specified proxy must be a dict')
+    if not isinstance(filenames, (pd.Series, pd.DataFrame)):
+        raise TypeError(f'filenames must be of type pd.Series or pd.DataFrame. Got: {type(filenames)}')
+    if not isinstance(parameter, str):
+        raise TypeError(f'parameter must be a string. Got: {type(parameter)}')
+    if not isinstance(proxy, (dict, type(None))):
+        raise TypeError(f'If specified, proxy must be a dict. Got: {type(proxy)}')
 
     _supported_params(parameter)
     if isinstance(filenames, pd.DataFrame):
         filenames = pd.Series(filenames.squeeze())
-    assert len(filenames) > 0, "At least 1 filename must be passed"
+    if not len(filenames) > 0:
+        raise ValueError("At least 1 filename must be passed")
 
     # Define the path to the cache directory
     cache_dir = os.path.join(os.path.expanduser("~"),
@@ -395,9 +402,10 @@ def to_datetime_index(parameter, ndbc_data):
             Dataframe with NDBC date columns removed, and datetime index
     '''
 
-    assert isinstance(parameter, str), 'parameter must be a string'
-    assert isinstance(
-        ndbc_data, pd.DataFrame), 'ndbc_data must be of type pd.DataFrame'
+    if not isinstance(parameter, str):
+        raise TypeError(f'parameter must be a string. Got: {type(parameter)}')
+    if not isinstance(ndbc_data, pd.DataFrame):
+        raise TypeError(f'ndbc_data must be of type pd.DataFrame. Got: {type(ndbc_data)}')
 
     df_datetime = ndbc_data.copy(deep=True)
     df_datetime['date'], ndbc_date_cols = dates_to_datetime(
@@ -437,9 +445,10 @@ def dates_to_datetime(data, return_date_cols=False, return_as_dataframe=False):
         List of the DataFrame columns headers for dates as provided by
         NDBC
     '''
-    assert isinstance(data, pd.DataFrame), 'data must be of type pd.DataFrame'
-    assert isinstance(return_date_cols,
-                      bool), 'return_date_cols must be of type bool'
+    if not isinstance(data, pd.DataFrame):
+        raise TypeError(f'data must be of type pd.DataFrame. Got: {type(data)}')
+    if not isinstance(return_date_cols,bool):
+        raise TypeError(f'return_date_cols must be of type bool. Got: {type(return_date_cols)}')
 
     df = data.copy(deep=True)
     cols = df.columns.values.tolist()
@@ -514,9 +523,12 @@ def _date_string_to_datetime(df, columns, year_fmt):
     df: DataFrame
         The passed df with a new column ['date'] with the datetime format
     '''
-    assert isinstance(df, pd.DataFrame), 'df must be of type pd.DataFrame'
-    assert isinstance(columns, list), 'Columns must be a list'
-    assert isinstance(year_fmt, str), 'year_fmt must be a string'
+    if not isinstance(df, pd.DataFrame):
+        raise TypeError(f'df must be of type pd.DataFrame. Got: {type(df)}')
+    if not isinstance(columns, list):
+       raise TypeError(f'columns must be a list. Got: {type(columns)}')
+    if not isinstance(year_fmt, str):
+        raise TypeError(f'year_fmt must be a string. Got: {type(year_fmt)}')
 
     # Convert to str and zero pad
     for key in columns:
@@ -572,7 +584,8 @@ def parameter_units(parameter=''):
         Dictionary of parameter units
     '''
 
-    assert isinstance(parameter, str), 'parameter must be a string'
+    if not isinstance(parameter, str):
+        raise TypeError(f'parameter must be a string. Got: {type(parameter)}')
 
     if parameter == 'adcp':
         units = {'DEP01': 'm',
@@ -750,7 +763,8 @@ def _supported_params(parameter):
     msg: bool
         Whether the parameter is supported.
     '''
-    assert isinstance(parameter, str), 'parameter must be a string'
+    if not isinstance(parameter, str):
+        raise TypeError(f'parameter must be a string. Got: {type(parameter)}')
     supported = True
     supported_params = [
         'swden',
@@ -838,8 +852,10 @@ def request_directional_data(buoy, year):
         Dataset containing the five parameter data indexed by frequency
         and date.
     """
-    assert isinstance(buoy, str), 'buoy must be a string'
-    assert isinstance(year, int), 'year must be an int'
+    if not isinstance(buoy, str):
+        raise TypeError(f'buoy must be a string. Got: {type(buoy)}')
+    if not isinstance(year, int):
+        raise TypeError(f'year must be an int. Got: {type(year)}')
 
     directional_parameters = ['swden', 'swdir', 'swdir2', 'swr1', 'swr2']
 
@@ -936,15 +952,21 @@ def _create_spectrum(data, frequencies, directions, name, units):
         DataArray containing the spectrum values indexed by frequency
         and wave direction.
     """
-    assert isinstance(data, np.ndarray), 'data must be an array'
-    assert isinstance(frequencies, np.ndarray), 'frequencies must be an array'
-    assert isinstance(directions, np.ndarray), 'directions must be an array'
-    assert isinstance(name, str), 'name must be a string'
-    assert isinstance(units, str), 'units must be a string'
+    if not isinstance(data, np.ndarray):
+        raise TypeError(f'data must be of type np.ndarray. Got: {type(data)}')
+    if not isinstance(frequencies, np.ndarray):
+        raise TypeError(f'frequencies must be of type np.ndarray. Got: {type(frequencies)}')
+    if not isinstance(directions, np.ndarray):
+        raise TypeError(f'directions must be of type np.ndarray. Got: {type(directions)}')
+    if not isinstance(name, str):
+        raise TypeError(f'name must be of type string. Got: {type(name)}')
+    if not isinstance(units, str):
+        raise TypeError(f'units must be of type string. Got: {type(units)}')
 
     msg = (f'data has wrong shape {data.shape}, ' +
            f'expected {(len(frequencies), len(directions))}')
-    assert data.shape == (len(frequencies), len(directions)), msg
+    if not data.shape == (len(frequencies), len(directions)):
+        raise ValueError(msg)
 
     direction_attrs = {
         'units': 'deg',
@@ -994,8 +1016,10 @@ def create_spread_function(data, directions):
         DataArray containing the spread function values indexed by
         frequency and wave direction.
     """
-    assert isinstance(data, xr.Dataset), 'data must be a Dataset'
-    assert isinstance(directions, np.ndarray), 'directions must be an array'
+    if not isinstance(data, xr.Dataset):
+        raise TypeError(f'data must be of type xr.Dataset. Got: {type(data)}')
+    if not isinstance(directions, np.ndarray):
+        raise TypeError(f'directions must be of type np.ndarray. Got: {type(directions)}')
 
     r1 = data['swr1'].data.reshape(-1, 1)
     r2 = data['swr2'].data.reshape(-1, 1)
@@ -1037,8 +1061,10 @@ def create_directional_spectrum(data, directions):
         DataArray containing the spectrum values indexed by frequency
         and wave direction.
     """
-    assert isinstance(data, xr.Dataset), 'data must be a Dataset'
-    assert isinstance(directions, np.ndarray), 'directions must be an array'
+    if not isinstance(data, xr.Dataset):
+        raise TypeError(f'data must be of type xr.Dataset. Got: {type(data)}')
+    if not isinstance(directions, np.ndarray):
+        raise TypeError(f'directions must be of type np.ndarray. Got: {type(directions)}')
 
     spread = create_spread_function(data, directions).values
     omnidirectional_spectrum = data['swden'].data.reshape(-1, 1)
diff --git a/mhkit/wave/io/swan.py b/mhkit/wave/io/swan.py
index c71a1a514..27aae18c2 100644
--- a/mhkit/wave/io/swan.py
+++ b/mhkit/wave/io/swan.py
@@ -21,8 +21,10 @@ def read_table(swan_file):
     metaDict: Dictionary
         Dictionary of metaData
     '''
-    assert isinstance(swan_file, str), 'swan_file must be of type str'
-    assert isfile(swan_file)==True, f'File not found: {swan_file}'
+    if not isinstance(swan_file, str):
+        raise TypeError(f'swan_file must be of type str. Got: {type(swan_file)}')
+    if not isfile(swan_file):
+        raise ValueError(f'File not found: {swan_file}')
     
     f = open(swan_file,'r')
     header_line_number = 4
@@ -62,8 +64,10 @@ def read_block(swan_file):
     metaDict: Dictionary
         Dictionary of metaData dependent on file type    
     '''
-    assert isinstance(swan_file, str), 'swan_file must be of type str'
-    assert isfile(swan_file)==True, f'File not found: {swan_file}'
+    if not isinstance(swan_file, str):
+        raise TypeError(f'swan_file must be of type str. Got: {type(swan_file)}')
+    if not isfile(swan_file):
+        raise ValueError(f'File not found: {swan_file}')
     
     extension = swan_file.split('.')[1].lower()
     if extension == 'mat':
@@ -92,8 +96,10 @@ def _read_block_txt(swan_file):
     metaDict: Dictionary
         Dictionary of metaData dependent on file type    
     '''
-    assert isinstance(swan_file, str), 'swan_file must be of type str'
-    assert isfile(swan_file)==True, f'File not found: {swan_file}'
+    if not isinstance(swan_file, str):
+        raise TypeError(f'swan_file must be of type str. Got: {type(swan_file)}')
+    if not isfile(swan_file):
+        raise ValueError(f'File not found: {swan_file}')
     
     f = open(swan_file) 
     runLines=[]
@@ -174,8 +180,10 @@ def _read_block_mat(swan_file):
     dataDict: Dictionary
         Dictionary of DataFrame of swan output variables
     '''
-    assert isinstance(swan_file, str), 'swan_file must be of type str'
-    assert isfile(swan_file)==True, f'File not found: {swan_file}'
+    if not isinstance(swan_file, str):
+        raise TypeError(f'swan_file must be of type str. Got: {type(swan_file)}')
+    if not isfile(swan_file):
+        raise ValueError(f'File not found: {swan_file}')
         
     dataDict = loadmat(swan_file, struct_as_record=False, squeeze_me=True)
     removeKeys = ['__header__', '__version__', '__globals__']
@@ -200,7 +208,8 @@ def _parse_line_metadata(line):
     metaDict: Dictionary
         Dictionary of variable metadata
     '''
-    assert isinstance(line, str), 'line must be of type str'
+    if not isinstance(line, str):
+        raise TypeError(f'line must be of type str. Got: {type(line)}')
     
     metaDict={}        
     meta=re.sub('\s+', " ", line.replace(',', ' ').strip('% \n').replace('**', 'vars:'))
@@ -239,19 +248,20 @@ def dictionary_of_block_to_table(dictionary_of_DataFrames, names=None):
         DataFrame with columns x,y,values where values = Dictionary.keys()
         or names        
     '''
-    assert isinstance(dictionary_of_DataFrames, dict), (
-                          'dictionary_of_DataFrames must be of type Dict')
-    assert bool(dictionary_of_DataFrames), 'dictionary_of_DataFrames is empty'
+    if not isinstance(dictionary_of_DataFrames, dict):
+        raise TypeError(f'dictionary_of_DataFrames must be of type dict. Got: {type(dictionary_of_DataFrames)}')
+    if not bool(dictionary_of_DataFrames):
+        raise ValueError(f'dictionary_of_DataFrames is empty. Got: {dictionary_of_DataFrames}')
     for key in dictionary_of_DataFrames:  
-        assert isinstance(dictionary_of_DataFrames[key],pd.DataFrame), (
-                          f'Dictionary key:{key} must be of type pd.DataFrame')
+        if not isinstance(dictionary_of_DataFrames[key],pd.DataFrame):
+            raise TypeError(f'Dictionary key:{key} must be of type pd.DataFrame. Got: {type(dictionary_of_DataFrames[key])}')
     if not isinstance(names, type(None)):
-        assert isinstance(names, list), (
-                'If specified names must be of type list')         
-        assert all([isinstance(elm, str) for elm in names]), (
-                'If specified all elements in names must be of type string')
-        assert len(names) == len(dictionary_of_DataFrames), (
-                'If specified names must the same length as dictionary_of_DataFrames')
+        if not isinstance(names, list):
+            raise TypeError(f'If specified, names must be of type list. Got: {type(names)}')
+        if not all([isinstance(elm, str) for elm in names]):
+                raise ValueError(f'If specified, all elements in names must be of type string. Got: {names}')
+        if not len(names) == len(dictionary_of_DataFrames):
+            raise ValueError('If specified, names must the same length as dictionary_of_DataFrames')
     
     if names == None:
         variables = [var for var in  dictionary_of_DataFrames.keys() ]
@@ -284,8 +294,10 @@ def block_to_table(data, name='values'):
     table: DataFrame
         DataFrame with columns x,y,values           
     '''
-    assert isinstance(data,pd.DataFrame), 'data must be of type pd.DataFrame'
-    assert isinstance(name, str), 'Name must be of type str'
+    if not isinstance(data,pd.DataFrame):
+        raise TypeError(f'data must be of type pd.DataFrame. Got: {type(data)}')
+    if not isinstance(name, str):
+        raise TypeError(f'If specified, name must be of type str. Got: {type(name)}')
     
     table = data.unstack().reset_index(name=name)
     table = table.rename(columns={'level_0':'x', 'level_1': 'y'})
diff --git a/mhkit/wave/performance.py b/mhkit/wave/performance.py
index 2b96809a9..80404fd87 100644
--- a/mhkit/wave/performance.py
+++ b/mhkit/wave/performance.py
@@ -23,8 +23,10 @@ def capture_length(P, J):
     L: numpy array or pandas Series
         Capture length [m]
     """
-    assert isinstance(P, (np.ndarray, pd.Series)), 'P must be of type np.ndarray or pd.Series'
-    assert isinstance(J, (np.ndarray, pd.Series)), 'J must be of type np.ndarray or pd.Series'
+    if not isinstance(P, (np.ndarray, pd.Series)):
+        raise TypeError(f'P must be of type np.ndarray or pd.Series. Got: {type(P)}')
+    if not isinstance(J, (np.ndarray, pd.Series)):
+        raise TypeError(f'J must be of type np.ndarray or pd.Series. Got: {type(J)}')
 
     L = P/J
 
@@ -49,7 +51,8 @@ def statistics(X):
     stats: pandas Series
         Statistics
     """
-    assert isinstance(X, (np.ndarray, pd.Series)), 'X must be of type np.ndarray or pd.Series'
+    if not isinstance(X, (np.ndarray, pd.Series)):
+        raise TypeError(f'X must be of type np.ndarray or pd.Series. Got: {type(X)}')
 
     stats = pd.Series(X).describe()
     stats['std'] = _std_ddof1(X)
@@ -128,12 +131,18 @@ def capture_length_matrix(Hm0, Te, L, statistic, Hm0_bins, Te_bins):
         equal to Te_bins
 
     """
-    assert isinstance(Hm0, (np.ndarray, pd.Series)), 'Hm0 must be of type np.ndarray or pd.Series'
-    assert isinstance(Te, (np.ndarray, pd.Series)), 'Te must be of type np.ndarray or pd.Series'
-    assert isinstance(L, (np.ndarray, pd.Series)), 'L must be of type np.ndarray or pd.Series'
-    assert isinstance(statistic, (str, types.FunctionType)), 'statistic must be of type str or callable'
-    assert isinstance(Hm0_bins, np.ndarray), 'Hm0_bins must be of type np.ndarray'
-    assert isinstance(Te_bins, np.ndarray), 'Te_bins must be of type np.ndarray'
+    if not isinstance(Hm0, (np.ndarray, pd.Series)):
+        raise TypeError(f'Hm0 must be of type np.ndarray or pd.Series. Got: {type(Hm0)}')
+    if not isinstance(Te, (np.ndarray, pd.Series)):
+        raise TypeError(f'Te must be of type np.ndarray or pd.Series. Got: {type(Te)}')
+    if not isinstance(L, (np.ndarray, pd.Series)):
+        raise TypeError(f'L must be of type np.ndarray or pd.Series. Got: {type(L)}')
+    if not isinstance(statistic, (str, types.FunctionType)):
+        raise TypeError(f'statistic must be of type str or callable. Got: {type(statistic)}')
+    if not isinstance(Hm0_bins, np.ndarray):
+        raise TypeError(f'Hm0_bins must be of type np.ndarray. Got: {type(Hm0_bins)}')
+    if not isinstance(Te_bins, np.ndarray):
+        raise TypeError(f'Te_bins must be of type np.ndarray. Got: {type(Te_bins)}')
 
     LM = _performance_matrix(Hm0, Te, L, statistic, Hm0_bins, Te_bins)
 
@@ -168,12 +177,18 @@ def wave_energy_flux_matrix(Hm0, Te, J, statistic, Hm0_bins, Te_bins):
         equal to Te_bins
 
     """
-    assert isinstance(Hm0, (np.ndarray, pd.Series)), 'Hm0 must be of type np.ndarray or pd.Series'
-    assert isinstance(Te, (np.ndarray, pd.Series)), 'Te must be of type np.ndarray or pd.Series'
-    assert isinstance(J, (np.ndarray, pd.Series)), 'J must be of type np.ndarray or pd.Series'
-    assert isinstance(statistic, (str, callable)), 'statistic must be of type str or callable'
-    assert isinstance(Hm0_bins, np.ndarray), 'Hm0_bins must be of type np.ndarray'
-    assert isinstance(Te_bins, np.ndarray), 'Te_bins must be of type np.ndarray'
+    if not isinstance(Hm0, (np.ndarray, pd.Series)):
+        raise TypeError(f'Hm0 must be of type np.ndarray or pd.Series. Got: {type(Hm0)}')
+    if not isinstance(Te, (np.ndarray, pd.Series)):
+        raise TypeError(f'Te must be of type np.ndarray or pd.Series. Got: {type(Te)}')
+    if not isinstance(J, (np.ndarray, pd.Series)):
+        raise TypeError(f'J must be of type np.ndarray or pd.Series. Got: {type(J)}')
+    if not isinstance(statistic, (str, callable)):
+        raise TypeError(f'statistic must be of type str or callable. Got: {type(statistic)}')
+    if not isinstance(Hm0_bins, np.ndarray):
+        raise TypeError(f'Hm0_bins must be of type np.ndarray. Got: {type(Hm0_bins)}')
+    if not isinstance(Te_bins, np.ndarray):
+        raise TypeError(f'Te_bins must be of type np.ndarray. Got: {type(Te_bins)}')
 
     JM = _performance_matrix(Hm0, Te, J, statistic, Hm0_bins, Te_bins)
 
@@ -197,8 +212,10 @@ def power_matrix(LM, JM):
         Power matrix
 
     """
-    assert isinstance(LM, pd.DataFrame), 'LM must be of type pd.DataFrame'
-    assert isinstance(JM, pd.DataFrame), 'JM must be of type pd.DataFrame'
+    if not isinstance(LM, pd.DataFrame):
+        raise TypeError(f'LM must be of type pd.DataFrame. Got: {type(LM)}')
+    if not isinstance(JM, pd.DataFrame):
+        raise TypeError(f'JM must be of type pd.DataFrame. Got: {type(JM)}')
 
     PM = LM*JM
 
@@ -221,8 +238,10 @@ def mean_annual_energy_production_timeseries(L, J):
         Mean annual energy production
 
     """
-    assert isinstance(L, (np.ndarray, pd.Series)), 'L must be of type np.ndarray or pd.Series'
-    assert isinstance(J, (np.ndarray, pd.Series)), 'J must be of type np.ndarray or pd.Series'
+    if not isinstance(L, (np.ndarray, pd.Series)):
+        raise TypeError(f'L must be of type np.ndarray or pd.Series. Got: {type(L)}')
+    if not isinstance(J, (np.ndarray, pd.Series)):
+        raise TypeError(f'J must be of type np.ndarray or pd.Series. Got: {type(J)}')
 
     T = 8766 # Average length of a year (h)
     n = len(L)
@@ -251,11 +270,15 @@ def mean_annual_energy_production_matrix(LM, JM, frequency):
         Mean annual energy production
 
     """
-    assert isinstance(LM, pd.DataFrame), 'LM must be of type pd.DataFrame'
-    assert isinstance(JM, pd.DataFrame), 'JM must be of type pd.DataFrame'
-    assert isinstance(frequency, pd.DataFrame), 'frequency must be of type pd.DataFrame'
-    assert LM.shape == JM.shape == frequency.shape, 'LM, JM, and frequency must be of the same size'
-    #assert frequency.sum().sum() == 1
+    if not isinstance(LM, pd.DataFrame):
+        raise TypeError(f'LM must be of type pd.DataFrame. Got: {type(LM)}')
+    if not isinstance(JM, pd.DataFrame):
+        raise TypeError(f'JM must be of type pd.DataFrame. Got: {type(JM)}')
+    if not isinstance(frequency, pd.DataFrame):
+        raise TypeError(f'frequency must be of type pd.DataFrame. Got: {type(frequency)}')
+    if not LM.shape == JM.shape == frequency.shape:
+        raise ValueError('LM, JM, and frequency must be of the same size')
+    #if not frequency.sum().sum() == 1
 
     T = 8766 # Average length of a year (h)
     maep = T * np.nansum(LM * JM * frequency)
@@ -309,13 +332,20 @@ def power_performance_workflow(S, h, P, statistic, frequency_bins=None, deep=Fal
     maep_matrix: float
         Mean annual energy production
     """
-    assert isinstance(S, (pd.DataFrame,pd.Series)), 'S must be of type pd.DataFrame or pd.Series'
-    assert isinstance(h, (int,float)), 'h must be of type int or float'
-    assert isinstance(P, (np.ndarray, pd.Series)), 'P must be of type np.ndarray or pd.Series'
-    assert isinstance(deep, bool), 'deep must be of type bool'
-    assert isinstance(rho, (int,float)), 'rho must be of type int or float'
-    assert isinstance(g, (int,float)), 'g must be of type int or float'
-    assert isinstance(ratio, (int,float)), 'ratio must be of type int or float'
+    if not isinstance(S, (pd.DataFrame,pd.Series)):
+        raise TypeError(f'S must be of type pd.DataFrame or pd.Series. Got: {type(S)}')
+    if not isinstance(h, (int,float)):
+        raise TypeError(f'h must be of type int or float. Got: {type(h)}')
+    if not isinstance(P, (np.ndarray, pd.Series)):
+        raise TypeError(f'P must be of type np.ndarray or pd.Series. Got: {type(P)}')
+    if not isinstance(deep, bool):
+        raise TypeError(f'deep must be of type bool. Got: {type(deep)}')
+    if not isinstance(rho, (int,float)):
+        raise TypeError(f'rho must be of type int or float. Got: {type(rho)}')
+    if not isinstance(g, (int,float)):
+        raise TypeError(f'g must be of type int or float. Got: {type(g)}')
+    if not isinstance(ratio, (int,float)):
+        raise TypeError(f'ratio must be of type int or float. Got: {type(ratio)}')
 
     # Compute the enegy periods from the spectra data
     Te = wave.resource.energy_period(S, frequency_bins=frequency_bins)
diff --git a/mhkit/wave/resource.py b/mhkit/wave/resource.py
index 5e6f54790..f0fd27b0f 100644
--- a/mhkit/wave/resource.py
+++ b/mhkit/wave/resource.py
@@ -36,15 +36,22 @@ def elevation_spectrum(eta, sample_rate, nnft, window='hann',
 
     # TODO: Add confidence intervals, equal energy frequency spacing, and NDBC
     #       frequency spacing
-    # TODO: may need an assert for the length of nnft- signal.welch breaks when nfft is too short
+    # TODO: may need to raise an error for the length of nnft- signal.welch breaks when nfft is too short
     # TODO: check for uniform sampling
-    assert isinstance(eta, pd.DataFrame), 'eta must be of type pd.DataFrame'
-    assert isinstance(sample_rate, (float,int)), 'sample_rate must be of type int or float'
-    assert isinstance(nnft, int), 'nnft must be of type int'
-    assert isinstance(window, str), 'window must be of type str'
-    assert isinstance(detrend, bool), 'detrend must be of type bool'
-    assert nnft > 0, 'nnft must be > 0'
-    assert sample_rate > 0, 'sample_rate must be > 0'
+    if not isinstance(eta, pd.DataFrame):
+        raise TypeError(f'eta must be of type pd.DataFrame. Got: {type(eta)}')
+    if not isinstance(sample_rate, (float,int)):
+        raise TypeError(f'sample_rate must be of type int or float. Got: {type(sample_rate)}')
+    if not isinstance(nnft, int):
+        raise TypeError(f'nnft must be of type int. Got: {type(nnft)}')
+    if not isinstance(window, str):
+        raise TypeError(f'window must be of type str. Got: {type(window)}')
+    if not isinstance(detrend, bool):
+        raise TypeError(f'detrend must be of type bool. Got: {type(detrend)}')
+    if not nnft > 0:
+        raise ValueError(f'nnft must be > 0. Got: {nnft}')
+    if not sample_rate > 0:
+        raise ValueError(f'sample_rate must be > 0. Got: {sample_rate}')
 
     S = pd.DataFrame()
     for col in eta.columns:
@@ -83,9 +90,12 @@ def pierson_moskowitz_spectrum(f, Tp, Hs):
         f = np.array(f)
     except:
         pass
-    assert isinstance(f, np.ndarray), 'f must be of type np.ndarray'
-    assert isinstance(Tp, (int,float)), 'Tp must be of type int or float'
-    assert isinstance(Hs, (int,float)), 'Hs must be of type int or float'
+    if not isinstance(f, np.ndarray):
+        raise TypeError(f'f must be of type np.ndarray. Got: {type(f)}')
+    if not isinstance(Tp, (int,float)):
+        raise TypeError(f'Tp must be of type int or float. Got: {type(Tp)}')
+    if not isinstance(Hs, (int,float)):
+        raise TypeError(f'Hs must be of type int or float. Got: {type(Hs)}')
 
     f.sort()
     B_PM = (5/4)*(1/Tp)**4
@@ -133,11 +143,14 @@ def jonswap_spectrum(f, Tp, Hs, gamma=None):
         f = np.array(f)
     except:
         pass
-    assert isinstance(f, np.ndarray), 'f must be of type np.ndarray'
-    assert isinstance(Tp, (int,float)), 'Tp must be of type int or float'
-    assert isinstance(Hs, (int,float)), 'Hs must be of type int or float'
-    assert isinstance(gamma, (int,float, type(None))), \
-        'gamma must be of type int or float'
+    if not isinstance(f, np.ndarray):
+        raise TypeError(f'f must be of type np.ndarray. Got: {type(f)}')
+    if not isinstance(Tp, (int,float)):
+        raise TypeError(f'Tp must be of type int or float. Got: {type(Tp)}')
+    if not isinstance(Hs, (int,float)):
+        raise TypeError(f'Hs must be of type int or float. Got: {type(Hs)}')
+    if not isinstance(gamma, (int,float, type(None))):
+        raise TypeError(f'If specified, gamma must be of type int or float. Got: {type(gamma)}')
 
     f.sort()
     B_PM = (5/4)*(1/Tp)**4
@@ -215,41 +228,48 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
 
     """
     time_index = np.array(time_index)
-    assert isinstance(S, pd.DataFrame), 'S must be of type pd.DataFrame'
-    assert isinstance(time_index, np.ndarray), ('time_index must be of type'
-            'np.ndarray')
-    assert isinstance(seed, (type(None),int)), 'seed must be of type int'
-    assert isinstance(frequency_bins, (type(None), np.ndarray, pd.DataFrame)),(
-            "frequency_bins must be of type None, np.ndarray, or pd,DataFrame")
-    assert isinstance(phases, (type(None), np.ndarray, pd.DataFrame)), (
-            'phases must be of type None, np.ndarray, or pd,DataFrame')
-    assert isinstance(method, str)
+    if not isinstance(S, pd.DataFrame):
+        raise TypeError(f'S must be of type pd.DataFrame. Got: {type(S)}')
+    if not isinstance(time_index, np.ndarray):
+        raise TypeError(f'time_index must be of type np.ndarray. Got: {type(time_index)}')
+    if not isinstance(seed, (type(None), int)):
+        raise TypeError(f'If specified, seed must be of type int. Got: {type(seed)}')
+    if not isinstance(frequency_bins, (type(None), np.ndarray, pd.DataFrame)):
+        raise TypeError(f'If specified, frequency_bins must be of type np.ndarray, or pd.DataFrame. Got: {type(frequency_bins)}')
+    if not isinstance(phases, (type(None), np.ndarray, pd.DataFrame)):
+        raise TypeError(f'If specified, phases must be of type np.ndarray, or pd.DataFrame. Got: {type(phases)}')
+    if not isinstance(method, str):
+        raise TypeError(f'method must be of type str. Got: {type(method)}')
 
     if frequency_bins is not None:
-        assert frequency_bins.squeeze().shape == (S.squeeze().shape[0],),(
-            'shape of frequency_bins must match shape of S')
+        if not frequency_bins.squeeze().shape == (S.squeeze().shape[0],):
+            raise ValueError('shape of frequency_bins must match shape of S')
     if phases is not None:
-        assert phases.squeeze().shape == S.squeeze().shape,(
-            'shape of phases must match shape of S')
+        if not phases.squeeze().shape == S.squeeze().shape:
+            raise ValueError('shape of phases must match shape of S')
         
     if method is not None:
-        assert method == 'ifft' or method == 'sum_of_sines',(
-            f"unknown method {method}, options are 'ifft' or 'sum_of_sines'")
+        if not (method == 'ifft' or method == 'sum_of_sines'):
+            raise ValueError(f"Method must be 'ifft' or 'sum_of_sines'. Got: {method}")
         
     if method == 'ifft':
-        assert S.index.values[0] == 0, ('ifft method must have zero frequency defined')        
+        if not S.index.values[0] == 0:
+            raise ValueError(f'ifft method must have zero frequency defined. Lowest frequency is: {S.index.values[0]}')        
 
     f = pd.Series(S.index)
     f.index = f
     if frequency_bins is None:
         delta_f = f.values[1]-f.values[0]
-        assert np.allclose(f.diff()[1:], delta_f)
+        if not np.allclose(f.diff()[1:], delta_f):
+            raise ValueError('Frequency bins are not evenly spaced. ' +
+                             "Define 'frequency_bins' or create a constant " +
+                             'frequency spacing for S.')
     elif isinstance(frequency_bins, np.ndarray):
         delta_f = pd.Series(frequency_bins, index=S.index)
         method = 'sum_of_sines'
     elif isinstance(frequency_bins, pd.DataFrame):
-        assert len(frequency_bins.columns) == 1, ('frequency_bins must only'
-                'contain 1 column')
+        if not len(frequency_bins.columns) == 1:
+            raise ValueError('frequency_bins must only contain 1 column')
         delta_f = frequency_bins.squeeze()
         method = 'sum_of_sines'
 
@@ -313,8 +333,10 @@ def frequency_moment(S, N, frequency_bins=None):
     m: pandas DataFrame
         Nth Frequency Moment indexed by S.columns
     """
-    assert isinstance(S, (pd.Series,pd.DataFrame)), 'S must be of type pd.DataFrame or pd.Series'
-    assert isinstance(N, int), 'N must be of type int'
+    if not isinstance(S, (pd.Series,pd.DataFrame)):
+        raise TypeError(f'S must be of type pd.DataFrame or pd.Series. Got: {type(S)}')
+    if not isinstance(N, int):
+        raise TypeError(f'N must be of type int. Got: {type(N)}')
 
     # Eq 8 in IEC 62600-101
     spec = S[S.index > 0] # omit frequency of 0
@@ -326,8 +348,8 @@ def frequency_moment(S, N, frequency_bins=None):
         delta_f[0] = f[1]-f[0]
     else:
 
-        assert isinstance(frequency_bins, (np.ndarray,pd.Series,pd.DataFrame)),(
-         'frequency_bins must be of type np.ndarray or pd.Series')
+        if not isinstance(frequency_bins, (np.ndarray,pd.Series,pd.DataFrame)):
+            raise TypeError(f'frequency_bins must be of type np.ndarray, pd.Series, or pd.DataFrame. Got: {type(frequency_bins)}')
         delta_f = pd.Series(frequency_bins)
 
     delta_f.index = f
@@ -358,7 +380,8 @@ def significant_wave_height(S, frequency_bins=None):
     Hm0: pandas DataFrame
         Significant wave height [m] index by S.columns
     """
-    assert isinstance(S, (pd.Series,pd.DataFrame)), 'S must be of type pd.DataFrame or pd.Series'
+    if not isinstance(S, (pd.Series,pd.DataFrame)):
+        raise TypeError(f'S must be of type pd.DataFrame or pd.Series. Got: {type(S)}')
 
     # Eq 12 in IEC 62600-101
 
@@ -384,7 +407,8 @@ def average_zero_crossing_period(S,frequency_bins=None):
     Tz: pandas DataFrame
         Average zero crossing period [s] indexed by S.columns
     """
-    assert isinstance(S, pd.DataFrame), 'S must be of type pd.DataFrame'
+    if not isinstance(S, pd.DataFrame):
+        raise TypeError(f'S must be of type pd.DataFrame. Got: {type(S)}')
 
     # Eq 15 in IEC 62600-101
     m0 = frequency_moment(S,0,frequency_bins=frequency_bins).squeeze() # convert to Series for calculation
@@ -413,7 +437,8 @@ def average_crest_period(S,frequency_bins=None):
         Average wave period [s] indexed by S.columns
 
     """
-    assert isinstance(S, pd.DataFrame), 'S must be of type pd.DataFrame'
+    if not isinstance(S, pd.DataFrame):
+        raise TypeError(f'S must be of type pd.DataFrame. Got: {type(S)}')
 
     m2 = frequency_moment(S,2,frequency_bins=frequency_bins).squeeze() # convert to Series for calculation
     m4 = frequency_moment(S,4,frequency_bins=frequency_bins).squeeze()
@@ -440,7 +465,8 @@ def average_wave_period(S,frequency_bins=None):
     Tm: pandas DataFrame
         Mean wave period [s] indexed by S.columns
     """
-    assert isinstance(S, pd.DataFrame), 'S must be of type pd.DataFrame'
+    if not isinstance(S, pd.DataFrame):
+        raise TypeError(f'S must be of type pd.DataFrame. Got: {type(S)}')
 
     m0 = frequency_moment(S,0,frequency_bins=frequency_bins).squeeze() # convert to Series for calculation
     m1 = frequency_moment(S,1,frequency_bins=frequency_bins).squeeze()
@@ -465,7 +491,8 @@ def peak_period(S):
     Tp: pandas DataFrame
         Wave peak period [s] indexed by S.columns
     """
-    assert isinstance(S, pd.DataFrame), 'S must be of type pd.DataFrame'
+    if not isinstance(S, pd.DataFrame):
+        raise TypeError(f'S must be of type pd.DataFrame. Got: {type(S)}')
 
     # Eq 14 in IEC 62600-101
     fp = S.idxmax(axis=0) # Hz
@@ -493,7 +520,8 @@ def energy_period(S,frequency_bins=None):
         Wave energy period [s] indexed by S.columns
     """
 
-    assert isinstance(S, (pd.Series,pd.DataFrame)), 'S must be of type pd.DataFrame or pd.Series'
+    if not isinstance(S, (pd.Series,pd.DataFrame)):
+        raise TypeError(f'S must be of type pd.DataFrame or pd.Series. Got: {type(S)}')
 
     mn1 = frequency_moment(S,-1,frequency_bins=frequency_bins).squeeze() # convert to Series for calculation
     m0  = frequency_moment(S,0,frequency_bins=frequency_bins).squeeze()
@@ -525,7 +553,8 @@ def spectral_bandwidth(S,frequency_bins=None):
     e: pandas DataFrame
         Spectral bandwidth [s] indexed by S.columns
     """
-    assert isinstance(S, pd.DataFrame), 'S must be of type pd.DataFrame'
+    if not isinstance(S, pd.DataFrame):
+        raise TypeError(f'S must be of type pd.DataFrame. Got: {type(S)}')
 
     m2 = frequency_moment(S,2,frequency_bins=frequency_bins).squeeze() # convert to Series for calculation
     m0 = frequency_moment(S,0,frequency_bins=frequency_bins).squeeze()
@@ -553,7 +582,8 @@ def spectral_width(S,frequency_bins=None):
     v: pandas DataFrame
         Spectral width [m] indexed by S.columns
     """
-    assert isinstance(S, pd.DataFrame), 'S must be of type pd.DataFrame'
+    if not isinstance(S, pd.DataFrame):
+        raise TypeError(f'S must be of type pd.DataFrame. Got: {type(S)}')
 
     mn2 = frequency_moment(S,-2,frequency_bins=frequency_bins).squeeze() # convert to Series for calculation
     m0 = frequency_moment(S,0,frequency_bins=frequency_bins).squeeze()
@@ -594,12 +624,18 @@ def energy_flux(S, h, deep=False, rho=1025, g=9.80665, ratio=2):
     J: pandas DataFrame
         Omni-directional wave energy flux [W/m] indexed by S.columns
     """
-    assert isinstance(S, (pd.Series,pd.DataFrame)), 'S must be of type pd.DataFrame or pd.Series'
-    assert isinstance(h, (int,float)), 'h must be of type int or float'
-    assert isinstance(deep, bool), 'deep must be of type bool'
-    assert isinstance(rho, (int,float)), 'rho must be of type int or float'
-    assert isinstance(g, (int,float)), 'g must be of type int or float'
-    assert isinstance(ratio, (int,float)), 'ratio must be of type int or float'
+    if not isinstance(S, (pd.Series,pd.DataFrame)):
+        raise TypeError(f'S must be of type pd.DataFrame or pd.Series. Got: {type(S)}')
+    if not isinstance(h, (int,float)):
+        raise TypeError(f'h must be of type int or float. Got: {type(h)}')
+    if not isinstance(deep, bool):
+        raise TypeError(f'deep must be of type bool. Got: {type(deep)}')
+    if not isinstance(rho, (int,float)):
+        raise TypeError(f'rho must be of type int or float. Got: {type(rho)}')
+    if not isinstance(g, (int,float)):
+        raise TypeError(f'g must be of type int or float. Got: {type(g)}')
+    if not isinstance(ratio, (int,float)):
+        raise TypeError(f'ratio must be of type int or float. Got: {type(ratio)}')
 
     if deep:
         # Eq 8 in IEC 62600-100, deep water simpilification
@@ -661,8 +697,10 @@ def energy_period_to_peak_period(Te, gamma):
     Tp: float or array
         Spectral peak period [s]
     """
-    assert isinstance(Te, (float, np.ndarray)), 'Te must be a float or a ndarray'
-    assert isinstance(gamma, (float, int)), 'gamma must be of type float or int'
+    if not isinstance(Te, (float, np.ndarray)):
+        raise TypeError(f'Te must be a float or a ndarray. Got: {type(Te)}')
+    if not isinstance(gamma, (float, int)):
+        raise TypeError(f'gamma must be of type float or int. Got: {type(gamma)}')
 
     factor = 0.8255 + 0.03852*gamma - 0.005537*gamma**2 + 0.0003154*gamma**3
 
@@ -694,12 +732,16 @@ def wave_celerity(k, h, g=9.80665, depth_check=False, ratio=2):
     """
     if isinstance(k, pd.DataFrame):
         k = k.squeeze()
-
-    assert isinstance(k, pd.Series), 'S must be of type pd.Series'
-    assert isinstance(h, (int,float)), 'h must be of type int or float'
-    assert isinstance(g, (int,float)), 'g must be of type int or float'
-    assert isinstance(depth_check, bool), 'depth_check must be of type bool'
-    assert isinstance(ratio, (int,float)), 'ratio must be of type int or float'
+    if not isinstance(k, (pd.Series, pd.DataFrame)):
+        raise TypeError(f'k must be of type pd.Series or pd.DataFrame. Got: {type(k)}')
+    if not isinstance(h, (int,float)):
+        raise TypeError(f'h must be of type int or float. Got: {type(h)}')
+    if not isinstance(g, (int,float)):
+        raise TypeError(f'g must be of type int or float. Got: {type(g)}')
+    if not isinstance(depth_check, bool):
+        raise TypeError(f'depth_check must be of type bool. Got: {type(depth_check)}')
+    if not isinstance(ratio, (int,float)):
+        raise TypeError(f'ratio must be of type int or float. Got: {type(ratio)}')
 
     f = k.index
     k = k.values
@@ -749,6 +791,9 @@ def wave_length(k):
     l: float or array
         Wave length [m] indexed by frequency
     """
+    if not isinstance(k, (int, float, list, np.ndarray, pd.DataFrame, pd.Series)):
+        raise TypeError(f'k must be of type int, float, list, np.ndarray, pd.DataFrame, or pd.Series. Got: {type(k)}')
+    
     if isinstance(k, (int, float, list)):
         k = np.array(k)
     elif isinstance(k, pd.DataFrame):
@@ -756,8 +801,6 @@ def wave_length(k):
     elif isinstance(k, pd.Series):
         k = k.values
 
-    assert isinstance(k, np.ndarray), 'k must be array-like'
-
     l = 2*np.pi/k
 
     return l
@@ -790,10 +833,14 @@ def wave_number(f, h, rho=1025, g=9.80665):
         f = np.atleast_1d(np.array(f))
     except:
         pass
-    assert isinstance(f, np.ndarray), 'f must be of type np.ndarray'
-    assert isinstance(h, (int,float)), 'h must be of type int or float'
-    assert isinstance(rho, (int,float)), 'rho must be of type int or float'
-    assert isinstance(g, (int,float)), 'g must be of type int or float'
+    if not isinstance(f, np.ndarray):
+        raise TypeError(f'f must be of type np.ndarray. Got: {type(f)}')
+    if not isinstance(h, (int,float)):
+        raise TypeError(f'h must be of type int or float. Got: {type(h)}')
+    if not isinstance(rho, (int,float)):
+        raise TypeError(f'rho must be of type int or float. Got: {type(rho)}')
+    if not isinstance(g, (int,float)):
+        raise TypeError(f'g must be of type int or float. Got: {type(g)}')
 
     w = 2*np.pi*f # angular frequency
     xi = w/np.sqrt(g/h) # note: =h*wa/sqrt(h*g/h)
@@ -811,7 +858,8 @@ def func(kk):
         w = w[mask]
 
         k, info, ier, mesg = _fsolve(func, k0_mask, full_output=True)
-        assert ier == 1, 'Wave number not found. ' + mesg
+        if not ier == 1:
+            raise ValueError('Wave number not found. ' + mesg)
         k0[mask] = k
 
     k = pd.DataFrame(k0, index=f, columns=['k'])
@@ -845,7 +893,11 @@ def depth_regime(l, h, ratio=2):
     depth_reg: boolean or boolean array
         Boolean True if deep water, False otherwise
     '''
-
+    if not isinstance(l, (int, float, list, np.ndarray, pd.DataFrame, pd.Series)):
+        raise TypeError(f'l must be of type int, float, list, np.ndarray, pd.DataFrame, or pd.Series. Got: {type(l)}')
+    if not isinstance(h, (int, float)):
+        raise TypeError(f'h must be of type int or float. Got: {type(h)}')
+    
     if isinstance(l, (int, float, list)):
         l = np.array(l)
     elif isinstance(l, pd.DataFrame):
@@ -853,9 +905,6 @@ def depth_regime(l, h, ratio=2):
     elif isinstance(l, pd.Series):
         l = l.values
 
-    assert isinstance(l, (np.ndarray)), "l must be array-like"
-    assert isinstance(h, (int, float)), "h must be of type int or float"
-
     depth_reg = h/l > ratio
 
     return  depth_reg


Test patch:
diff --git a/mhkit/tests/dolfyn/test_rotate_adv.py b/mhkit/tests/dolfyn/test_rotate_adv.py
index 26ab76d3c..c67f42a2a 100644
--- a/mhkit/tests/dolfyn/test_rotate_adv.py
+++ b/mhkit/tests/dolfyn/test_rotate_adv.py
@@ -167,7 +167,6 @@ def test_rotate_warnings(self):
             rotate2(warn2, 'earth')
         with self.assertRaises(Exception):
             set_inst2head_rotmat(warn3, np.eye(3))
-            set_inst2head_rotmat(warn3, np.eye(3))
         with self.assertRaises(Exception):
             set_inst2head_rotmat(warn4, np.eye(3))
 
diff --git a/mhkit/tests/wave/io/test_ndbc.py b/mhkit/tests/wave/io/test_ndbc.py
index dc1d7a6c7..610dc2977 100644
--- a/mhkit/tests/wave/io/test_ndbc.py
+++ b/mhkit/tests/wave/io/test_ndbc.py
@@ -128,7 +128,7 @@ def test_ndbc_request_data_from_dataframe(self):
         assert_frame_equal(self.swden, ndbc_data['1996'])
 
     def test_ndbc_request_data_filenames_length(self):
-        with self.assertRaises(AssertionError):
+        with self.assertRaises(ValueError):
             wave.io.ndbc.request_data('swden', pd.Series(dtype=float))
 
     def test_ndbc_to_datetime_index(self):
@@ -214,7 +214,7 @@ def test_plot_directional_spectrum(self):
             self.directional_data, directions)
         wave.graphics.plot_directional_spectrum(
             spectrum,
-            min=0.0,
+            color_level_min=0.0,
             fill=True,
             nlevels=6,
             name="Elevation Variance",