diff --git a/docs/sphinx/source/whatsnew/v0.16.2.rst b/docs/sphinx/source/whatsnew/v0.16.2.rst index f704525a12..9f58e5c971 100644 --- a/docs/sphinx/source/whatsnew/v0.16.2.rst +++ b/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`) diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 38bf122b91..55f9b21bc1 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -1156,9 +1156,9 @@ 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. @@ -1166,10 +1166,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 ------- @@ -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 diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py index b0d9863ad5..ae1c8c14e0 100644 --- a/tests/test_irradiance.py +++ b/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