Skip to content

"Numpy-isation" of irradiance decomposition functions #1456

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

Open
wants to merge 16 commits into
base: main
Choose a base branch
from
Open

"Numpy-isation" of irradiance decomposition functions #1456

Changes from 9 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
72af610
add a todo in disc function
alexisstgc May 10, 2022
e37c4bb
vectorize disc
alexisstgc May 12, 2022
a3f3f8e
vectorize _delta_kt_prime_dirint
alexisstgc May 12, 2022
273ee88
put shape as function argument in _delte_kt_prime_dirint
alexisstgc May 12, 2022
bb51082
vectorize _temp_dew_dirint
alexisstgc May 12, 2022
6632758
improve docstrings
alexisstgc May 12, 2022
f941da5
vectorize _dirint_coeffs
alexisstgc May 12, 2022
e0c8d9a
vectorize _dirint_bins
alexisstgc May 12, 2022
105c49f
vectorize dirint
alexisstgc May 12, 2022
5ab066a
avoid dividing by zero in dirindex
alexisstgc May 12, 2022
17bee9b
add PVLIB_IRR_DTYPE env variable to choose precision
alexisstgc May 19, 2022
b8631e4
add getenv import
alexisstgc May 19, 2022
09c1b42
fix keyword argument call in _dirint_coeffs
alexisstgc May 19, 2022
201e2b0
stickler code quality
alexisstgc May 20, 2022
a357465
stickler code quality #2
alexisstgc May 20, 2022
5bf9049
stickler
alexisstgc May 20, 2022
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
115 changes: 76 additions & 39 deletions pvlib/irradiance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1394,11 +1394,13 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
--------
dirint
"""
# TODO: this function should be vectorized
# this is the I0 calculation from the reference
# SSC uses solar constant = 1367.0 (checked 2018 08 15)
I0 = get_extra_radiation(datetime_or_doy, 1370., 'spencer')

# Considering the extra radiation is only time dependent, broadcast it to ghi's dimensions
I0 = np.broadcast_to(I0, ghi.shape).astype(np.float32)
Copy link
Member

Choose a reason for hiding this comment

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

I don't believe that all pvlib users would prefer to use float32 - this seems to be a specific requirement of your use case

Copy link
Member

Choose a reason for hiding this comment

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

Is there any downside to using float32 here? Is I0 known with sufficient accuracy to justify float64?


kt = clearness_index(ghi, solar_zenith, I0, min_cos_zenith=min_cos_zenith,
max_clearness_index=1)

Expand Down Expand Up @@ -1543,7 +1545,7 @@ def dirint(ghi, solar_zenith, times, pressure=101325., use_delta_kt_prime=True,
Global Horizontal to Direct Normal Insolation", Technical Report No.
SERI/TR-215-3087, Golden, CO: Solar Energy Research Institute, 1987.
"""

shape = ghi.shape
disc_out = disc(ghi, solar_zenith, times, pressure=pressure,
min_cos_zenith=min_cos_zenith, max_zenith=max_zenith)
airmass = disc_out['airmass']
Expand All @@ -1552,10 +1554,10 @@ def dirint(ghi, solar_zenith, times, pressure=101325., use_delta_kt_prime=True,
kt_prime = clearness_index_zenith_independent(
kt, airmass, max_clearness_index=1)
delta_kt_prime = _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime,
times)
w = _temp_dew_dirint(temp_dew, times)
shape)
w = _temp_dew_dirint(temp_dew, shape)

