Add option to control airmass/perez enhancement factor in clearsky.ineichen (#459)

* add perez_enhancement option

* pass kwargs through location.get_clearsky to ineichen

* hard code test values to avoid test failures

* refactor test_irradiance.py input data into fixtures

* refactor location tests to avoid tight cs dependence

* mock/refactor all tests except losses

* fix typo in test

* update whatsnew

* update test ineichen values for perez_enhancement=False

* update test_clearsky values

* add test for perez enhancement

* update tests for py2.7-3.5

* add requires_scipy for min env test

* update comment to point to correct refs

Author: wholmgren

docs/sphinx/source/whatsnew/v0.6.0.rst

@@ -11,6 +11,12 @@ API Changes
   support from PVSystem.pvwatts_losses. Enables custom losses specification
   in ModelChain calculations. (:issue:`484`)
 * removed irradiance parameter from ModelChain.run_model and ModelChain.prepare_inputs
+* Add ``perez_enhancement`` keyword argument to clearsky.ineichen to control
+  whether or not the "perez enhancement factor" is applied. The enhancement
+  factor was always applied until now. Now it is turned off by default. The
+  enhancement factor can yield unphysical results, especially for latitudes
+  closer to the poles and especially in the winter months. It may yield
+  improved results under other conditions. (:issue:`435`)
 
 
 Enhancements
@@ -81,8 +87,11 @@ Documentation
 Testing
 ~~~~~~~
 * Add pytest-mock dependency
-* Use pytest-mock to ensure that PVSystem methods call corresponding functions
-  correctly. Removes implicit dependence on precise return values of functions
+* Use pytest-mock to ensure that PVSystem and ModelChain methods call
+  corresponding functions correctly. Removes implicit dependence on precise
+  return values of some function/methods. (:issue:`394`)
+* Additional test refactoring to limit test result dependence to a single
+  function per test. (:issue:`394`)
 * Use pytest-mock to ensure that ModelChain DC model is set up correctly.
 
 

pvlib/clearsky.py

@@ -16,7 +16,7 @@
 
 
 def ineichen(apparent_zenith, airmass_absolute, linke_turbidity,
-             altitude=0, dni_extra=1364.):
+             altitude=0, dni_extra=1364., perez_enhancement=False):
     '''
     Determine clear sky GHI, DNI, and DHI from Ineichen/Perez model.
 
@@ -47,6 +47,12 @@ def ineichen(apparent_zenith, airmass_absolute, linke_turbidity,
         Extraterrestrial irradiance. The units of ``dni_extra``
         determine the units of the output.
 
+    perez_enhancement : bool, default False
+        Controls if the Perez enhancement factor should be applied.
+        Setting to True may produce spurious results for times when
+        the Sun is near the horizon and the airmass is high.
+        See https://github.com/pvlib/pvlib-python/issues/435
+
     Returns
     -------
     clearsky : DataFrame (if Series input) or OrderedDict of arrays
@@ -79,28 +85,9 @@ def ineichen(apparent_zenith, airmass_absolute, linke_turbidity,
         ISES Solar World Congress, June 2003. Goteborg, Sweden.
     '''
 
-    # Dan's note on the TL correction: By my reading of the publication
-    # on pages 151-157, Ineichen and Perez introduce (among other
-    # things) three things. 1) Beam model in eqn. 8, 2) new turbidity
-    # factor in eqn 9 and appendix A, and 3) Global horizontal model in
-    # eqn. 11. They do NOT appear to use the new turbidity factor (item
-    # 2 above) in either the beam or GHI models. The phrasing of
-    # appendix A seems as if there are two separate corrections, the
-    # first correction is used to correct the beam/GHI models, and the
-    # second correction is used to correct the revised turibidity
-    # factor. In my estimation, there is no need to correct the
-    # turbidity factor used in the beam/GHI models.
-
-    # Create the corrected TL for TL < 2
-    # TLcorr = TL;
-    # TLcorr(TL < 2) = TLcorr(TL < 2) - 0.25 .* (2-TLcorr(TL < 2)) .^ (0.5);
-
-    # This equation is found in Solar Energy 73, pg 311. Full ref: Perez
-    # et. al., Vol. 73, pp. 307-317 (2002). It is slightly different
-    # than the equation given in Solar Energy 73, pg 156. We used the
-    # equation from pg 311 because of the existence of known typos in
-    # the pg 156 publication (notably the fh2-(TL-1) should be fh2 *
-    # (TL-1)).
+    # ghi is calculated using either the equations in [1] by setting
+    # perez_enhancement=False (default behavior) or using the model
+    # in [2] by setting perez_enhancement=True.
 
     # The NaN handling is a little subtle. The AM input is likely to
     # have NaNs that we'll want to map to 0s in the output. However, we
