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

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

@@ -14,7 +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
 ~~~~~~~~~~~~
@@ -50,3 +53,4 @@ Contributors
 * Eesh Saxena (:ghuser:`eeshsaxena`)
 * Karl Hill (:ghuser:`karlhillx`)
 * Yonry Zhu (:ghuser:`yonryzhu`)
+* Mark Campanelli (:ghuser:`markcampanelli`)

pvlib/irradiance.py

@@ -1042,10 +1042,15 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     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
     -------
@@ -1125,6 +1130,11 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
             * poa_circumsolar
             * poa_horizon
 
+    Raises
+    ------
+    FloatingPointError
+        If dni is zero when dhi is not zero and not NaN.
+
 
     References
     ----------
@@ -1147,12 +1157,21 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     kappa = 1.041  # for solar_zenith in radians
     z = np.radians(solar_zenith)  # convert to radians
 
-    # delta is the sky's "brightness"
+    # delta is the sky's "brightness", NaN airmass ok, assumes dni_extra > 0
     delta = dhi * airmass / dni_extra
 
-    # epsilon is the sky's "clearness"
-    with np.errstate(invalid='ignore'):
-        eps = ((dhi + dni) / dhi + kappa * (z ** 3)) / (1 + kappa * (z ** 3))
+    # epsilon is the sky's "clearness". Preserves NaNs for dni or dhi.
+    # Assumes:
+    # - dni >=0 and dhi >= 0.
+    # - dni->0^+ faster than dhi->0^+.
+    irr_ratio = np.zeros(np.broadcast(dni, dhi).shape)
+    with np.errstate(divide="raise"):
+        irr_ratio = np.divide(
+            dni, dhi, out=irr_ratio, where=np.logical_not(
+                np.logical_and(dhi == 0, dni == 0)
+            )
+        )
+    eps = 1 + irr_ratio / (1 + kappa * (z ** 3))
 
     # numpy indexing below will not work with a Series
     if isinstance(eps, pd.Series):
@@ -1192,17 +1211,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
 
     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(sky_diffuse, pd.Series):
-        sky_diffuse[np.isnan(airmass)] = 0
+        sky_diffuse[airmass_nan_idx] = 0
     else:
-        sky_diffuse = np.where(np.isnan(airmass), 0, sky_diffuse)
+        sky_diffuse = np.where(airmass_nan_idx, 0, sky_diffuse)
 
     if return_components:
         diffuse_components = OrderedDict()

tests/test_irradiance.py

@@ -7,8 +7,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 (
@@ -287,7 +288,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)
@@ -393,6 +393,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