diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst index b13165335f..92c2d48589 100644 --- a/docs/sphinx/source/whatsnew/v0.6.1.rst +++ b/docs/sphinx/source/whatsnew/v0.6.1.rst @@ -21,6 +21,14 @@ API Changes * Added keyword argument ``horizon`` to :func:`~pvlib.solarposition.pyephem` and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'``. (:issue:`588`) +* Add `max_airmass` keyword argument to :py:func:`pvlib.irradiance.disc`. + Default value (`max_airmass=12`) is consistent with polynomial fit in + original paper describing the model. This change may result in different + output of functions that use the `disc` *Kn* calculation for times when + input zenith angles approach 90 degrees. This includes + :py:func:`pvlib.irradiance.dirint` and :py:func:`pvlib.irradiance.dirindex` + when `min_cos_zenith` and `max_zenith` kwargs are used, as well as + :py:func:`pvlib.irradiance.gti_dirint`. (:issue:`450`) * Changed key names for `components` returned from :py:func:`pvlib.clearsky.detect_clearsky`. (:issue:`596`) * Changed function name from `pvlib.solarposition.get_rise_set_transit` diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 524dd106bf..4855187ea4 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -1295,7 +1295,7 @@ def _kt_kt_prime_factor(airmass): def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325, - min_cos_zenith=0.065, max_zenith=87): + min_cos_zenith=0.065, max_zenith=87, max_airmass=12): """ Estimate Direct Normal Irradiance from Global Horizontal Irradiance using the DISC model. @@ -1338,6 +1338,11 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325, Maximum value of zenith to allow in DNI calculation. DNI will be set to 0 for times with zenith values greater than `max_zenith`. + max_airmass : numeric, default 12 + Maximum value of the airmass to allow in Kn calculation. + Default value (12) comes from range over which Kn was fit + to airmass in the original paper. + Returns ------- output : OrderedDict or DataFrame @@ -1376,7 +1381,7 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325, if pressure is not None: am = atmosphere.get_absolute_airmass(am, pressure) - Kn = _disc_kn(kt, am) + Kn, am = _disc_kn(kt, am, max_airmass=max_airmass) dni = Kn * I0 bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0) @@ -1393,16 +1398,29 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325, return output -def _disc_kn(clearness_index, airmass): +def _disc_kn(clearness_index, airmass, max_airmass=12): """ Calculate Kn for `disc` + + Parameters + ---------- + clearness_index : numeric + airmass : numeric + max_airmass : float + airmass > max_airmass is set to max_airmass before being used + in calculating Kn. + + Returns + ------- + Kn : numeric + am : numeric + airmass used in the calculation of Kn. am <= max_airmass. """ # short names for equations kt = clearness_index am = airmass - # consider adding - # am = np.maximum(am, 12) # GH 450 + am = np.minimum(am, max_airmass) # GH 450 # powers of kt will be used repeatedly, so compute only once kt2 = kt * kt # about the same as kt ** 2 @@ -1423,7 +1441,7 @@ def _disc_kn(clearness_index, airmass): Knc = 0.866 - 0.122*am + 0.0121*am**2 - 0.000653*am**3 + 1.4e-05*am**4 Kn = Knc - delta_kn - return Kn + return Kn, am def dirint(ghi, solar_zenith, times, pressure=101325., use_delta_kt_prime=True, @@ -1934,7 +1952,7 @@ def _gti_dirint_lt_90(poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth, # calculate kt and DNI from GTI kt = clearness_index(poa_global_i, aoi, I0) # kt from Marion eqn 2 - disc_dni = np.maximum(_disc_kn(kt, airmass) * I0, 0) + disc_dni = np.maximum(_disc_kn(kt, airmass)[0] * I0, 0) kt_prime = clearness_index_zenith_independent(kt, airmass) # dirint DNI in Marion eqn 3 dni = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith, @@ -2017,7 +2035,7 @@ def _gti_dirint_gte_90(poa_global, aoi, solar_zenith, solar_azimuth, airmass = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966') airmass = atmosphere.get_absolute_airmass(airmass, pressure) kt = kt_prime_gte_90 * _kt_kt_prime_factor(airmass) - disc_dni = np.maximum(_disc_kn(kt, airmass) * I0, 0) + disc_dni = np.maximum(_disc_kn(kt, airmass)[0] * I0, 0) dni_gte_90 = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith, False, temp_dew) diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index 0e1f1ece09..cd1c5b168d 100644 --- a/pvlib/test/test_irradiance.py +++ b/pvlib/test/test_irradiance.py @@ -404,28 +404,57 @@ def test_disc_min_cos_zenith_max_zenith(): times = pd.DatetimeIndex(['2016-07-19 06:11:00'], tz='America/Phoenix') out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times) expected = pd.DataFrame(np.array( - [[0.00000000e+00, 1.16046346e-02, 3.63954476e+01]]), + [[0.00000000e+00, 1.16046346e-02, 12.0]]), columns=columns, index=times) assert_frame_equal(out, expected) + # max_zenith and/or max_airmass keep these results reasonable out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times, min_cos_zenith=0) expected = pd.DataFrame(np.array( - [[0.00000000e+00, 1.0, 3.63954476e+01]]), + [[0.00000000e+00, 1.0, 12.0]]), columns=columns, index=times) assert_frame_equal(out, expected) + # still get reasonable values because of max_airmass=12 limit out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times, max_zenith=100) expected = pd.DataFrame(np.array( - [[6.68577449e+03, 1.16046346e-02, 3.63954476e+01]]), + [[0., 1.16046346e-02, 12.0]]), columns=columns, index=times) assert_frame_equal(out, expected) + # still get reasonable values because of max_airmass=12 limit out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times, min_cos_zenith=0, max_zenith=100) expected = pd.DataFrame(np.array( - [[7.21238390e+03, 1.00000000e+00, 3.63954476e+01]]), + [[277.50185968, 1.0, 12.0]]), + columns=columns, index=times) + assert_frame_equal(out, expected) + + # max_zenith keeps this result reasonable + out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times, + min_cos_zenith=0, max_airmass=100) + expected = pd.DataFrame(np.array( + [[0.00000000e+00, 1.0, 36.39544757]]), + columns=columns, index=times) + assert_frame_equal(out, expected) + + # allow zenith to be close to 90 and airmass to be infinite + # and we get crazy values + out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times, + max_zenith=100, max_airmass=100) + expected = pd.DataFrame(np.array( + [[6.68577449e+03, 1.16046346e-02, 3.63954476e+01]]), + columns=columns, index=times) + assert_frame_equal(out, expected) + + # allow min cos zenith to be 0, zenith to be close to 90, + # and airmass to be very big and we get even higher DNI values + out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times, + min_cos_zenith=0, max_zenith=100, max_airmass=100) + expected = pd.DataFrame(np.array( + [[7.21238390e+03, 1., 3.63954476e+01]]), columns=columns, index=times) assert_frame_equal(out, expected) @@ -497,13 +526,18 @@ def test_dirint_min_cos_zenith_max_zenith(): expected = pd.Series([0.0, 0.0], index=times, name='dni') assert_series_equal(out, expected) - out = irradiance.dirint(ghi, solar_zenith, times, max_zenith=100) - expected = pd.Series([862.198, 848.387], index=times, name='dni') + out = irradiance.dirint(ghi, solar_zenith, times, max_zenith=90) + expected = pd.Series([0.0, 0.0], index=times, name='dni') + assert_series_equal(out, expected, check_less_precise=True) + + out = irradiance.dirint(ghi, solar_zenith, times, min_cos_zenith=0, + max_zenith=90) + expected = pd.Series([0.0, 144.264507], index=times, name='dni') assert_series_equal(out, expected, check_less_precise=True) out = irradiance.dirint(ghi, solar_zenith, times, min_cos_zenith=0, max_zenith=100) - expected = pd.Series([147655.5994, 3749.8542], index=times, name='dni') + expected = pd.Series([0.0, 144.264507], index=times, name='dni') assert_series_equal(out, expected, check_less_precise=True) @@ -670,13 +704,13 @@ def test_dirindex_min_cos_zenith_max_zenith(): assert_series_equal(out, expected) out = irradiance.dirindex(ghi, ghi_clearsky, dni_clearsky, solar_zenith, - times, max_zenith=100) - expected = pd.Series([0., 5.], index=times) + times, max_zenith=90) + expected = pd.Series([nan, nan], index=times) assert_series_equal(out, expected) out = irradiance.dirindex(ghi, ghi_clearsky, dni_clearsky, solar_zenith, times, min_cos_zenith=0, max_zenith=100) - expected = pd.Series([0., 5.], index=times) + expected = pd.Series([nan, 5.], index=times) assert_series_equal(out, expected)