Added support for array input to spectral.jonswap

docs/whatsnew.rst

@@ -36,7 +36,7 @@ New functions/methods
 Bug fixes
 ^^^^^^^^^
 
-* Store comments from SWAN spectrum files as single string intead of a
+* Store comments from SWAN spectrum files as single string instead of a
   list of strings as the scipy netcdf I/O cannot cope with lists of
   strings.
 
@@ -48,7 +48,12 @@ Bug fixes
 
 * Do not squeeze energy matrix as that may cause conflicts with
   dimensions of length unity. Instead, reshape the energy matrix to
-  the approriate size.
+  the appropriate size.
+
+* Fix a datatype bug in spectral.jonswap() that raised `ValueError:
+  Integers to negative integer powers are not allowed` in the case Hm0 or Tp
+  were integers. Also modify spectral.jonswap() to accept arrays as inputs
+  for Hm0 and Tp. Only one-dimensional arrays for now. - Caio Stringari
 
 Tests
 ^^^^^

oceanwaves/spectral.py

@@ -4,6 +4,9 @@
 
 from oceanwaves.utils import *
 
+import logging
+logger = logging.getLogger(__name__)
+
 
 def jonswap(f, Hm0, Tp, gamma=3.3, sigma_low=.07, sigma_high=.09,
             g=9.81, method='yamaguchi', normalize=True):
