DRAFT: Handle dni/dhi=0/0 in irradiance.perez without NaNs

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

@@ -14,6 +14,10 @@ Deprecations
 
 Bug fixes
 ~~~~~~~~~
+* Fix :py:func:`pvlib.irradiance.perez` to return 0 when both DHI and
+  DNI are zero. (:issue:`2801`, :pull:`2808`)
+* Fix :py:func:`pvlib.irradiance.perez` to return NaN when
+  airmass is NaN and either DHI or DNI are NaN. (:issue:`2801`, :pull:`2808`)
 
 
 Enhancements
@@ -43,3 +47,4 @@ Maintenance
 Contributors
 ~~~~~~~~~~~~
 
+* Mark Campanelli (:ghuser:`markcampanelli`)

pvlib/irradiance.py

@@ -1156,20 +1156,25 @@ def king(surface_tilt, dhi, ghi, solar_zenith):
     return poa_sky_diffuse
 
 
-def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
-          solar_zenith, solar_azimuth, airmass,
-          model='allsitescomposite1990', return_components=False):
+def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, solar_zenith,
+          solar_azimuth, airmass, model='allsitescomposite1990',
+          return_components=False):
     '''
     Determine diffuse irradiance from the sky on a tilted surface using
     one of the Perez models.
 
     Perez models determine the diffuse irradiance from the sky (ground
     reflected irradiance is not included in this algorithm) on a tilted
     surface using the surface tilt angle, surface azimuth angle, diffuse
-    horizontal irradiance, direct normal irradiance, extraterrestrial
-    irradiance, sun zenith angle, sun azimuth angle, and relative (not
-    pressure-corrected) airmass. Optionally a selector may be used to
-    use any of Perez's model coefficient sets.
+    horizontal irradiance (DHI), direct normal irradiance (DNI),
+    extraterrestrial irradiance, sun zenith angle, sun azimuth angle, and
+    relative (not pressure-corrected) airmass. Optionally a selector may be
+    used to use any of Perez's model coefficient sets. It is expected that if
+    DHI is zero, then DNI is also zero, otherwise a FloatingPointError is
+    raised due to a division of a nonzero (and not NaN) value by zero. It is
+    also expected that extraterrestrial irradiance is positive. If airmass
+    is NaN, then the total and all components are zero if they should not
+    otherwise be NaN because DHI or DNI was NaN.
 
     Warning
     -------
@@ -1202,11 +1207,13 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     solar_azimuth : numeric
         Solar azimuth angle. See :term:`solar_azimuth`. [°]
 
-    airmass : numeric
-        Relative (not pressure-corrected) airmass values. If AM is a
-        DataFrame it must be of the same size as all other DataFrame
+    airmass : numeric or None
+        Relative (not pressure-corrected) airmass values. If ``airmass`` is a
+        DataFrame, then it must be of the same size as all other DataFrame
         inputs. AM must be >=0 (careful using the 1/cos(z) model of AM
-        generation). [unitless]
+        generation). If None, then the 'kastenyoung1989' model is used, which
+        corresponds to the calibration of the F1 and F2 coefficients in the
+        original work by Perez. [unitless]
 
     model : string, default 'allsitescomposite1990'
         A string which selects the desired set of Perez coefficients. If
@@ -1252,6 +1259,11 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
             * poa_horizon: The portion of sky diffuse irradiance on a tilted
               plane from the horizon. [Wm⁻²]
 
+    Raises
+    ------
+    FloatingPointError
+        If dni is zero when dhi is not zero and not NaN.
+
 
     References
     ----------
@@ -1271,31 +1283,20 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
        Perez Diffuse Radiation Model". SAND88-7030
     '''
 
-    kappa = 1.041  # for solar_zenith in radians
-    z = np.radians(solar_zenith)  # convert to radians
+    # FIXME Verify _calc_zeta assumes:
+    # - dni >=0 and dhi >= 0.
+    # - dni->0^+ faster than dhi->0^+.
 