@@ -119,8 +106,12 @@ def ineichen(apparent_zenith, airmass_absolute, linke_turbidity,
     cg1 = 5.09e-05 * altitude + 0.868
     cg2 = 3.92e-05 * altitude + 0.0387
 
-    ghi = (np.exp(-cg2*airmass_absolute*(fh1 + fh2*(tl - 1))) *
-           np.exp(0.01*airmass_absolute**1.8))
+    ghi = np.exp(-cg2*airmass_absolute*(fh1 + fh2*(tl - 1)))
+
+    # https://github.com/pvlib/pvlib-python/issues/435
+    if perez_enhancement:
+        ghi *= np.exp(0.01*airmass_absolute**1.8)
+
     # use fmax to map airmass nans to 0s. multiply and divide by tl to
     # reinsert tl nans
     ghi = cg1 * dni_extra * cos_zenith * tl / tl * np.fmax(ghi, 0)

pvlib/location.py

@@ -218,7 +218,7 @@ def get_clearsky(self, times, model='ineichen', solar_position=None,
 
             cs = clearsky.ineichen(apparent_zenith, airmass_absolute,
                                    linke_turbidity, altitude=self.altitude,
-                                   dni_extra=dni_extra)
+                                   dni_extra=dni_extra, **kwargs)
         elif model == 'haurwitz':
             cs = clearsky.haurwitz(apparent_zenith)
         elif model == 'simplified_solis':

pvlib/test/test_clearsky.py

@@ -3,6 +3,7 @@
 from collections import OrderedDict
 
 import numpy as np
+from numpy import nan
 import pandas as pd
 import pytz
 
@@ -20,35 +21,66 @@
 
 
 def test_ineichen_series():
-    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_localized, tus.latitude,
-                                                 tus.longitude)
-    am = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
-    am = atmosphere.absoluteairmass(am, atmosphere.alt2pres(tus.altitude))
+    times = pd.date_range(start='2014-06-24', end='2014-06-25', freq='3h',
+                          tz='America/Phoenix')
+    apparent_zenith = pd.Series(np.array(
+        [124.0390863, 113.38779941, 82.85457044, 46.0467599, 10.56413562,
+         34.86074109, 72.41687122, 105.69538659, 124.05614124]),
+        index=times)
+    am = pd.Series(np.array(
+        [nan, nan, 6.97935524, 1.32355476, 0.93527685,
+         1.12008114, 3.01614096, nan, nan]),
+        index=times)
     expected = pd.DataFrame(np.
-        array([[    0.        ,     0.        ,     0.        ],
-               [    0.        ,     0.        ,     0.        ],
-               [   91.12492792,   321.16092181,    51.17628184],
-               [  716.46580533,   888.90147035,    99.5050056 ],
-               [ 1053.42066043,   953.24925854,   116.32868969],
-               [  863.54692781,   922.06124712,   106.95536561],
-               [  271.06382274,   655.44925241,    73.05968071],
-               [    0.        ,     0.        ,     0.        ],
-               [    0.        ,     0.        ,     0.        ]]),
+        array([[   0.        ,    0.        ,    0.        ],
+               [   0.        ,    0.        ,    0.        ],
+               [  65.49426624,  321.16092181,   25.54562017],
+               [ 704.6968125 ,  888.90147035,   87.73601277],
+               [1044.1230677 ,  953.24925854,  107.03109696],
+               [ 853.02065704,  922.06124712,   96.42909484],
+               [ 251.99427693,  655.44925241,   53.9901349 ],
+               [   0.        ,    0.        ,    0.        ],
+               [   0.        ,    0.        ,    0.        ]]),
                             columns=['ghi', 'dni', 'dhi'],
-                            index=times_localized)
+                            index=times)
 