@@ -13,9 +16,9 @@ def jonswap(f, Hm0, Tp, gamma=3.3, sigma_low=.07, sigma_high=.09,
     ----------
     f : numpy.ndarray
         Array of frequencies
-    Hm0 : float
+    Hm0 : float, numpy.ndarray
         Required zeroth order moment wave height
-    Tp : float
+    Tp : float, numpy.ndarray
         Required peak wave period
     gamma : float
         JONSWAP peak-enhancement factor (default: 3.3)
@@ -29,14 +32,42 @@ def jonswap(f, Hm0, Tp, gamma=3.3, sigma_low=.07, sigma_high=.09,
         Method to compute alpha (default: yamaguchi)
     normalize : bool
         Normalize resulting spectrum to match ``Hm0``
-    
+
     Returns
     -------
     E : numpy.ndarray
-        Array of shape ``f`` with wave energy densities 
-    
+        Array of shape ``f, Hm0.shape`` with wave energy densities
+
     '''
 
+    # C Stringari - 04/06/2018
+    # check input data types to avoid the following error:
+    # ValueError: Integers to negative integer powers are not allowed.
+
+    # raise an warning if the frequency array starts with zero. if the
+    # user gives an array with zeros, the output will be inf at that
+    # frequency
+    if 0.0 in f:
+        logger.warn('Frequency array contains zeros.')
+
+    # get the input dtypes and promote to float, if needed
+    f = ensure_float(f)
+    Hm0 = ensure_float(Hm0)
+    Tp = ensure_float(Tp)
+
+    # check shapes of Hm0 and Tp, raise an error if the don't match
+    if isinstance(Hm0, np.ndarray):
+        if isinstance(Tp, np.ndarray):
+            if Hm0.shape != Tp.shape:
+                raise ValueError("Dimensions of Hm0 and Tp should match.")
+
+    # This is a very naive implementation to deal with array inputs,
+    # but will work if Hm0 and Tp are vectors.
+    if isinstance(Hm0, np.ndarray):
+        f = f[:, np.newaxis].repeat(len(Hm0), axis=1)
+        Hm0 = Hm0[np.newaxis, :].repeat(len(f), axis=0)
+        Tp = Tp[np.newaxis, :].repeat(len(f), axis=0)
+
     # Pierson-Moskowitz
     if method.lower() == 'yamaguchi':
         alpha = 1. / (.06533 * gamma ** .8015 + .13467) / 16.
@@ -45,17 +76,18 @@ def jonswap(f, Hm0, Tp, gamma=3.3, sigma_low=.07, sigma_high=.09,
     else:
         raise ValueError('Unknown method: %s' % method)
 
-    E_pm  = alpha * Hm0**2 * Tp**-4 * f**-5 * np.exp(-1.25 * (Tp * f)**-4)
-        
+    E_pm = alpha * Hm0**2 * Tp**-4 * f**-5 * np.exp(-1.25 * (Tp * f)**-4)
+
     # JONSWAP
     sigma = np.ones(f.shape) * sigma_low
     sigma[f > 1./Tp] = sigma_high
 
-    E_js = E_pm * gamma**np.exp(-0.5 * (Tp * f - 1)**2. / sigma**2.);
-    
+    E_js = E_pm * gamma**np.exp(-0.5 * (Tp * f - 1)**2. / sigma**2.)
+
     if normalize:
-        E_js *= Hm0**2. / (16. * trapz_and_repeat(E_js, f, axis=-1))
-
+        # axis=0 seems to work fine with all kinds of inputs
+        E_js *= Hm0**2. / (16. * trapz_and_repeat(E_js, f, axis=0))
+
     return E_js
 
 
@@ -80,11 +112,10 @@ def directional_spreading(theta, theta_peak=0., s=20., units='deg',
     -------
     p_theta : numpy.ndarray
        Array of directional weights
-
     '''
-    
+
     from math import gamma
-    
+
     theta = np.asarray(theta, dtype=np.float)
 
     # convert units to radians
@@ -97,19 +128,49 @@ def directional_spreading(theta, theta_peak=0., s=20., units='deg',
         raise ValueError('Unknown units: %s')
 
     # compute directional spreading
-    #A1 = (2.**s) * (gamma(s / 2 + 1))**2. / (np.pi * gamma(s + 1))
-    #p_theta = A1 * np.maximum(0., np.cos(theta - theta_peak))
+    # A1 = (2.**s) * (gamma(s / 2 + 1))**2. / (np.pi * gamma(s + 1))
+    # p_theta = A1 * np.maximum(0., np.cos(theta - theta_peak))
     p_theta = np.maximum(0., np.cos(theta - theta_peak))**s
 
     # convert to original units
     if units.lower().startswith('deg'):
         theta = np.degrees(theta)
         theta_peak = np.degrees(theta_peak)
         p_theta = np.degrees(p_theta)
-    
+
     # normalize directional spreading
     if normalize:
         p_theta /= trapz_and_repeat(p_theta, theta - theta_peak, axis=-1)
 
     return p_theta
 
+
+def ensure_float(var):
+    '''
+    Auxiliary function to detect and fix dtypes, if needed
+
+    Parameters
+    ----------
+    var : anything
+        Array to check
+
+    Returns
+    -------
+    var : numpy.ndarray
+       The same as the input but either as a float or an array of floats
+    '''
+    # fist, if it's a list, convert to a numpy.ndarray
+    if isinstance(var, list):
+        var = np.array(var)
+    # if it's an np.ndarray(), make sure it's a array of floats
+    elif isinstance(var, np.ndarray):
+        if var.dtype != np.dtype('float'):
+            var = var.astype(np.float)
+    # if it's a float, well, do nothing
+    elif isinstance(var, float):
+        var = var
+    # unknown data type
+    else:
+        raise ValueError("Data type could not be detected automatically.")
+    # allways return a float or array of floats
+    return var

tests/test_spectral.py

@@ -5,6 +5,7 @@
 FREQ = np.arange(0, 1, .05)[1:]
 THETA = np.arange(0, 180, 5)
 
+
 def test_jonswap_1():
     '''Test computation of jonswap spectrum following Yamaguchi'''
     E = jonswap(FREQ, Hm0=1., Tp=4., method='yamaguchi',
@@ -20,6 +21,14 @@ def test_jonswap_2():
     assert_almost_equals(FREQ[np.argmax(E)], 1/4.)
 
 
+def test_jonswap_with_arrays():
+    '''Test computation of jonswap spectrum using arrays as input'''
+    E = jonswap(FREQ, Hm0=np.array([1.0, 1.0]), Tp=np.array([4.0, 4.0]),
+                method='goda', normalize=True)
+    assert_equal(E.ndim, 2)
+    assert_equal(E.shape[1], 2)
+
+
 @raises(ValueError)
 def test_jonswap_exception1():
     '''Test exception if unsupported spectrum method is requested'''