diff --git a/docs/sphinx/source/whatsnew/v0.15.3.rst b/docs/sphinx/source/whatsnew/v0.15.3.rst index 7124d28cae..7b1d2dbb15 100644 --- a/docs/sphinx/source/whatsnew/v0.15.3.rst +++ b/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`) diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 9db6da3aa8..23307f9328 100644 --- a/pvlib/irradiance.py +++ b/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() diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py index acc8495d3b..94057b0ed7 100644 --- a/tests/test_irradiance.py +++ b/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