fix erbs near zenith=90 (#683)

* fix erbs near zenith=90

* document return type change

* fix copy/paste in test

* fix docs

* revert accidental changes to test_clearsky

Author: wholmgren

docs/sphinx/source/whatsnew.rst

@@ -6,6 +6,7 @@ What's New
 
 These are new features and improvements of note in each release.
 
+.. include:: whatsnew/v0.6.2.rst
 .. include:: whatsnew/v0.6.1.rst
 .. include:: whatsnew/v0.6.0.rst
 .. include:: whatsnew/v0.5.2.rst

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

@@ -14,20 +14,31 @@ date will require Python 3. (:issue:`501`)
 
 API Changes
 ~~~~~~~~~~~
-
+* :py:func:`~pvlib.irradiance.erbs` *doy* argument changed to
+  *datetime_or_doy* to be consistent with allowed types and similar
+  functions (:py:func:`~pvlib.irradiance.disc`,
+  :py:func:`~pvlib.irradiance.get_extra_radiation`).
+* :py:func:`~pvlib.irradiance.erbs` DataFrame vs. OrderedDict return
+  behavior now determined by type of *datetime_or_doy* instead of
+  *ghi* or *zenith*.
+* Added *min_cos_zenith* and *max_zenith* keyword arguments to
+  :py:func:`~pvlib.irradiance.erbs`. (:issue:`681`)
 
 Enhancements
 ~~~~~~~~~~~~
-* Add US CRN data reader to `pvlib.iotools`.
-* Add SOLRAD data reader to `pvlib.iotools`.
-* Add EPW data reader to `pvlib.iotools`. (:issue:`591`)
+* Add US CRN data reader to :ref:`iotools`.
+* Add SOLRAD data reader to :ref:`iotools`.
+* Add EPW data reader to :ref:`iotools`. (:issue:`591`)
 
 Bug fixes
 ~~~~~~~~~
 * Compatibility with pandas 0.24 deprecations. (:issue:`659`)
-* pvwatts_ac raised ZeroDivisionError when called with scalar pdc=0
-  and a RuntimeWarning for array(0) input. Now correctly returns 0s
-  of the appropriate type. (:issue:`675`)
+* :py:func:`~pvlib.pvsystem.pvwatts_ac` raised ``ZeroDivisionError``
+  when called with scalar ``pdc=0``
+  and a ``RuntimeWarning`` for ``array(0)`` input. Now correctly returns
+  0s of the appropriate type. (:issue:`675`)
+* Fixed :py:func:`~pvlib.irradiance.erbs` behavior when zenith is
+  near 90 degrees. (:issue:`681`)
 
 
 Testing

pvlib/iotools/epw.py

@@ -70,8 +70,7 @@ def read_epw(filename, coerce_year=None):
     longitude         Float   site longitude
     TZ                Float   UTC offset
     altitude          Float   site elevation
-    ===============   ======  ========================================
-
+    ===============   ======  =========================================
 
 
     =============================       ==============================================================================================================================================================

pvlib/irradiance.py

@@ -2152,7 +2152,7 @@ def _gti_dirint_gte_90_kt_prime(aoi, solar_zenith, solar_azimuth, times,
     return kt_prime_gte_90
 
 
-def erbs(ghi, zenith, doy):
+def erbs(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     r"""
     Estimate DNI and DHI from GHI using the Erbs model.
 
@@ -2179,8 +2179,15 @@ def erbs(ghi, zenith, doy):
         Global horizontal irradiance in W/m^2.
     zenith: numeric
         True (not refraction-corrected) zenith angles in decimal degrees.
-    doy: scalar, array or DatetimeIndex
-        The day of the year.
+    datetime_or_doy : int, float, array, pd.DatetimeIndex
+        Day of year or array of days of year e.g.
+        pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
+    min_cos_zenith : numeric, default 0.065
+        Minimum value of cos(zenith) to allow when calculating global
+        clearness index `kt`. Equivalent to zenith = 86.273 degrees.
+    max_zenith : numeric, default 87
+        Maximum value of zenith to allow in DNI calculation. DNI will be
+        set to 0 for times with zenith values greater than `max_zenith`.
 
     Returns
     -------
@@ -2205,14 +2212,10 @@ def erbs(ghi, zenith, doy):
     disc
     """
 
-    dni_extra = get_extra_radiation(doy)
+    dni_extra = get_extra_radiation(datetime_or_doy)
 
-    # This Z needs to be the true Zenith angle, not apparent,
-    # to get extraterrestrial horizontal radiation)
-    i0_h = dni_extra * tools.cosd(zenith)
-
-    kt = ghi / i0_h
-    kt = np.maximum(kt, 0)
+    kt = clearness_index(ghi, zenith, dni_extra, min_cos_zenith=min_cos_zenith,
+                         max_clearness_index=1)
 
     # For Kt <= 0.22, set the diffuse fraction
     df = 1 - 0.09*kt
@@ -2229,14 +2232,18 @@ def erbs(ghi, zenith, doy):
     dhi = df * ghi
 
     dni = (ghi - dhi) / tools.cosd(zenith)
+    bad_values = (zenith > max_zenith) | (ghi < 0) | (dni < 0)
+    dni = np.where(bad_values, 0, dni)
+    # ensure that closure relationship remains valid
+    dhi = np.where(bad_values, ghi, dhi)
 
     data = OrderedDict()
     data['dni'] = dni
     data['dhi'] = dhi
     data['kt'] = kt
 
-    if isinstance(dni, pd.Series):
-        data = pd.DataFrame(data)
+    if isinstance(datetime_or_doy, pd.DatetimeIndex):
+        data = pd.DataFrame(data, index=datetime_or_doy)
 
     return data
 

pvlib/test/test_irradiance.py

@@ -620,21 +620,59 @@ def test_gti_dirint():
 
 
 def test_erbs():
-    ghi = pd.Series([0, 50, 1000, 1000])
-    zenith = pd.Series([120, 85, 10, 10])
-    doy = pd.Series([1, 1, 1, 180])
-    expected = pd.DataFrame(np.
-        array([[ -0.00000000e+00,   0.00000000e+00,  -0.00000000e+00],
-               [  9.67127061e+01,   4.15709323e+01,   4.05715990e-01],
-               [  7.94187742e+02,   2.17877755e+02,   7.18119416e-01],
-               [  8.42358014e+02,   1.70439297e+02,   7.68919470e-01]]),
-        columns=['dni', 'dhi', 'kt'])
+    index = pd.DatetimeIndex(['20190101']*3 + ['20190620'])
+    ghi = pd.Series([0, 50, 1000, 1000], index=index)
+    zenith = pd.Series([120, 85, 10, 10], index=index)
+    expected = pd.DataFrame(np.array(
+        [[0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
+         [9.67192672e+01, 4.15703604e+01, 4.05723511e-01],
+         [7.94205651e+02, 2.17860117e+02, 7.18132729e-01],
+         [8.42001578e+02, 1.70790318e+02, 7.68214312e-01]]),
+        columns=['dni', 'dhi', 'kt'], index=index)
 
-    out = irradiance.erbs(ghi, zenith, doy)
+    out = irradiance.erbs(ghi, zenith, index)
 
     assert_frame_equal(np.round(out, 0), np.round(expected, 0))
 
 
+def test_erbs_min_cos_zenith_max_zenith():
+    # map out behavior under difficult conditions with various
+    # limiting kwargs settings
+    columns = ['dni', 'dhi', 'kt']
+    times = pd.DatetimeIndex(['2016-07-19 06:11:00'], tz='America/Phoenix')
+
+    # max_zenith keeps these results reasonable
+    out = irradiance.erbs(ghi=1.0, zenith=89.99999,
+                          datetime_or_doy=times, min_cos_zenith=0)
+    expected = pd.DataFrame(np.array(
+        [[0., 1., 1.]]),
+        columns=columns, index=times)
+    assert_frame_equal(out, expected)
+
+    # 4-5 9s will produce bad behavior without max_zenith limit
+    out = irradiance.erbs(ghi=1.0, zenith=89.99999,
+                          datetime_or_doy=times, max_zenith=100)
+    expected = pd.DataFrame(np.array(
+        [[6.00115286e+03, 9.98952601e-01, 1.16377640e-02]]),
+        columns=columns, index=times)
+    assert_frame_equal(out, expected)
+
+    # 1-2 9s will produce bad behavior without either limit
+    out = irradiance.erbs(ghi=1.0, zenith=89.99, datetime_or_doy=times,
+                          min_cos_zenith=0, max_zenith=100)
+    expected = pd.DataFrame(np.array(
+        [[4.78419761e+03, 1.65000000e-01, 1.00000000e+00]]),
+        columns=columns, index=times)
+    assert_frame_equal(out, expected)
+
+    # check default behavior under hardest condition
+    out = irradiance.erbs(ghi=1.0, zenith=90, datetime_or_doy=times)
+    expected = pd.DataFrame(np.array(
+        [[0., 1., 0.01163776]]),
+        columns=columns, index=times)
+    assert_frame_equal(out, expected)
+
+
 def test_erbs_all_scalar():
     ghi = 1000
     zenith = 10