Skip to content

add max_airmass=12 kwarg to disc #626

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Dec 11, 2018
Merged

add max_airmass=12 kwarg to disc #626

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
2514d08
add max_airmass=12 kwarg to disc
wholmgren Nov 26, 2018
1b8d39f
more tests, clean up
wholmgren Nov 26, 2018
1a44e32
Merge branch 'master' into discamlim
wholmgren Dec 11, 2018
74bf083
fix bad whatsnew merge
wholmgren Dec 11, 2018
fe3ff01
more merge fixing
wholmgren Dec 11, 2018
7d72e34
move limit so it helps in gti_dirint. reword whatsnew
wholmgren Dec 11, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions docs/sphinx/source/whatsnew/v0.6.1.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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`
Expand Down
34 changes: 26 additions & 8 deletions pvlib/irradiance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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)
Expand Down
54 changes: 44 additions & 10 deletions pvlib/test/test_irradiance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)


Expand Down Expand Up @@ -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)


Expand Down