-    out = clearsky.ineichen(ephem_data['apparent_zenith'], am, 3)
+    out = clearsky.ineichen(apparent_zenith, am, 3)
+    assert_frame_equal(expected, out)
+
+
+def test_ineichen_series_perez_enhancement():
+    times = pd.date_range(start='2014-06-24', end='2014-06-25', freq='3h',
+                          tz='America/Phoenix')
+    apparent_zenith = pd.Series(np.array(
+        [124.0390863, 113.38779941, 82.85457044, 46.0467599, 10.56413562,
+         34.86074109, 72.41687122, 105.69538659, 124.05614124]),
+        index=times)
+    am = pd.Series(np.array(
+        [nan, nan, 6.97935524, 1.32355476, 0.93527685,
+         1.12008114, 3.01614096, nan, nan]),
+        index=times)
+    expected = pd.DataFrame(np.
+        array([[   0.        ,    0.        ,    0.        ],
+               [   0.        ,    0.        ,    0.        ],
+               [  91.1249279 ,  321.16092171,   51.17628184],
+               [ 716.46580547,  888.9014706 ,   99.50500553],
+               [1053.42066073,  953.24925905,  116.3286895 ],
+               [ 863.54692748,  922.06124652,  106.9553658 ],
+               [ 271.06382275,  655.44925213,   73.05968076],
+               [   0.        ,    0.        ,    0.        ],
+               [   0.        ,    0.        ,    0.        ]]),
+                            columns=['ghi', 'dni', 'dhi'],
+                            index=times)
+
+    out = clearsky.ineichen(apparent_zenith, am, 3, perez_enhancement=True)
     assert_frame_equal(expected, out)
 
 
 def test_ineichen_scalar_input():
     expected = OrderedDict()
-    expected['ghi'] = 1048.592893113678
+    expected['ghi'] = 1038.159219
     expected['dni'] = 942.2081860378344
-    expected['dhi'] = 120.6989665520498
+    expected['dhi'] = 110.26529293612793
 
     out = clearsky.ineichen(10., 1., 3.)
     for k, v in expected.items():
@@ -74,9 +106,9 @@ def test_ineichen_nans():
     expected['dni'] = np.full(length, np.nan)
     expected['dhi'] = np.full(length, np.nan)
 
-    expected['ghi'][length-1] = 1053.205472
-    expected['dni'][length-1] = 946.352797
-    expected['dhi'][length-1] = 121.2299
+    expected['ghi'][length-1] = 1042.72590228
+    expected['dni'][length-1] = 946.35279683
+    expected['dhi'][length-1] = 110.75033088
 
     out = clearsky.ineichen(apparent_zenith, airmass_absolute,
                             linke_turbidity, dni_extra=dni_extra)
