Skip to content

Incidence Angle Modifier zero instead of np.nan (#338) #339

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 16 commits into from
Aug 1, 2017
Merged

Incidence Angle Modifier zero instead of np.nan (#338) #339

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
24d8fbc
changed np.nan to 0 in functions ashraeiam, physicaliam, sapm_aoi_loss
nextkaufmann Jun 14, 2017
acacf3a
updated docstrings
nextkaufmann Jun 14, 2017
6167d7f
updated tests to account for changes in functions
nextkaufmann Jun 14, 2017
d92f8fb
replace 0-degree aoi in physicaliam with 1e-06, analogous to the calc…
nextkaufmann Jun 14, 2017
fbbff59
updated tests to account for input angles of np.nan
nextkaufmann Jun 16, 2017
51340e2
updated test for ashraeiam to check for possible negative values for …
nextkaufmann Jun 16, 2017
990acf1
updated docstring with info about nan-inputs resulting in nan.
nextkaufmann Jun 16, 2017
3d3e8b5
updated variables-and-symbols-page with allowable range of aoi
nextkaufmann Jul 3, 2017
31461b1
changed whatsnew for v0.5
nextkaufmann Jul 3, 2017
28c5e17
applied @adriesse suggestion to split calculation of tau
nextkaufmann Jul 4, 2017
23e67e1
Merge branch 'master' into modelchain_zero_instead_of_npnan
jkfm Jul 4, 2017
a607a10
forgot to delete unused line
nextkaufmann Jul 4, 2017
a329a8a
typo :-(
nextkaufmann Jul 4, 2017
fd593a2
typo in variables_styles_rules
nextkaufmann Jul 18, 2017
6ca99af
implemented code suggestion by @adriesse
nextkaufmann Jul 18, 2017
a52dba0
Merge remote-tracking branch 'origin/modelchain_zero_instead_of_npnan…
nextkaufmann Jul 18, 2017
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
29 changes: 29 additions & 0 deletions docs/sphinx/source/whatsnew/v0.5.0.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
.. _whatsnew_0500:

v0.5.0 ()
---------


Bug fixes
~~~~~~~~~


Enhancements
~~~~~~~~~~~~


API Changes
~~~~~~~~~~~

* Changes calculation of the Incidence Angle Modifier to return 0 instead of np.nan for angles >= 90°.
This improves the calculation of effective irradiance close to sunrise and sunset. (:issue:`338`)


Documentation
~~~~~~~~~~~~~


Contributors
~~~~~~~~~~~~

* Johannes Kaufmann
2 changes: 1 addition & 1 deletion pvlib/data/variables_style_rules.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ dni;direct normal irradiance
dni_extra;direct normal irradiance at top of atmosphere (extraterrestrial)
dhi;diffuse horizontal irradiance
ghi;global horizontal irradiance
aoi;angle of incidence
aoi;angle of incidence between :math:`90\deg` and :math:`90\deg`
aoi_projection;cos(aoi)
airmass;airmass
airmass_relative;relative airmass
Expand Down
68 changes: 41 additions & 27 deletions pvlib/pvsystem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -700,7 +700,7 @@ def ashraeiam(aoi, b=0.05):
----------
aoi : numeric
The angle of incidence between the module normal vector and the
sun-beam vector in degrees.
sun-beam vector in degrees. Angles of nan will result in nan.

b : float, default 0.05
A parameter to adjust the modifier as a function of angle of
Expand All @@ -712,7 +712,7 @@ def ashraeiam(aoi, b=0.05):
The incident angle modifier calculated as 1-b*(sec(aoi)-1) as
described in [2,3].

Returns nan for all abs(aoi) >= 90 and for all IAM values that
Returns zeros for all abs(aoi) >= 90 and for all IAM values that
would be less than 0.

References
Expand All @@ -735,7 +735,7 @@ def ashraeiam(aoi, b=0.05):

iam = 1 - b*((1/np.cos(np.radians(aoi)) - 1))

iam = np.where(np.abs(aoi) >= 90, np.nan, iam)
iam = np.where(np.abs(aoi) >= 90, 0, iam)
iam = np.maximum(0, iam)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is the np.maximum line still necessary?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is, because for angles close to +/-90° the function would return a negative value for the iam (the intercept for b=0.05 is at 87.21°: https://www.desmos.com/calculator/y8wvzrx3b9).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, the np.maximum line is needed


if isinstance(iam, pd.Series):
Expand Down Expand Up @@ -764,7 +764,8 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
----------
aoi : numeric
The angle of incidence between the module normal vector and the
sun-beam vector in degrees.
sun-beam vector in degrees. Angles of 0 are replaced with 1e-06
to ensure non-nan results. Angles of nan will result in nan.

n : numeric, default 1.526
The effective index of refraction (unitless). Reference [1]
Expand Down Expand Up @@ -814,27 +815,40 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
spa
ashraeiam
'''
zeroang = 1e-06

aoi = np.where(aoi == 0, zeroang, aoi)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why did you add this line? I understand what it does, but not why we need it now when we didn't previously.

It's not clear to me that we actually need zeroang at all. I think its only purpose is to avoid infs, but maybe those be more cleanly handled with np.where or similar.

Copy link
Contributor Author

@jkfm jkfm Jun 16, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I understood it, if the input angle is 0, tau will be nan, which is why we substitute it with zeroang.
As iam is defined as tau(theta) / tau(0°), we need zeroang anyway to calculate tau(0°). That's why I moved the variable definition to the beginning of the function.

I guess we could define tau0 as a constant in the function and handle aoi == 0 with np.where after the calculations are done, but I feel like the current function is more intuitive to understand.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should keep zeroang and for input theta < zeroang, physicaliam should return tau0. Passing theta=0 through will return nan rather than zero. Line 821 computes tau with the leading factor exp(-KL/cos(theta)) which is why we should return tau0 for small theta.

tau0 is computed awkwardly. I'd replace line 829 - 833 with the limit value which we should have implemented in Matlab:
tau0 = exp(-KL)*(1 - ((1-n)/(1+n))^2)
We used zeroang instead in Matlab which is probably why its in pvlib-python

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok, now I see that you've changed the test so that an input of 0 returns 1 instead of np.nan. I should have caught this when I wrote the test.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am surprised that the model would not work for an angle of exactly zero degrees. Could this have something to do with the equation discrepancy mentioned in the comments further above?

Copy link
Contributor Author

@jkfm jkfm Jul 3, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think it has to do with the function discrepancy. The unmodified function for theta_r is not defined for angles > ~40°.

The tau(theta)-function in turn simply has a discontinuity for theta = 0. It seems this was ignored in the paper...


# angle of reflection
thetar_deg = tools.asind(1.0 / n*(tools.sind(aoi)))

tau = (np.exp(- 1.0 * (K*L / tools.cosd(thetar_deg))) *
((1 - 0.5*((((tools.sind(thetar_deg - aoi)) ** 2) /
((tools.sind(thetar_deg + aoi)) ** 2) +
((tools.tand(thetar_deg - aoi)) ** 2) /
((tools.tand(thetar_deg + aoi)) ** 2))))))
# reflectance and transmittance for normal incidence light
rho_zero = ((1-n) / (1+n)) ** 2
tau_zero = np.exp(-K*L)

zeroang = 1e-06
# reflectance for parallel and perpendicular polarized light
rho_para = (tools.tand(thetar_deg - aoi) /
tools.tand(thetar_deg + aoi)) ** 2
rho_perp = (tools.sind(thetar_deg - aoi) /
tools.sind(thetar_deg + aoi)) ** 2

# transmittance for non-normal light
tau = np.exp(-K*L / tools.cosd(thetar_deg))

thetar_deg0 = tools.asind(1.0 / n*(tools.sind(zeroang)))
# iam is ratio of non-normal to normal incidence transmitted light
# after deducting the reflected portion of each
iam = ((1 - (rho_para + rho_perp) / 2) / (1 - rho_zero) * tau / tau_zero)

tau0 = (np.exp(- 1.0 * (K*L / tools.cosd(thetar_deg0))) *
((1 - 0.5*((((tools.sind(thetar_deg0 - zeroang)) ** 2) /
((tools.sind(thetar_deg0 + zeroang)) ** 2) +
((tools.tand(thetar_deg0 - zeroang)) ** 2) /
((tools.tand(thetar_deg0 + zeroang)) ** 2))))))
# angles near zero produce nan, but iam is defined as one
small_angle = 1e-06
iam = np.where(np.abs(aoi) < small_angle, 1.0, iam)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use
tau0 = (np.exp(- 1.0 * (K*L)) * (1 - (1 - n)^2/(1 + n)^2)
instead of the expression on line 826. I'd switch to tau0 when thetar_deg<zeroing.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with your suggestions, especially after having consulted my D&B book. I would suggest calculating reflectance and absorption separately before combining them in order to clarify the code.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I applied both your suggestions correctly. Could you take a look?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It wasn't quite correct, so I thought it might be easier to communicate it this way:

    # angle of refraction
    thetar_deg = tools.asind(tools.sind(aoi) / n)

    # reflectance and transmittance for normal incidence light
    rho_zero = ((1-n) / (1+n)) ** 2
    tau_zero = np.exp(-K*L)

    # reflectance for parallel and perpendicular polarized light
    rho_para = (tools.tand(thetar_deg - aoi) /
                tools.tand(thetar_deg + aoi)) ** 2
    rho_perp = (tools.sind(thetar_deg - aoi) /
                tools.sind(thetar_deg + aoi)) ** 2

    # transmittance for non-normal light
    tau = np.exp(-K*L / tools.cosd(thetar_deg))

    # iam is ratio of non-normal to normal incidence transmitted light
    # after deducting the reflected portion of each
    iam = ((1 - (rho_para + rho_perp) / 2) / (1 - rho_zero)
          * tau / tau_zero )


    # angles near zero produce nan, but iam is defined as one
    small_angle = 1e-06
    iam = np.where(np.abs(aoi) < small_angle, 1.0, iam)

    # angles at 90 degrees can produce tiny negative values, which should be zero
    # this is a result of calculation precision rather than the physical model
    iam = np.where(iam <  0, 0, iam)    

    # for light coming from behind the plane, none can enter the module
    iam = np.where(aoi > 90, 0, iam)

    # negative angles are not logical
    iam = np.where(aoi <  0, np.nan, iam)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the detailed explanation! So does this refer to reference 2 in the Docstring (I don't have access so cannot check)?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ref [1] simply describes what is in ref [2].

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I implemented your suggestion except for setting negative angles to np.nan, to keep consistency with the other IAM functions for now...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consistency has some appeal...

iam = tau / tau0
# angles at 90 degrees can produce tiny negative values, which should be zero
# this is a result of calculation precision rather than the physical model
iam = np.where(iam < 0, 0, iam)

iam = np.where((np.abs(aoi) >= 90) | (iam < 0), np.nan, iam)
# for light coming from behind the plane, none can enter the module
iam = np.where(aoi > 90, 0, iam)

if isinstance(aoi, pd.Series):
iam = pd.Series(iam, index=aoi.index)
Expand Down Expand Up @@ -1465,7 +1479,7 @@ def sapm_aoi_loss(aoi, module, upper=None):
----------
aoi : numeric
Angle of incidence in degrees. Negative input angles will return
nan values.
zeros.

module : dict-like
A dict, Series, or DataFrame defining the SAPM performance
Expand Down Expand Up @@ -1507,7 +1521,7 @@ def sapm_aoi_loss(aoi, module, upper=None):

aoi_loss = np.polyval(aoi_coeff, aoi)
aoi_loss = np.clip(aoi_loss, 0, upper)
aoi_loss = np.where(aoi < 0, np.nan, aoi_loss)
aoi_loss = np.where(aoi < 0, 0, aoi_loss)

if isinstance(aoi, pd.Series):
aoi_loss = pd.Series(aoi_loss, aoi.index)
Expand Down Expand Up @@ -2084,8 +2098,8 @@ def adrinverter(v_dc, p_dc, inverter, vtol=0.10):
See Notes for required keys.

vtol : numeric, default 0.1
A unit-less fraction that determines how far the efficiency model is
allowed to extrapolate beyond the inverter's normal input voltage
A unit-less fraction that determines how far the efficiency model is
allowed to extrapolate beyond the inverter's normal input voltage
operating range. 0.0 <= vtol <= 1.0

Returns
Expand All @@ -2109,21 +2123,21 @@ def adrinverter(v_dc, p_dc, inverter, vtol=0.10):
Column Description
======= ============================================================
p_nom The nominal power value used to normalize all power values,
typically the DC power needed to produce maximum AC power
typically the DC power needed to produce maximum AC power
output, (W).

v_nom The nominal DC voltage value used to normalize DC voltage
values, typically the level at which the highest efficiency
v_nom The nominal DC voltage value used to normalize DC voltage
values, typically the level at which the highest efficiency
is achieved, (V).

pac_max The maximum AC output power value, used to clip the output
pac_max The maximum AC output power value, used to clip the output
if needed, (W).

ce_list This is a list of 9 coefficients that capture the influence
of input voltage and power on inverter losses, and thereby
efficiency.

p_nt ac-power consumed by inverter at night (night tare) to
p_nt ac-power consumed by inverter at night (night tare) to
maintain circuitry required to sense PV array voltage, (W).
======= ============================================================

Expand Down
26 changes: 13 additions & 13 deletions pvlib/test/test_pvsystem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,41 +94,41 @@ def test_systemdef_dict():

@needs_numpy_1_10
def test_ashraeiam():
thetas = np.linspace(-90, 90, 9)
thetas = np.array([-90. , -67.5, -45. , -22.5, 0. , 22.5, 45. , 67.5, 89., 90. , np.nan])
iam = pvsystem.ashraeiam(thetas, .05)
expected = np.array([ nan, 0.9193437 , 0.97928932, 0.99588039, 1. ,
0.99588039, 0.97928932, 0.9193437 , nan])
expected = np.array([ 0, 0.9193437 , 0.97928932, 0.99588039, 1. ,
0.99588039, 0.97928932, 0.9193437 , 0, 0, np.nan])
assert_allclose(iam, expected, equal_nan=True)


@needs_numpy_1_10
def test_PVSystem_ashraeiam():
module_parameters = pd.Series({'b': 0.05})
system = pvsystem.PVSystem(module_parameters=module_parameters)
thetas = np.linspace(-90, 90, 9)
thetas = np.array([-90. , -67.5, -45. , -22.5, 0. , 22.5, 45. , 67.5, 89., 90. , np.nan])
iam = system.ashraeiam(thetas)
expected = np.array([ nan, 0.9193437 , 0.97928932, 0.99588039, 1. ,
0.99588039, 0.97928932, 0.9193437 , nan])
expected = np.array([ 0, 0.9193437 , 0.97928932, 0.99588039, 1. ,
0.99588039, 0.97928932, 0.9193437 , 0, 0, np.nan])
assert_allclose(iam, expected, equal_nan=True)


@needs_numpy_1_10
def test_physicaliam():
thetas = np.linspace(-90, 90, 9)
thetas = np.array([-90. , -67.5, -45. , -22.5, 0. , 22.5, 45. , 67.5, 90. , np.nan])
iam = pvsystem.physicaliam(thetas, 1.526, 0.002, 4)
expected = np.array([ nan, 0.8893998 , 0.98797788, 0.99926198, nan,
0.99926198, 0.98797788, 0.8893998 , nan])
expected = np.array([ 0, 0.8893998, 0.98797788, 0.99926198, 1,
0.99926198, 0.98797788, 0.8893998, 0, np.nan])
assert_allclose(iam, expected, equal_nan=True)


@needs_numpy_1_10
def test_PVSystem_physicaliam():
module_parameters = pd.Series({'K': 4, 'L': 0.002, 'n': 1.526})
system = pvsystem.PVSystem(module_parameters=module_parameters)
thetas = np.linspace(-90, 90, 9)
thetas = np.array([-90. , -67.5, -45. , -22.5, 0. , 22.5, 45. , 67.5, 90. , np.nan])
iam = system.physicaliam(thetas)
expected = np.array([ nan, 0.8893998 , 0.98797788, 0.99926198, nan,
0.99926198, 0.98797788, 0.8893998 , nan])
expected = np.array([ 0, 0.8893998 , 0.98797788, 0.99926198, 1,
0.99926198, 0.98797788, 0.8893998 , 0, np.nan])
assert_allclose(iam, expected, equal_nan=True)


Expand Down Expand Up @@ -239,7 +239,7 @@ def test_PVSystem_sapm_spectral_loss(sapm_module_params):
@pytest.mark.parametrize('aoi,expected', [
(45, 0.9975036250000002),
(np.array([[-30, 30, 100, np.nan]]),
np.array([[np.nan, 1.007572, 0, np.nan]])),
np.array([[0, 1.007572, 0, np.nan]])),
(pd.Series([80]), pd.Series([0.597472]))
])
def test_sapm_aoi_loss(sapm_module_params, aoi, expected):
Expand Down