note and tests for disc relative vs absolute airmass (#625)

Author: wholmgren

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

@@ -44,6 +44,8 @@ Enhancements
   (:issue:`604`)
 * Change :py:func:`pvlib.pvsystem.sapm_spectral_loss` to avoid numpy warning.
 * Add warning message when :py:func:`pvlib.spa` is reloaded.
+* Add option for :py:func:`pvlib.irradiance.disc` to use relative airmass
+  by supplying `pressure=None`. (:issue:`449`)
 
 
 Bug fixes
@@ -69,3 +71,4 @@ Contributors
 * Ben Ellis (:ghuser:`bhellis725`)
 * Cliff Hansen (:ghuser:`cwhanse`)
 * Mark Mikofski (:ghuser:`mikofski`)
+* Anton Driesse (:ghuser:`adriesse`)

pvlib/irradiance.py

@@ -1306,6 +1306,13 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
 
     The pvlib implementation limits the clearness index to 1.
 
+    The original report describing the DISC model [1]_ uses the
+    relative airmass rather than the absolute (pressure-corrected)
+    airmass. However, the NREL implementation of the DISC model [2]_
+    uses absolute airmass. PVLib Matlab also uses the absolute airmass.
+    pvlib python defaults to absolute airmass, but the relative airmass
+    can be used by supplying `pressure=None`.
+
     Parameters
     ----------
     ghi : numeric
@@ -1319,8 +1326,9 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
         Day of year or array of days of year e.g.
         pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
 
-    pressure : numeric, default 101325
-        Site pressure in Pascal.
+    pressure : None or numeric, default 101325
+        Site pressure in Pascal. If None, relative airmass is used
+        instead of absolute (pressure-corrected) airmass.
 
     min_cos_zenith : numeric, default 0.065
         Minimum value of cos(zenith) to allow when calculating global
@@ -1344,15 +1352,13 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
 
     References
     ----------
-    [1] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
-    Global Horizontal to Direct Normal Insolation", Technical
-    Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research
-    Institute, 1987.
+    .. [1] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
+       Global Horizontal to Direct Normal Insolation", Technical
+       Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research
+       Institute, 1987.
 
-    [2] J.W. "Fourier series representation of the position of the sun".
-    Found at:
-    http://www.mail-archive.com/sundial@uni-koeln.de/msg01050.html on
-    January 12, 2012
+    .. [2] Maxwell, E. "DISC Model", Excel Worksheet.
+       https://www.nrel.gov/grid/solar-resource/disc.html
 
     See Also
     --------
@@ -1367,7 +1373,8 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
                          max_clearness_index=1)
 
     am = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
-    am = atmosphere.get_absolute_airmass(am, pressure)
+    if pressure is not None:
+        am = atmosphere.get_absolute_airmass(am, pressure)
 
     Kn = _disc_kn(kt, am)
     dni = Kn * I0

pvlib/test/test_irradiance.py

@@ -359,17 +359,23 @@ def test_poa_components(irrad_data, ephem_data, dni_et, relative_airmass):
     assert_frame_equal(out, expected)
 
 
-def test_disc_value():
+@pytest.mark.parametrize('pressure,expected', [
+    (93193,  [[830.46567,   0.79742,   0.93505],
+              [676.09497,   0.63776,   3.02102]]),
+    (None,   [[868.72425,   0.79742,   1.01664],
+              [680.66679,   0.63776,   3.28463]]),
+    (101325, [[868.72425,   0.79742,   1.01664],
+              [680.66679,   0.63776,   3.28463]])
+])
+def test_disc_value(pressure, expected):
+    # see GH 449 for pressure=None vs. 101325.
     columns = ['dni', 'kt', 'airmass']
     times = pd.DatetimeIndex(['2014-06-24T1200', '2014-06-24T1800'],
                              tz='America/Phoenix')
     ghi = pd.Series([1038.62, 254.53], index=times)
     zenith = pd.Series([10.567, 72.469], index=times)
-    pressure = 93193.
     out = irradiance.disc(ghi, zenith, times, pressure=pressure)
-    expected_values = np.array(
-        [[830.46567,   0.79742,   0.93505],
-         [676.09497,   0.63776,   3.02102]])
+    expected_values = np.array(expected)
     expected = pd.DataFrame(expected_values, columns=columns, index=times)
     # check the pandas dataframe. check_less_precise is weird
     assert_frame_equal(out, expected, check_less_precise=True)