Merge pull request #98 from wholmgren/refactor-turbidity

Refactor clearsky.py

Author: wholmgren

docs/sphinx/source/whatsnew/v0.2.2.txt

@@ -10,11 +10,15 @@ Enhancements
 ~~~~~~~~~~~~
 
 * Adds Python 3.5 compatibility (:issue:`87`)
+* Moves the Linke turbidity lookup into ``clearsky.lookup_linke_turbidity``.
+  The API for ``clearsky.ineichen`` remains the same. (:issue:`95`)
 
 
 Bug fixes
 ~~~~~~~~~
 
+* Fixes an import and KeyError in the IPython notebook tutorials
+  (:issue:`94`).
 * Uses the ``logging`` module properly by replacing ``format``
   calls with ``args``. This results in a 5x speed increase for
   ``tracking.singleaxis`` (:issue:`89`).
@@ -24,4 +28,5 @@ Contributors
 ~~~~~~~~~~~~
 
 * Will Holmgren
+* jetheurer
 * dacoex

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

@@ -1,59 +1,137 @@
 import logging
 pvl_logger = logging.getLogger('pvlib')
 
-import datetime
-
 import numpy as np
 import pandas as pd
 
 from nose.tools import raises
-
 from numpy.testing import assert_almost_equal
+from pandas.util.testing import assert_frame_equal, assert_series_equal
 
 from pvlib.location import Location
 from pvlib import clearsky
 from pvlib import solarposition
 
 # setup times and location to be tested.
-times = pd.date_range(start=datetime.datetime(2014,6,24), 
-                      end=datetime.datetime(2014,6,26), freq='1Min')
-
 tus = Location(32.2, -111, 'US/Arizona', 700)
-
+times = pd.date_range(start='2014-06-24', end='2014-06-25', freq='3h')
 times_localized = times.tz_localize(tus.tz)
 
 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
-    clearsky.ineichen(times, tus)
+    expected = pd.DataFrame(np.array([[0.,0.,0.],
+                                      [0.,0.,0.],
+                                      [40.53660309,302.47614235,78.1470311],
+                                      [98.88372629,865.98938602,699.93403875],
+                                      [122.57870881,931.83716051,1038.62116584],
+                                      [109.30270612,899.88002304,847.68806472],
+                                      [64.25699595,629.91187925,254.53048144],
+                                      [0.,0.,0.],
+                                      [0.,0.,0.]]),
+                            columns=['dhi', 'dni', 'ghi'],
+                            index=times_localized)
+    out = clearsky.ineichen(times, tus)
+    assert_frame_equal(expected, out)
 
+
 def test_ineichen_supply_linke():
-    clearsky.ineichen(times, tus, linke_turbidity=3)
+    expected = pd.DataFrame(np.array([[0.,0.,0.],
+                                      [0.,0.,0.],
+                                      [40.18673553,322.0649964,80.23287692],
+                                      [95.14405816,876.49507151,703.48596755],
+                                      [118.45873721,939.81653473,1042.34531752],
+                                      [105.36671577,909.113377,851.3283881],
+                                      [61.91607984,647.40869542,257.47471759],
+                                      [0.,0.,0.],
+                                      [0.,0.,0.]]),
+                            columns=['dhi', 'dni', 'ghi'],
+                            index=times_localized)
+    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='pyephem')
+                      solarposition_method='ephemeris')
+
 
 def test_ineichen_airmass():
-    clearsky.ineichen(times, tus, linke_turbidity=3,
-                      airmass_model='simple')
+    expected = pd.DataFrame(np.array([[0.,0.,0.],
+                                      [0.,0.,0.],
+                                      [41.70761136,293.72203458,78.22953786],
+                                      [95.20590465,876.1650047,703.31872722],
+                                      [118.46089555,939.8078753,1042.33896321],
+                                      [105.39577655,908.97804342,851.24640259],
+                                      [62.35382269,642.91022293,256.55363539],
+                                      [0.,0.,0.],
+                                      [0.,0.,0.]]),
+                            columns=['dhi', 'dni', 'ghi'],
+                            index=times_localized)
+    out = clearsky.ineichen(times, tus, linke_turbidity=3,
+                            airmass_model='simple')
+    assert_frame_equal(expected, out)
+
+
+def test_lookup_linke_turbidity():
+    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)
+
+
+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_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
 
-# 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
+    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)