-    # delta is the sky's "brightness"
-    delta = dhi * airmass / dni_extra
-
-    # epsilon is the sky's "clearness"
-    # np.divide so a Python-scalar dhi=0 yields inf/nan like the array
-    # path (handled below via digitize) instead of raising ZeroDivisionError,
-    # which np.errstate cannot suppress for native scalar division.
-    with np.errstate(invalid='ignore', divide='ignore'):
-        eps = (np.divide(dhi + dni, dhi) + kappa * (z ** 3)) \
-            / (1 + kappa * (z ** 3))
-
-    # numpy indexing below will not work with a Series
-    if isinstance(eps, pd.Series):
-        eps = eps.values
+    delta = _calc_delta(dhi, dni_extra, solar_zenith, airmass)
+    zeta, z = _calc_zeta(dhi, dni, solar_zenith)
 
     # Perez et al define clearness bins according to the following
     # rules. 1 = overcast ... 8 = clear (these names really only make
     # sense for small zenith angles, but...) these values will
     # eventually be used as indicies for coeffecient look ups
-    ebin = np.digitize(eps, (0., 1.065, 1.23, 1.5, 1.95, 2.8, 4.5, 6.2))
+    ebin = np.digitize(zeta, (0., 1.065, 1.23, 1.5, 1.95, 2.8, 4.5, 6.2))
     ebin = np.array(ebin)  # GH 642
-    ebin[np.isnan(eps)] = 0
+    ebin[np.isnan(zeta)] = 0
 
     # correct for 0 indexing in coeffecient lookup
     # later, ebin = -1 will yield nan coefficients
@@ -1305,7 +1306,7 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     # in a subfunction to clean up the code.
     F1c, F2c = _get_perez_coefficients(model)
 
-    # results in invalid eps (ebin = -1) being mapped to nans
+    # results in invalid zeta (ebin = -1) being mapped to nans
     nans = np.array([np.nan, np.nan, np.nan])
     F1c = np.vstack((F1c, nans))
     F2c = np.vstack((F2c, nans))
@@ -1323,17 +1324,23 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     B = np.maximum(B, tools.cosd(85))
 
     # Calculate Diffuse POA from sky dome
-    term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))
-    term2 = F1 * A / B
-    term3 = F2 * tools.sind(surface_tilt)
+    term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))  # isotropic
+    term2 = F1 * A / B  # circumsolar
+    term3 = F2 * tools.sind(surface_tilt)  # horizon
 
     poa_sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0)
 
+    # Use NaN airmass to coerce to zero values that are not otherwise NaN.
+    airmass_nan_idx = np.logical_and(
+        np.isnan(airmass), np.logical_not(
+            np.logical_or(np.isnan(dhi), np.isnan(dni))
+        )
+    )
     # we've preserved the input type until now, so don't ruin it!
     if isinstance(poa_sky_diffuse, pd.Series):
-        poa_sky_diffuse[np.isnan(airmass)] = 0
+        poa_sky_diffuse[airmass_nan_idx] = 0
     else:
-        poa_sky_diffuse = np.where(np.isnan(airmass), 0, poa_sky_diffuse)
+        poa_sky_diffuse = np.where(airmass_nan_idx, 0, poa_sky_diffuse)
 
     if return_components:
         diffuse_components = {
@@ -1361,12 +1368,13 @@ def _calc_delta(dhi, dni_extra, solar_zenith, airmass=None):
     Compute the delta parameter, which represents sky dome "brightness"
     in the Perez and Perez-Driesse models.
 
-    Helper function for perez_driesse transposition.
+    Helper function for perez and perez_driesse transposition.
     '''
     if airmass is None:
         # use the same airmass model as in the original perez work
-        airmass = atmosphere.get_relative_airmass(solar_zenith,
-                                                  'kastenyoung1989')
+        airmass = atmosphere.get_relative_airmass(
+            solar_zenith, 'kastenyoung1989'
+        )
 
     max_airmass = atmosphere.get_relative_airmass(90, 'kastenyoung1989')
     airmass = np.where(solar_zenith >= 90, max_airmass, airmass)
@@ -1377,9 +1385,9 @@ def _calc_delta(dhi, dni_extra, solar_zenith, airmass=None):
 def _calc_zeta(dhi, dni, zenith):
     '''
     Compute the zeta parameter, which represents sky dome "clearness"
-    in the Perez-Driesse model.
+    in the Perez and Perez-Driesse model.
 
-    Helper function for perez_driesse transposition.
+    Helper function for perez and perez_driesse transposition.
     '''
     dhi = np.asarray(dhi)
     dni = np.asarray(dni)
@@ -1393,10 +1401,10 @@ def _calc_zeta(dhi, dni, zenith):
 
     # then apply the kappa correction in a manner analogous to eq. 7
     kappa = 1.041
-    kterm = kappa * np.radians(zenith) ** 3
-    zeta = zeta / (1 - kterm * (zeta - 1))
+    z = np.radians(zenith)
+    zeta = zeta / (1 - kappa * z ** 3 * (zeta - 1))
 
-    return zeta
+    return zeta, z
 
 
 def _f(i, j, zeta):
@@ -1470,9 +1478,11 @@ def perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
 
     airmass : numeric, optional
         Relative (not pressure-corrected) airmass values. If ``airmass`` is a
-        DataFrame it must be of the same size as all other DataFrame
+        DataFrame, then it must be of the same size as all other DataFrame
         inputs. AM must be >=0 (careful using the 1/cos(z) model of AM
-        generation). [unitless]
+        generation). If not provided or None, then the kastenyoung1989 model
+        is used, which corresponds to the calibration of the F1 and F2
+        coefficients in the original work by Perez. [unitless]
 
     return_components: bool (optional, default=False)
         Flag used to decide whether to return the calculated diffuse components
@@ -1529,9 +1539,7 @@ def perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct., 2023
 
     delta = _calc_delta(dhi, dni_extra, solar_zenith, airmass)
-    zeta = _calc_zeta(dhi, dni, solar_zenith)
-
-    z = np.radians(solar_zenith)
+    zeta, z = _calc_zeta(dhi, dni, solar_zenith)
 
     F1 = _f(1, 1, zeta) + _f(1, 2, zeta) * delta + _f(1, 3, zeta) * z
     F2 = _f(2, 1, zeta) + _f(2, 2, zeta) * delta + _f(2, 3, zeta) * z

tests/test_irradiance.py

@@ -6,8 +6,9 @@
 import pandas as pd
 
 import pytest
-from numpy.testing import (assert_almost_equal,
-                           assert_allclose)
+from numpy.testing import (
+    assert_almost_equal, assert_allclose, assert_equal, assert_raises
+)
 from pvlib import irradiance, albedo
 
 from .conftest import (
@@ -350,7 +351,6 @@ def test_perez_driesse_airmass(irrad_data, ephem_data, dni_et):
     out = irradiance.perez_driesse(40, 180, irrad_data['dhi'], dni,
                                    dni_et, ephem_data['apparent_zenith'],
                                    ephem_data['azimuth'], airmass=None)
-    print(out)
     expected = pd.Series(np.array(
         [0.,   29.991,  np.nan,   47.397]),
         index=irrad_data.index)
@@ -456,6 +456,133 @@ def test_perez_negative_horizon():
     assert_series_equal(sum_components, expected_for_sum, check_less_precise=2)
 
 
+def test_perez_zero_dhi_and_dni_scalar():
+    # Divides zero by zero.
+    args = (20, 180, 0.0, 0.0, 1366.1, 89.96, 256.28, 37.32)
+
+    out = irradiance.perez(*args)
+    poa_sky_diffuse_expected = 0.0
+    assert_equal(out, poa_sky_diffuse_expected)
+
+    out = irradiance.perez(*args, return_components=True)
+    expected = {
+        "poa_sky_diffuse": poa_sky_diffuse_expected,
+        "poa_isotropic": 0.0,
+        "poa_circumsolar": 0.0,
+        "poa_horizon": 0.0,
+    }
+    assert len(out) == len(expected)
+    for key in expected.keys():
+        assert_equal(out[key], expected[key], err_msg=key)
+
+    assert_equal(
+        out["poa_sky_diffuse"],
+        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
+    )
+
+
+def test_perez_array_dhi_and_dni_combos():
+    # Divides zero and non-zero by zero and various NaN division combos.
+    args = (
+        20, 180,
+        np.array([0.0, 10.0, np.nan, np.nan, 0.0, 100.0, np.nan]),
+        np.array([0.0, 0.0, 0.0, 100.0, np.nan, np.nan, np.nan]),
+        1366.1, 89.96, 256.28, 37.32
+    )
+
+    out = irradiance.perez(*args)
+    poa_sky_diffuse_expected = np.array(
+        [0.0, 9.424924186619206, np.nan, np.nan, np.nan, np.nan, np.nan]
+    )
+    assert_allclose(out, poa_sky_diffuse_expected)
+
+    out = irradiance.perez(*args, return_components=True)
+    expected = {
+        "poa_sky_diffuse": poa_sky_diffuse_expected,
+        "poa_isotropic": np.array(
+            [0.0, 9.162258932459126, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+        "poa_circumsolar": np.array(
+            [0.0, 0.5187450944545264, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+        "poa_horizon": np.array(
+            [0.0, -0.2560798402944465, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+    }
+    assert len(out) == len(expected)
+    for key in expected.keys():
+        assert_allclose(out[key], expected[key], err_msg=key)
+
+    assert_almost_equal(
+        out["poa_sky_diffuse"],
+        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
+    )
+
+
+def test_perez_array_dhi_and_dni_combos_nan_airmass():
+    # Divides zero and non-zero by zero and various NaN division combos, when
+    # airmass is NaN (e.g., for sun below horizon).
+    args = (
+        20, 180,
+        np.array([0.0, 10.0, np.nan, np.nan, 0.0, 100.0, np.nan]),
+        np.array([0.0, 0.0, 0.0, 100.0, np.nan, np.nan, np.nan]),
+        1366.1, 91, 256.28, np.nan
+    )
+
+    out = irradiance.perez(*args)
+    poa_sky_diffuse_expected = np.array(
+        [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
+    )
+    assert_allclose(out, poa_sky_diffuse_expected)
+
+    out = irradiance.perez(*args, return_components=True)
+    expected = {
+        "poa_sky_diffuse": poa_sky_diffuse_expected,
+        "poa_isotropic": np.array(
+            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+        "poa_circumsolar": np.array(
+            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+        "poa_horizon": np.array(
+            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+    }
+    assert len(out) == len(expected)
+    for key in expected.keys():
+        assert_allclose(out[key], expected[key], err_msg=key)
+
+    assert_almost_equal(
+        out["poa_sky_diffuse"],
+        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
+    )
+
+
+def test_perez_zero_dhi_nonzero_dni_scalar():
+    # Divides nonzero by zero.
+    args = (20, 180, 0.0, 100.0, 1366.1, 89.96, 256.28, 37.32)
+
+    with assert_raises(FloatingPointError):
+        irradiance.perez(*args)
+
+    with assert_raises(FloatingPointError):
+        irradiance.perez(*args, return_components=True)
+
+
+def test_perez_zero_dhi_nonzero_dni_array():
+    # Divides nonzero by zero.
+    args = (
+        20, 180, np.array([0.0, 10.0, np.nan]), 100.0, 1366.1, 89.96,
+        256.28, 37.32
+    )
+
+    with assert_raises(FloatingPointError):
+        irradiance.perez(*args)
+
+    with assert_raises(FloatingPointError):
+        irradiance.perez(*args, return_components=True)
+
+
 def test_perez_arrays(irrad_data, ephem_data, dni_et, relative_airmass):
     dni = irrad_data['dni'].copy()
     dni.iloc[2] = np.nan