dirint_coeffs = _dirint_coeffs(times, kt_prime, solar_zenith, w,
dirint_coeffs = _dirint_coeffs(kt_prime, solar_zenith, w,
delta_kt_prime)

# Perez eqn 5
Expand All @@ -1581,51 +1583,75 @@ def _dirint_from_dni_ktprime(dni, kt_prime, solar_zenith, use_delta_kt_prime,
return dni_dirint


def _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime, times):
def _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime, shape):
"""
Calculate delta_kt_prime (Perez eqn 2 and eqn 3), or return a default value
for use with :py:func:`_dirint_bins`.

Parameters
----------
kt_prime : Zenith-independent clearness index
use_delta_kt_prime : Boolean flag whether to use calculate the stability
index or to use the default value
shape : Shape of the input data

Returns
----------
delta_kt_prime : Stability index
"""
if use_delta_kt_prime:
# Perez eqn 2
kt_next = kt_prime.shift(-1)
kt_previous = kt_prime.shift(1)
# replace nan with values that implement Perez Eq 3 for first and last
# positions. Use kt_previous and kt_next to handle series of length 1
kt_next.iloc[-1] = kt_previous.iloc[-1]
kt_previous.iloc[0] = kt_next.iloc[0]
delta_kt_prime = 0.5 * ((kt_prime - kt_next).abs().add(
(kt_prime - kt_previous).abs(),
fill_value=0))
kt_diffp1 = np.abs(np.diff(kt_prime, axis=-1))
kt_diffm1 = np.flip(np.abs(np.diff(np.flip(kt_prime, axis=-1), axis=-1)), axis=-1)

kt_next = np.empty(shape, dtype=np.float32)
# work only on last dimension
kt_next[..., :-1] = kt_diffp1
kt_next[..., -1] = kt_diffm1[..., -1]

kt_previous = np.empty(shape, dtype=np.float32)
# work only on last dimension
kt_previous[..., 1:] = kt_diffm1
kt_previous[..., 0] = kt_diffp1[..., 0]

delta_kt_prime = 0.5 * (np.abs(kt_prime - kt_next) + np.abs(kt_prime - kt_previous))
else:
# do not change unless also modifying _dirint_bins
delta_kt_prime = pd.Series(-1, index=times)
delta_kt_prime = np.full(shape, -1)
return delta_kt_prime


def _temp_dew_dirint(temp_dew, times):
def _temp_dew_dirint(temp_dew, shape):
"""
Calculate precipitable water from surface dew point temp (Perez eqn 4),
or return a default value for use with :py:func:`_dirint_bins`.

Parameters
----------
temp_dew : Surface dew point temperature array
shape : Shape of the input data

Returns
----------
w : Precipitable water estimated from surface dew-point temperature
"""
if temp_dew is not None:
# Perez eqn 4
w = pd.Series(np.exp(0.07 * temp_dew - 0.075), index=times)
w = np.exp(0.07 * temp_dew - 0.075)
else:
# do not change unless also modifying _dirint_bins
w = pd.Series(-1, index=times)
w = np.full(shape, -1)
return w


def _dirint_coeffs(times, kt_prime, solar_zenith, w, delta_kt_prime):
def _dirint_coeffs(kt_prime, solar_zenith, w, delta_kt_prime):
"""
Determine the DISC to DIRINT multiplier `dirint_coeffs`.

dni = disc_out['dni'] * dirint_coeffs

Parameters
----------
times : pd.DatetimeIndex
kt_prime : Zenith-independent clearness index
solar_zenith : Solar zenith angle
w : precipitable water estimated from surface dew-point temperature
Expand All @@ -1635,31 +1661,37 @@ def _dirint_coeffs(times, kt_prime, solar_zenith, w, delta_kt_prime):
-------
dirint_coeffs : array-like
"""
kt_prime_bin, zenith_bin, w_bin, delta_kt_prime_bin = \
_dirint_bins(times, kt_prime, solar_zenith, w, delta_kt_prime)
kt_prime_bin, zenith_bin, w_bin, delta_kt_prime_bin = _dirint_bins(
kt_prime=kt_prime,
solar_zenith=solar_zenith,
w=w,
delta_kt_prime=delta_kt_prime,
)

# get the coefficients
coeffs = _get_dirint_coeffs()
coeffs = _get_dirint_coeffs().astype(np.float32)

# subtract 1 to account for difference between MATLAB-style bin
# assignment and Python-style array lookup.
dirint_coeffs = coeffs[kt_prime_bin-1, zenith_bin-1,
delta_kt_prime_bin-1, w_bin-1]

# convert unassigned bins to nan
dirint_coeffs = np.where((kt_prime_bin == 0) | (zenith_bin == 0) |
(w_bin == 0) | (delta_kt_prime_bin == 0),
np.nan, dirint_coeffs)
dirint_coeffs = coeffs[kt_prime_bin - 1, zenith_bin - 1,
delta_kt_prime_bin - 1, w_bin - 1]

