From 2514d0881066b07ff4804e5e5abf0f6c54580591 Mon Sep 17 00:00:00 2001 From: Will Holmgren Date: Sun, 25 Nov 2018 17:10:34 -0700 Subject: [PATCH 1/5] add max_airmass=12 kwarg to disc --- docs/sphinx/source/whatsnew/v0.6.1.rst | 5 +++++ pvlib/irradiance.py | 11 ++++++---- pvlib/test/test_irradiance.py | 28 +++++++++++++++----------- 3 files changed, 28 insertions(+), 16 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst index 8bce89ac3f..6e671aeffb 100644 --- a/docs/sphinx/source/whatsnew/v0.6.1.rst +++ b/docs/sphinx/source/whatsnew/v0.6.1.rst @@ -24,6 +24,11 @@ API Changes * Changed function name from :py:func:`pvlib.solarposition.get_rise_set_transit` (deprecated) to :py:func:`pvlib.solarposition.sun_rise_set_transit_spa. `sun_rise_set_transit_spa` requires time input to be localized to the specified latitude/longitude. (:issue:`316') +* 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. Output of + :py:func:`pvlib.irradiance.dirint` and `:py:func:`pvlib.irradiance.dirindex`, + are different when `min_cos_zenith` and `max_zenith` kwargs are used. Enhancements diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index eb63371b81..c0a1c7ad5f 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. @@ -1330,6 +1330,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 @@ -1369,6 +1374,7 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325, am = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966') am = atmosphere.get_absolute_airmass(am, pressure) + am = np.minimum(am, max_airmass) # GH 450 Kn = _disc_kn(kt, am) dni = Kn * I0 @@ -1394,9 +1400,6 @@ def _disc_kn(clearness_index, airmass): kt = clearness_index am = airmass - # consider adding - # am = np.maximum(am, 12) # GH 450 - # powers of kt will be used repeatedly, so compute only once kt2 = kt * kt # about the same as kt ** 2 kt3 = kt2 * kt # 5-10x faster than kt ** 3 diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index 5bde6a262a..0912aedb9a 100644 --- a/pvlib/test/test_irradiance.py +++ b/pvlib/test/test_irradiance.py @@ -372,8 +372,7 @@ def test_disc_value(): [676.09497, 0.63776, 3.02102]]) expected = pd.DataFrame(expected_values, columns=columns, index=times) # check the pandas dataframe. check_less_precise is weird - assert_frame_equal(out, expected, check_less_precise=True) - # use np.assert_allclose to check values more clearly + assert_frame_equal(out, expected, check_less_precise=True) # use np.assert_allclose to check values more clearly assert_allclose(out.values, expected_values, atol=1e-5) @@ -398,19 +397,19 @@ 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) 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) out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times, - max_zenith=100) + max_zenith=100, max_airmass=100) expected = pd.DataFrame(np.array( [[6.68577449e+03, 1.16046346e-02, 3.63954476e+01]]), columns=columns, index=times) @@ -419,7 +418,7 @@ def test_disc_min_cos_zenith_max_zenith(): 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) @@ -491,13 +490,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) @@ -664,13 +668,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) From 1b8d39f00f98cb6ab5fc4240006d51b6770350ee Mon Sep 17 00:00:00 2001 From: Will Holmgren Date: Sun, 25 Nov 2018 17:24:35 -0700 Subject: [PATCH 2/5] more tests, clean up --- pvlib/test/test_irradiance.py | 36 ++++++++++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index 0912aedb9a..c8c054c3e3 100644 --- a/pvlib/test/test_irradiance.py +++ b/pvlib/test/test_irradiance.py @@ -372,7 +372,8 @@ def test_disc_value(): [676.09497, 0.63776, 3.02102]]) expected = pd.DataFrame(expected_values, columns=columns, index=times) # check the pandas dataframe. check_less_precise is weird - assert_frame_equal(out, expected, check_less_precise=True) # use np.assert_allclose to check values more clearly + assert_frame_equal(out, expected, check_less_precise=True) + # use np.assert_allclose to check values more clearly assert_allclose(out.values, expected_values, atol=1e-5) @@ -401,6 +402,7 @@ def test_disc_min_cos_zenith_max_zenith(): 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( @@ -408,13 +410,15 @@ def test_disc_min_cos_zenith_max_zenith(): 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, max_airmass=100) + 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( @@ -422,6 +426,32 @@ def test_disc_min_cos_zenith_max_zenith(): 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) + def test_dirint_value(): times = pd.DatetimeIndex(['2014-06-24T12-0700','2014-06-24T18-0700']) From 74bf0838e1a7449db860f8a9704eb2fd080fed03 Mon Sep 17 00:00:00 2001 From: Will Holmgren Date: Mon, 10 Dec 2018 20:38:55 -0700 Subject: [PATCH 3/5] fix bad whatsnew merge --- docs/sphinx/source/whatsnew/v0.6.1.rst | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst index 8df775cbfd..bd626fe0c3 100644 --- a/docs/sphinx/source/whatsnew/v0.6.1.rst +++ b/docs/sphinx/source/whatsnew/v0.6.1.rst @@ -19,17 +19,12 @@ API Changes they will be removed in v0.7. Use the new :py:func:`pvlib.iotools.read_tmy2` and :py:func:`pvlib.iotools.read_tmy3` instead. (:issue:`261`) * Added keyword argument ``horizon`` to :func:`~pvlib.solarposition.pyephem` - and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'`` -* Changed key names for `components` returned from :py:func:`pvlib.clearsky.detect_clearsky` -* Changed function name from :py:func:`pvlib.solarposition.get_rise_set_transit` (deprecated) - to :py:func:`pvlib.solarposition.sun_rise_set_transit_spa. `sun_rise_set_transit_spa` - requires time input to be localized to the specified latitude/longitude. (:issue:`316') + and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'``. * 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. Output of :py:func:`pvlib.irradiance.dirint` and `:py:func:`pvlib.irradiance.dirindex`, are different when `min_cos_zenith` and `max_zenith` kwargs are used. - and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'``. (:issue:`588`) * Changed key names for `components` returned from :py:func:`pvlib.clearsky.detect_clearsky`. (:issue:`596`) From fe3ff01d79b4b5b5e1e398bc0414aa7f3755ffca Mon Sep 17 00:00:00 2001 From: Will Holmgren Date: Mon, 10 Dec 2018 20:41:05 -0700 Subject: [PATCH 4/5] more merge fixing --- docs/sphinx/source/whatsnew/v0.6.1.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst index bd626fe0c3..2406e4608b 100644 --- a/docs/sphinx/source/whatsnew/v0.6.1.rst +++ b/docs/sphinx/source/whatsnew/v0.6.1.rst @@ -20,12 +20,13 @@ API Changes and :py:func:`pvlib.iotools.read_tmy3` instead. (:issue:`261`) * 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. Output of :py:func:`pvlib.irradiance.dirint` and `:py:func:`pvlib.irradiance.dirindex`, are different when `min_cos_zenith` and `max_zenith` kwargs are used. - (:issue:`588`) + (: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` From 7d72e34d4418b877ba5e7cc8c20fde871317b99d Mon Sep 17 00:00:00 2001 From: Will Holmgren Date: Mon, 10 Dec 2018 21:07:54 -0700 Subject: [PATCH 5/5] move limit so it helps in gti_dirint. reword whatsnew --- docs/sphinx/source/whatsnew/v0.6.1.rst | 10 ++++++---- pvlib/irradiance.py | 27 ++++++++++++++++++++------ 2 files changed, 27 insertions(+), 10 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst index 2406e4608b..92c2d48589 100644 --- a/docs/sphinx/source/whatsnew/v0.6.1.rst +++ b/docs/sphinx/source/whatsnew/v0.6.1.rst @@ -23,10 +23,12 @@ API Changes (: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. Output of - :py:func:`pvlib.irradiance.dirint` and `:py:func:`pvlib.irradiance.dirindex`, - are different when `min_cos_zenith` and `max_zenith` kwargs are used. - (:issue:`450`) + 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 46dc30d947..4855187ea4 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -1381,8 +1381,7 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325, if pressure is not None: am = atmosphere.get_absolute_airmass(am, pressure) - am = np.minimum(am, max_airmass) # GH 450 - 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) @@ -1399,14 +1398,30 @@ 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 + 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 kt3 = kt2 * kt # 5-10x faster than kt ** 3 @@ -1426,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, @@ -1937,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, @@ -2020,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)