refactor linke turbidity lookup

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

@@ -1,8 +1,6 @@
 import logging
 pvl_logger = logging.getLogger('pvlib')
 
-import datetime
-
 import numpy as np
 import pandas as pd
 
@@ -15,17 +13,12 @@
 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():
@@ -38,7 +31,7 @@ def test_ineichen_supply_linke():
 
 def test_ineichen_solpos():
     clearsky.ineichen(times, tus, linke_turbidity=3,
-                      solarposition_method='pyephem')
+                      solarposition_method='nrel_numpy')
 
 def test_ineichen_airmass():
     clearsky.ineichen(times, tus, linke_turbidity=3,
@@ -50,6 +43,9 @@ def test_ineichen_keys():
     assert 'dni' in clearsky_data.columns
     assert 'dhi' in clearsky_data.columns
 
+def test_lookup_linke_turbidity():
+    raise Exception
+
 # test the haurwitz clear sky implementation
 def test_haurwitz():
     clearsky.haurwitz(ephem_data['zenith'])