fix pvlib#263 and pvlib#262

* set minimum index as 0 not 1
* adjust _linearly_scale so that is increments are 5 minutes or 1/12 of an
 arcdegree from the centers of neighboring indices
* this fixes the out of bounds error for latitudes or longitudes at the
 limits of the range
* move _linearly_scale to be right after lookup_linke_turbidity
* add link to pdf on AOD and LT from MeteoTest
* raise IndexError if outputmatrix is more than half an index outside of
 the limits

Signed-off-by: Mark Mikofski <[email protected]>

Author: mikofski

pvlib/clearsky.py

@@ -5,6 +5,8 @@
 
 from __future__ import division
 
+import logging
+
 import os
 from collections import OrderedDict
 
@@ -13,6 +15,9 @@
 
 from pvlib import tools
 
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.DEBUG)
+
 
 def ineichen(apparent_zenith, airmass_absolute, linke_turbidity,
              altitude=0, dni_extra=1364.):
@@ -186,6 +191,12 @@ def lookup_linke_turbidity(time, latitude, longitude, filepath=None,
     # so divide the number from the file by 20 to get the
     # turbidity.
 
+    # The nodes of the grid are 5' (1/12=0.0833[arcdeg]) apart.
+    # From Section 8 of Aerosol optical depth and Linke turbidity climatology
+    # http://www.meteonorm.com/images/uploads/downloads/ieashc36_report_TL_AOD_climatologies.pdf
+    # 1st row: 89.9583 S, 2nd row: 89.875 S
+    # 1st column: 179.9583 W, 2nd column: 179.875 W
+
     try:
         import scipy.io
     except ImportError:
@@ -201,10 +212,10 @@ def lookup_linke_turbidity(time, latitude, longitude, filepath=None,
     linke_turbidity_table = mat['LinkeTurbidity']
 
     latitude_index = (
-        np.around(_linearly_scale(latitude, 90, -90, 1, 2160))
+        np.around(_linearly_scale(latitude, 90, -90, 0, 2160))
         .astype(np.int64))
     longitude_index = (
-        np.around(_linearly_scale(longitude, -180, 180, 1, 4320))
+        np.around(_linearly_scale(longitude, -180, 180, 0, 4320))
         .astype(np.int64))
 
     g = linke_turbidity_table[latitude_index][longitude_index]
@@ -232,6 +243,33 @@ def lookup_linke_turbidity(time, latitude, longitude, filepath=None,
     return linke_turbidity
 
 
+def _linearly_scale(inputmatrix, inputmin, inputmax, outputmin, outputmax):
+    """
+    Used by linke turbidity lookup function.
+    """
+    inputrange = inputmax - inputmin
+    outputrange = outputmax - outputmin
+    delta = outputrange/inputrange  # number of indices per input unit
+    inputmin = inputmin + 1.0 / delta / 2.0  # shift index to center
+    outputmax = outputmax - 1
+    outputmatrix = (inputmatrix - inputmin) * delta + outputmin
+    LOGGER.debug('outputmatrix: %g', outputmatrix)
+    # TODO: allow wrapping for longitude
+    err = IndexError('Input, %g, is out of range (%g, %g).' %
+                     (inputmatrix, inputmax - inputrange, inputmax))
+    if outputmatrix > outputmax:
+        if np.around(outputmatrix - outputmax, 1) <= 0.5:
+            outputmatrix = outputmax
+        else:
+            raise err
+    elif outputmatrix < outputmin:
+        if np.around(outputmin - outputmatrix, 1) <= 0.5:
+            outputmatrix = outputmin
+        else:
+            raise err
+    return outputmatrix
+
+
 def haurwitz(apparent_zenith):
     '''
     Determine clear sky GHI from Haurwitz model.
@@ -281,15 +319,6 @@ def haurwitz(apparent_zenith):
     return df_out
 
 
-def _linearly_scale(inputmatrix, inputmin, inputmax, outputmin, outputmax):
-    """ used by linke turbidity lookup function """
-
-    inputrange = inputmax - inputmin
-    outputrange = outputmax - outputmin
-    outputmatrix = (inputmatrix-inputmin) * outputrange/inputrange + outputmin
-    return outputmatrix
-
-
 def simplified_solis(apparent_elevation, aod700=0.1, precipitable_water=1.,
                      pressure=101325., dni_extra=1364.):
     """

pvlib/test/test_clearsky.py

@@ -1,5 +1,7 @@
 from collections import OrderedDict
 
+import logging
+
 import numpy as np
 import pandas as pd
 
@@ -14,6 +16,9 @@
 
 from conftest import requires_scipy
 
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.DEBUG)
+
 
 def test_ineichen_series():
     tus = Location(32.2, -111, 'US/Arizona', 700)
@@ -440,3 +445,46 @@ def test_simplified_solis_nans_series():
                                     precipitable_water, pressure, dni_extra)
 
     assert_frame_equal(expected, out)
+
+
+def test_linke_turbidity_corners():
+    """
+    Test Linke turbidity corners out of bounds
+    """
+    time = pd.DatetimeIndex('%d/1/2016' % (m + 1) for m in range(12))
+    # Northwest
+    assert np.allclose(
+        clearsky.lookup_linke_turbidity(time, 90, -180, interp_turbidity=False),
+        [1.9, 1.9, 1.9, 2.0, 2.05, 2.05, 2.1, 2.1, 2.0, 1.95, 1.9, 1.9])
+    # Southwest
+    assert np.allclose(
+        clearsky.lookup_linke_turbidity(time, -90, -180, interp_turbidity=False),
+        [1.35, 1.3, 1.45, 1.35, 1.35, 1.35, 1.35, 1.35, 1.35, 1.4, 1.4, 1.3])
+    # Northeast
+    assert np.allclose(
+        clearsky.lookup_linke_turbidity(time, 90, 180, interp_turbidity=False),
+        [1.9, 1.9, 1.9, 2.0, 2.05, 2.05, 2.1, 2.1, 2.0, 1.95, 1.9, 1.9])
+    # Southeast
+    assert np.allclose(
+        clearsky.lookup_linke_turbidity(time, -90, 180, interp_turbidity=False),
+        [1.35, 1.7, 1.35, 1.35, 1.35, 1.35, 1.35, 1.35, 1.35, 1.35, 1.35, 1.7])
+    try:
+        clearsky.lookup_linke_turbidity(time, 91, -122, interp_turbidity=False)
+    except IndexError as err:
+        LOGGER.exception(err)
+        assert isinstance(err, IndexError)
+    try:
+        clearsky.lookup_linke_turbidity(time, 38.2, 181, interp_turbidity=False)
+    except IndexError as err:
+        LOGGER.exception(err)
+        assert isinstance(err, IndexError)
+    try:
+        clearsky.lookup_linke_turbidity(time, -91, -122, interp_turbidity=False)
+    except IndexError as err:
+        LOGGER.exception(err)
+        assert isinstance(err, IndexError)
+    try:
+        clearsky.lookup_linke_turbidity(time, 38.2, -181, interp_turbidity=False)
+    except IndexError as err:
+        LOGGER.exception(err)
+        assert isinstance(err, IndexError)