# convert unassigned bins to 1, instead of putting nan originally
# whenever the dirint coeff is not known it should be preferred to fall back
# to disc values and therefore have a coeff of 1.
dirint_coeffs[
(kt_prime_bin == 0) | (zenith_bin == 0) |
(w_bin == 0) | (delta_kt_prime_bin == 0)
] = 1
return dirint_coeffs


def _dirint_bins(times, kt_prime, zenith, w, delta_kt_prime):
def _dirint_bins(kt_prime, zenith, w, delta_kt_prime):
"""
Determine the bins for the DIRINT coefficients.

Parameters
----------
times : pd.DatetimeIndex
kt_prime : Zenith-independent clearness index
zenith : Solar zenith angle
w : precipitable water estimated from surface dew-point temperature
Expand All @@ -1669,11 +1701,12 @@ def _dirint_bins(times, kt_prime, zenith, w, delta_kt_prime):
-------
tuple of kt_prime_bin, zenith_bin, w_bin, delta_kt_prime_bin
"""
shape = kt_prime.shape
# @wholmgren: the following bin assignments use MATLAB's 1-indexing.
# Later, we'll subtract 1 to conform to Python's 0-indexing.

# Create kt_prime bins
kt_prime_bin = pd.Series(0, index=times, dtype=np.int64)
kt_prime_bin = np.zeros(shape, dtype=np.int8)
kt_prime_bin[(kt_prime >= 0) & (kt_prime < 0.24)] = 1
kt_prime_bin[(kt_prime >= 0.24) & (kt_prime < 0.4)] = 2
kt_prime_bin[(kt_prime >= 0.4) & (kt_prime < 0.56)] = 3
Expand All @@ -1682,7 +1715,7 @@ def _dirint_bins(times, kt_prime, zenith, w, delta_kt_prime):
kt_prime_bin[(kt_prime >= 0.8) & (kt_prime <= 1)] = 6

# Create zenith angle bins
zenith_bin = pd.Series(0, index=times, dtype=np.int64)
zenith_bin = np.zeros(shape, dtype=np.int8)
zenith_bin[(zenith >= 0) & (zenith < 25)] = 1
zenith_bin[(zenith >= 25) & (zenith < 40)] = 2
zenith_bin[(zenith >= 40) & (zenith < 55)] = 3
Expand All @@ -1691,15 +1724,15 @@ def _dirint_bins(times, kt_prime, zenith, w, delta_kt_prime):
zenith_bin[(zenith >= 80)] = 6

# Create the bins for w based on dew point temperature
w_bin = pd.Series(0, index=times, dtype=np.int64)
w_bin = np.zeros(shape, dtype=np.int8)
w_bin[(w >= 0) & (w < 1)] = 1
w_bin[(w >= 1) & (w < 2)] = 2
w_bin[(w >= 2) & (w < 3)] = 3
w_bin[(w >= 3)] = 4
w_bin[(w == -1)] = 5

# Create delta_kt_prime binning.
delta_kt_prime_bin = pd.Series(0, index=times, dtype=np.int64)
delta_kt_prime_bin = np.zeros(shape, dtype=np.int8)
delta_kt_prime_bin[(delta_kt_prime >= 0) & (delta_kt_prime < 0.015)] = 1
delta_kt_prime_bin[(delta_kt_prime >= 0.015) &
(delta_kt_prime < 0.035)] = 2
Expand Down Expand Up @@ -1800,9 +1833,13 @@ def dirindex(ghi, ghi_clearsky, dni_clearsky, zenith, times, pressure=101325.,
min_cos_zenith=min_cos_zenith,
max_zenith=max_zenith)

dni_dirindex = dni_clearsky * dni_dirint / dni_dirint_clearsky
# Avoid dividing by zero with clearsky. Whenever the dni_dirint_clearsky is zero, the array defaults
# to dni_dirint (which logically should be also zero)
nonzero = dni_dirint_clearsky != 0
dni_dirindex = dni_dirint
dni_dirindex[nonzero] = dni_clearsky[nonzero] * dni_dirint[nonzero] / dni_dirint_clearsky[nonzero]

dni_dirindex[dni_dirindex < 0] = 0.
dni_dirindex[dni_dirindex < 0] = 0.0

return dni_dirindex

Expand Down