add the new lookup code. add more tests

Author: wholmgren

pvlib/clearsky.py

@@ -115,45 +115,9 @@ def ineichen(time, location, linke_turbidity=None,
 
 
     if linke_turbidity is None:
-        # The .mat file 'LinkeTurbidities.mat' contains a single 2160 x 4320 x 12
-        # matrix of type uint8 called 'LinkeTurbidity'. The rows represent global
-        # latitudes from 90 to -90 degrees; the columns represent global longitudes
-        # from -180 to 180; and the depth (third dimension) represents months of
-        # the year from January (1) to December (12). To determine the Linke
-        # turbidity for a position on the Earth's surface for a given month do the
-        # following: LT = LinkeTurbidity(LatitudeIndex, LongitudeIndex, month).  
-        # Note that the numbers within the matrix are 20 * Linke Turbidity, 
-        # so divide the number from the file by 20 to get the
-        # turbidity.
-        
-        try:
-            import scipy.io
-        except ImportError:
-            raise ImportError('The Linke turbidity lookup table requires scipy. ' +
-                              'You can still use clearsky.ineichen if you ' +
-                              'supply your own turbidities.')
-        
-        # consider putting this code at module level
-        this_path = os.path.dirname(os.path.abspath(__file__))
-        logger.debug('this_path=%s', this_path)
-        
-        mat = scipy.io.loadmat(os.path.join(this_path, 'data', 'LinkeTurbidities.mat'))
-        linke_turbidity = mat['LinkeTurbidity']
-        LatitudeIndex = np.round_(_linearly_scale(location.latitude,90,- 90,1,2160))
-        LongitudeIndex = np.round_(_linearly_scale(location.longitude,- 180,180,1,4320))
-        g = linke_turbidity[LatitudeIndex][LongitudeIndex]
-        if interp_turbidity:
-            logger.info('interpolating turbidity to the day')
-            g2 = np.concatenate([[g[-1]], g, [g[0]]]) # wrap ends around
-            days = np.linspace(-15, 380, num=14) # map day of year onto month (approximate)
-            LT = pd.Series(np.interp(time.dayofyear, days, g2), index=time)
-        else:
-            logger.info('using monthly turbidity')
-            ApplyMonth = lambda x:g[x[0]-1]
-            LT = pd.DataFrame(time.month, index=time)
-            LT = LT.apply(ApplyMonth, axis=1)
-        TL = LT / 20.
-        logger.info('using TL=\n%s', TL)
+        TL = lookup_linke_turbidity(time, location.latitude,
+                                    location.longitude,
+                                    interp_turbidity=interp_turbidity)
     else:
         TL = linke_turbidity
 
@@ -221,6 +185,87 @@ def ineichen(time, location, linke_turbidity=None,
     return df_out
 
 
+def lookup_linke_turbidity(time, latitude, longitude, filepath=None,
+                           interp_turbidity=True):
+    """
+    Look up the Linke Turibidity from the ``LinkeTurbidities.mat``
+    data file supplied with pvlib.
+
+    Parameters
+    ----------
+    time : pandas.DatetimeIndex
+
+    latitude : float
+
+    longitude : float
+
+    filepath : string
+        The path to the ``.mat`` file.
+
+    interp_turbidity : bool
+        If ``True``, interpolates the monthly Linke turbidity values
+        found in ``LinkeTurbidities.mat`` to daily values.
+
+    Returns
+    -------
+    turbidity : Series
+    """
+
+    # The .mat file 'LinkeTurbidities.mat' contains a single 2160 x 4320 x 12
+    # matrix of type uint8 called 'LinkeTurbidity'. The rows represent global
+    # latitudes from 90 to -90 degrees; the columns represent global longitudes
+    # from -180 to 180; and the depth (third dimension) represents months of
+    # the year from January (1) to December (12). To determine the Linke
+    # turbidity for a position on the Earth's surface for a given month do the
+    # following: LT = LinkeTurbidity(LatitudeIndex, LongitudeIndex, month).  
+    # Note that the numbers within the matrix are 20 * Linke Turbidity, 
+    # so divide the number from the file by 20 to get the
+    # turbidity. 
+    
+    try:
+        import scipy.io
+    except ImportError:
+        raise ImportError('The Linke turbidity lookup table requires scipy. ' +
+                          'You can still use clearsky.ineichen if you ' +
+                          'supply your own turbidities.')
+
+    if filepath is None:
+        pvlib_path = os.path.dirname(os.path.abspath(__file__))
+        filepath = os.path.join(pvlib_path, 'data', 'LinkeTurbidities.mat')
+
+    mat = scipy.io.loadmat(filepath)
+    linke_turbidity_table = mat['LinkeTurbidity']
+
+    latitude_index = np.around(_linearly_scale(latitude, 90, -90, 1, 2160))
+    longitude_index = np.around(_linearly_scale(longitude, -180, 180, 1, 4320))
+
+    g = linke_turbidity_table[latitude_index][longitude_index]
+
+    if interp_turbidity:
+        logger.info('interpolating turbidity to the day')
+        # Cata covers 1 year.
+        # Assume that data corresponds to the value at
+        # the middle of each month.
+        # This means that we need to add previous Dec and next Jan
+        # to the array so that the interpolation will work for 
+        # Jan 1 - Jan 15 and Dec 16 - Dec 31.
+        # Then we map the month value to the day of year value.
+        # This is approximate and could be made more accurate.
+        g2 = np.concatenate([[g[-1]], g, [g[0]]])
+        days = np.linspace(-15, 380, num=14)
+        linke_turbidity = pd.Series(np.interp(time.dayofyear, days, g2),
+                                    index=time)
+    else:
+        logger.info('using monthly turbidity')
+        apply_month = lambda x: g[x[0]-1]
+        linke_turbidity = pd.DataFrame(time.month, index=time)
+        linke_turbidity = linke_turbidity.apply(apply_month, axis=1)
+
+    linke_turbidity /= 20.
+
+    return linke_turbidity
+
+
 def haurwitz(apparent_zenith):
     '''
     Determine clear sky GHI from Haurwitz model.

pvlib/test/test_clearsky.py

@@ -6,7 +6,7 @@
 
 from nose.tools import raises
 from numpy.testing import assert_almost_equal
-from pandas.util.testing import assert_frame_equal
+from pandas.util.testing import assert_frame_equal, assert_series_equal
 
 from pvlib.location import Location
 from pvlib import clearsky
@@ -20,10 +20,8 @@
 ephem_data = solarposition.get_solarposition(times, tus)
 
 
-# test the ineichen clear sky model implementation in a few ways
-
 def test_ineichen_required():
-    # the clearsky function should lookup the linke turbidity on its own
+    # the clearsky function should call lookup_linke_turbidity by default
     # will fail without scipy
     expected = pd.DataFrame(np.array([[0.,0.,0.],
                                       [0.,0.,0.],
@@ -39,6 +37,7 @@ def test_ineichen_required():
     out = clearsky.ineichen(times, tus)
     assert_frame_equal(expected, out)
 
+
 def test_ineichen_supply_linke():
     expected = pd.DataFrame(np.array([[0.,0.,0.],
                                       [0.,0.,0.],
@@ -54,10 +53,12 @@ def test_ineichen_supply_linke():
     out = clearsky.ineichen(times, tus, linke_turbidity=3)
     assert_frame_equal(expected, out)
 
+
 def test_ineichen_solpos():
     clearsky.ineichen(times, tus, linke_turbidity=3,
                       solarposition_method='ephemeris')
 
+
 def test_ineichen_airmass():
     expected = pd.DataFrame(np.array([[0.,0.,0.],
                                       [0.,0.,0.],
@@ -74,19 +75,63 @@ def test_ineichen_airmass():
                             airmass_model='simple')
     assert_frame_equal(expected, out)
 
-def test_ineichen_keys():
-    clearsky_data = clearsky.ineichen(times, tus, linke_turbidity=3)
-    assert 'ghi' in clearsky_data.columns
-    assert 'dni' in clearsky_data.columns
-    assert 'dhi' in clearsky_data.columns
 
 def test_lookup_linke_turbidity():
-    raise Exception
+    times = pd.date_range(start='2014-06-24', end='2014-06-25',
+                          freq='12h', tz=tus.tz)
+    # expect same value on 2014-06-24 0000 and 1200, and
+    # diff value on 2014-06-25
+    expected = pd.Series(np.array([3.10126582, 3.10126582, 3.11443038]),
+                         index=times)
+    out = clearsky.lookup_linke_turbidity(times, tus.latitude, tus.longitude)
+    assert_series_equal(expected, out)
+
+
+def test_lookup_linke_turbidity_nointerp():
+    times = pd.date_range(start='2014-06-24', end='2014-06-25',
+                          freq='12h', tz=tus.tz)
+    # expect same value for all days
+    expected = pd.Series(np.array([3., 3., 3.]), index=times)
+    out = clearsky.lookup_linke_turbidity(times, tus.latitude, tus.longitude,
+                                          interp_turbidity=False)
+    assert_series_equal(expected, out)
 
-# test the haurwitz clear sky implementation
-def test_haurwitz():
-    clearsky.haurwitz(ephem_data['zenith'])
 
-def test_haurwitz_keys():
-    clearsky_data = clearsky.haurwitz(ephem_data['zenith'])
-    assert 'ghi' in clearsky_data.columns
+def test_lookup_linke_turbidity_months():
+    times = pd.date_range(start='2014-04-01', end='2014-07-01',
+                          freq='1M', tz=tus.tz)
+    expected = pd.Series(np.array([2.8943038, 2.97316456, 3.18025316]),
+                         index=times)
+    out = clearsky.lookup_linke_turbidity(times, tus.latitude,
+                                          tus.longitude)
+    assert_series_equal(expected, out)
+
+
+def test_lookup_linke_turbidity_nointerp_months():
+    times = pd.date_range(start='2014-04-10', end='2014-07-10',
+                          freq='1M', tz=tus.tz)
+    expected = pd.Series(np.array([2.85, 2.95, 3.]), index=times)
+    out = clearsky.lookup_linke_turbidity(times, tus.latitude, tus.longitude,
+                                          interp_turbidity=False)
+    assert_series_equal(expected, out)
+    # changing the dates shouldn't matter if interp=False
+    times = pd.date_range(start='2014-04-05', end='2014-07-05',
+                          freq='1M', tz=tus.tz)
+    out = clearsky.lookup_linke_turbidity(times, tus.latitude, tus.longitude,
+                                          interp_turbidity=False)
+    assert_series_equal(expected, out)
+
+
+def test_haurwitz():
+    expected = pd.DataFrame(np.array([[0.],
+                                      [0.],
+                                      [82.85934048],
+                                      [699.74514735],
+                                      [1016.50198354],
+                                      [838.32103769],
+                                      [271.90853863],
+                                      [0.],
+                                      [0.]]),
+                             columns=['ghi'], index=times_localized)
+    out = clearsky.haurwitz(ephem_data['zenith'])
+    assert_frame_equal(expected, out)