@@ -89,43 +121,43 @@ def test_ineichen_arrays():
     expected = OrderedDict()
 
     expected['ghi'] = (np.
-        array([[[ 1106.78342709,  1064.7691287 ,  1024.34972343],
-                [  847.84529406,   815.66047425,   784.69741345],
-                [  192.19092519,   184.89521884,   177.87646277]],
+        array([[[1095.77074798, 1054.17449885, 1014.15727338],
+                [ 839.40909243,  807.54451692,  776.88954373],
+                [ 190.27859353,  183.05548067,  176.10656239]],
 
-               [[  959.12310134,   775.2374976 ,   626.60692548],
-                [  734.73092205,   593.86637713,   480.00875328],
-                [  166.54997871,   134.61857872,   108.80915072]],
+               [[ 773.49041181,  625.19479557,  505.33080493],
+                [ 592.52803177,  478.92699901,  387.10585505],
+                [ 134.31520045,  108.56393694,   87.74977339]],
 
-               [[ 1026.15144142,   696.85030591,   473.22483724],
-                [  786.0776095 ,   533.81830453,   362.51125692],
-                [  178.18932781,   121.00678573,    82.17463061]]]))
+               [[ 545.9968869 ,  370.78162375,  251.79449885],
+                [ 418.25788117,  284.03520249,  192.88577665],
+                [  94.81136442,   64.38555328,   43.72365587]]]))
 
     expected['dni'] = (np.
-        array([[[ 1024.58284359,   942.20818604,   861.11344424],
-                [ 1024.58284359,   942.20818604,   861.11344424],
-                [ 1024.58284359,   942.20818604,   861.11344424]],
+        array([[[1014.38807396,  942.20818604,  861.11344424],
+                [1014.38807396,  942.20818604,  861.11344424],
+                [1014.38807396,  942.20818604,  861.11344424]],
 
-               [[  687.61305142,   419.14891162,   255.50098235],
-                [  687.61305142,   419.14891162,   255.50098235],
-                [  687.61305142,   419.14891162,   255.50098235]],
+               [[ 687.61305142,  419.14891162,  255.50098235],
+                [ 687.61305142,  419.14891162,  255.50098235],
+                [ 687.61305142,  419.14891162,  255.50098235]],
 
-               [[  458.62196014,   186.46177428,    75.80970012],
-                [  458.62196014,   186.46177428,    75.80970012],
-                [  458.62196014,   186.46177428,    75.80970012]]]))
+               [[ 458.62196014,  186.46177428,   75.80970012],
+                [ 458.62196014,  186.46177428,   75.80970012],
+                [ 458.62196014,  186.46177428,   75.80970012]]]))
 
     expected['dhi'] = (np.
-            array([[[  82.20058349,  122.56094266,  163.23627919],
-                    [  62.96930021,   93.88712907,  125.04624459],
-                    [  14.27398153,   21.28248435,   28.34568241]],
+        array([[[ 81.38267402, 111.96631281, 153.04382915],
+                [ 62.3427452 ,  85.77117175, 117.23837487],
+                [ 14.13195304,  19.44274618,  26.57578203]],
 
-                   [[ 271.51004993,  356.08858598,  371.10594313],
-                    [ 207.988765  ,  272.77968255,  284.28364554],
-                    [  47.14722539,   61.83413404,   64.44187075]],
+               [[ 85.87736039, 206.04588395, 249.82982258],
+                [ 65.78587472, 157.84030442, 191.38074731],
+                [ 14.91244713,  35.77949226,  43.38249342]],
 
-                   [[ 567.52948128,  510.38853163,  397.41513712],
-                    [ 434.75280544,  390.98029849,  304.4376574 ],
-                    [  98.5504602 ,   88.62803842,   69.01041434]]]))
+               [[ 87.37492676, 184.31984947, 175.98479873],
+                [ 66.93307711, 141.19719644, 134.81217714],
+                [ 15.17249681,  32.00680597,  30.5594396 ]]]))
 
     apparent_zenith = np.linspace(0, 80, 3)
     airmass_absolute = np.linspace(1, 10, 3)
@@ -142,7 +174,7 @@ def test_ineichen_arrays():
 
 def test_ineichen_dni_extra():
     expected = pd.DataFrame(
-        np.array([[ 1053.20547182,   946.35279683,   121.22990042]]),
+        np.array([[1042.72590228,  946.35279683,  110.75033088]]),
         columns=['ghi', 'dni', 'dhi'])
 
     out = clearsky.ineichen(10, 1, 3, dni_extra=pd.Series(1370))
@@ -151,7 +183,7 @@ def test_ineichen_dni_extra():
 
 def test_ineichen_altitude():
     expected = pd.DataFrame(
-        np.array([[ 1145.64245696,   994.95377835,   165.80426215]]),
+        np.array([[1134.24312405,  994.95377835,  154.40492924]]),
         columns=['ghi', 'dni', 'dhi'])
 
     out = clearsky.ineichen(10, 1, 3, altitude=pd.Series(2000))
@@ -526,18 +558,18 @@ def test_detect_clearsky_components(detect_clearsky_data):
     assert_series_equal(expected['Clear or not'], clear_samples,
                         check_dtype=False, check_names=False)
     assert isinstance(components, OrderedDict)
-    assert np.allclose(alpha, 0.95345573579557108)
+    assert np.allclose(alpha, 0.9633903181941296)
 
 
 @requires_scipy
 def test_detect_clearsky_iterations(detect_clearsky_data):
     expected, cs = detect_clearsky_data
-    alpha = 1.0348
+    alpha = 1.0448
     with pytest.warns(RuntimeWarning):
         clear_samples = clearsky.detect_clearsky(
             expected['GHI'], cs['ghi']*alpha, cs.index, 10, max_iterations=1)
-    assert (clear_samples[:'2012-04-01 10:39:00'] == True).all()
-    assert (clear_samples['2012-04-01 10:40:00':] == False).all()
+    assert (clear_samples[:'2012-04-01 10:41:00'] == True).all()
+    assert (clear_samples['2012-04-01 10:42:00':] == False).all()
     clear_samples = clearsky.detect_clearsky(
             expected['GHI'], cs['ghi']*alpha, cs.index, 10, max_iterations=20)
     assert_series_equal(expected['Clear or not'], clear_samples,