"""
The ``irradiance`` module contains functions for modeling global
horizontal irradiance, direct normal irradiance, diffuse horizontal
irradiance, and total irradiance under various conditions.
"""
import datetime
from collections import OrderedDict
from functools import partial
import numpy as np
import pandas as pd
from scipy.interpolate import splev
from scipy.optimize import bisect
from pvlib import atmosphere, solarposition, tools
import pvlib # used to avoid dni name collision in complete_irradiance
# see References section of get_ground_diffuse function
SURFACE_ALBEDOS = {'urban': 0.18,
'grass': 0.20,
'fresh grass': 0.26,
'soil': 0.17,
'sand': 0.40,
'snow': 0.65,
'fresh snow': 0.75,
'asphalt': 0.12,
'concrete': 0.30,
'aluminum': 0.85,
'copper': 0.74,
'fresh steel': 0.35,
'dirty steel': 0.08,
'sea': 0.06}
def _handle_extra_radiation_types(datetime_or_doy, epoch_year):
# This block will set the functions that can be used to convert the
# inputs to either day of year or pandas DatetimeIndex, and the
# functions that will yield the appropriate output type. It's
# complicated because there are many day-of-year-like input types,
# and the different algorithms need different types. Maybe you have
# a better way to do it.
if isinstance(datetime_or_doy, pd.DatetimeIndex):
to_doy = tools._pandas_to_doy # won't be evaluated unless necessary
def to_datetimeindex(x): return x # noqa: E306
to_output = partial(pd.Series, index=datetime_or_doy)
elif isinstance(datetime_or_doy, pd.Timestamp):
to_doy = tools._pandas_to_doy
to_datetimeindex = \
tools._datetimelike_scalar_to_datetimeindex
to_output = tools._scalar_out
elif isinstance(datetime_or_doy,
(datetime.date, datetime.datetime, np.datetime64)):
to_doy = tools._datetimelike_scalar_to_doy
to_datetimeindex = \
tools._datetimelike_scalar_to_datetimeindex
to_output = tools._scalar_out
elif np.isscalar(datetime_or_doy): # ints and floats of various types
def to_doy(x): return x # noqa: E306
to_datetimeindex = partial(tools._doy_to_datetimeindex,
epoch_year=epoch_year)
to_output = tools._scalar_out
else: # assume that we have an array-like object of doy
def to_doy(x): return x # noqa: E306
to_datetimeindex = partial(tools._doy_to_datetimeindex,
epoch_year=epoch_year)
to_output = tools._array_out
return to_doy, to_datetimeindex, to_output
[docs]def aoi_projection(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
"""
Calculates the dot product of the sun position unit vector and the surface
normal unit vector; in other words, the cosine of the angle of incidence.
Usage note: When the sun is behind the surface the value returned is
negative. For many uses negative values must be set to zero.
Input all angles in degrees.
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal.
surface_azimuth : numeric
Panel azimuth from north.
solar_zenith : numeric
Solar zenith angle.
solar_azimuth : numeric
Solar azimuth angle.
Returns
-------
projection : numeric
Dot product of panel normal and solar angle.
"""
projection = (
tools.cosd(surface_tilt) * tools.cosd(solar_zenith) +
tools.sind(surface_tilt) * tools.sind(solar_zenith) *
tools.cosd(solar_azimuth - surface_azimuth))
# GH 1185
projection = np.clip(projection, -1, 1)
try:
projection.name = 'aoi_projection'
except AttributeError:
pass
return projection
[docs]def aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
"""
Calculates the angle of incidence of the solar vector on a surface.
This is the angle between the solar vector and the surface normal.
Input all angles in degrees.
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal.
surface_azimuth : numeric
Panel azimuth from north.
solar_zenith : numeric
Solar zenith angle.
solar_azimuth : numeric
Solar azimuth angle.
Returns
-------
aoi : numeric
Angle of incidence in degrees.
"""
projection = aoi_projection(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth)
aoi_value = np.rad2deg(np.arccos(projection))
try:
aoi_value.name = 'aoi'
except AttributeError:
pass
return aoi_value
[docs]def poa_horizontal_ratio(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth):
"""
Calculates the ratio of the beam components of the plane of array
irradiance and the horizontal irradiance.
Input all angles in degrees.
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal.
surface_azimuth : numeric
Panel azimuth from north.
solar_zenith : numeric
Solar zenith angle.
solar_azimuth : numeric
Solar azimuth angle.
Returns
-------
ratio : numeric
Ratio of the plane of array irradiance to the horizontal plane
irradiance
"""
cos_poa_zen = aoi_projection(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth)
cos_solar_zenith = tools.cosd(solar_zenith)
# ratio of tilted and horizontal beam irradiance
ratio = cos_poa_zen / cos_solar_zenith
try:
ratio.name = 'poa_ratio'
except AttributeError:
pass
return ratio
[docs]def beam_component(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
dni):
"""
Calculates the beam component of the plane of array irradiance.
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal.
surface_azimuth : numeric
Panel azimuth from north.
solar_zenith : numeric
Solar zenith angle.
solar_azimuth : numeric
Solar azimuth angle.
dni : numeric
Direct Normal Irradiance
Returns
-------
beam : numeric
Beam component
"""
beam = dni * aoi_projection(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth)
beam = np.maximum(beam, 0)
return beam
[docs]def get_total_irradiance(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth,
dni, ghi, dhi, dni_extra=None, airmass=None,
albedo=0.25, surface_type=None,
model='isotropic',
model_perez='allsitescomposite1990'):
r"""
Determine total in-plane irradiance and its beam, sky diffuse and ground
reflected components, using the specified sky diffuse irradiance model.
.. math::
I_{tot} = I_{beam} + I_{sky diffuse} + I_{ground}
Sky diffuse models include:
* isotropic (default)
* klucher
* haydavies
* reindl
* king
* perez
* perez-driesse
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal. [degree]
surface_azimuth : numeric
Panel azimuth from north. [degree]
solar_zenith : numeric
Solar zenith angle. [degree]
solar_azimuth : numeric
Solar azimuth angle. [degree]
dni : numeric
Direct Normal Irradiance. [W/m2]
ghi : numeric
Global horizontal irradiance. [W/m2]
dhi : numeric
Diffuse horizontal irradiance. [W/m2]
dni_extra : numeric, optional
Extraterrestrial direct normal irradiance. [W/m2]
airmass : numeric, optional
Relative airmass (not adjusted for pressure). [unitless]
albedo : numeric, default 0.25
Ground surface albedo. [unitless]
surface_type : str, optional
Surface type. See :py:func:`~pvlib.irradiance.get_ground_diffuse` for
the list of accepted values.
model : str, default 'isotropic'
Irradiance model. Can be one of ``'isotropic'``, ``'klucher'``,
``'haydavies'``, ``'reindl'``, ``'king'``, ``'perez'``,
``'perez-driesse'``.
model_perez : str, default 'allsitescomposite1990'
Used only if ``model='perez'``. See :py:func:`~pvlib.irradiance.perez`.
Returns
-------
total_irrad : OrderedDict or DataFrame
Contains keys/columns ``'poa_global', 'poa_direct', 'poa_diffuse',
'poa_sky_diffuse', 'poa_ground_diffuse'``.
Notes
-----
Models ``'haydavies'``, ``'reindl'``, ``'perez'`` and ``'perez-driesse'``
require ``'dni_extra'``. Values can be calculated using
:py:func:`~pvlib.irradiance.get_extra_radiation`.
The ``'perez'`` and ``'perez-driesse'`` models require relative airmass
(``airmass``) as input. If ``airmass`` is not provided, it is calculated
using the defaults in :py:func:`~pvlib.atmosphere.get_relative_airmass`.
"""
poa_sky_diffuse = get_sky_diffuse(
surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
dni, ghi, dhi, dni_extra=dni_extra, airmass=airmass, model=model,
model_perez=model_perez)
poa_ground_diffuse = get_ground_diffuse(surface_tilt, ghi, albedo,
surface_type)
aoi_ = aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth)
irrads = poa_components(aoi_, dni, poa_sky_diffuse, poa_ground_diffuse)
return irrads
[docs]def get_sky_diffuse(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth,
dni, ghi, dhi, dni_extra=None, airmass=None,
model='isotropic',
model_perez='allsitescomposite1990'):
r"""
Determine in-plane sky diffuse irradiance component
using the specified sky diffuse irradiance model.
Sky diffuse models include:
* isotropic (default)
* klucher
* haydavies
* reindl
* king
* perez
* perez-driesse
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal. [degree]
surface_azimuth : numeric
Panel azimuth from north. [degree]
solar_zenith : numeric
Solar zenith angle. [degree]
solar_azimuth : numeric
Solar azimuth angle. [degree]
dni : numeric
Direct Normal Irradiance. [W/m2]
ghi : numeric
Global horizontal irradiance. [W/m2]
dhi : numeric
Diffuse horizontal irradiance. [W/m2]
dni_extra : numeric, optional
Extraterrestrial direct normal irradiance. [W/m2]
airmass : numeric, optional
Relative airmass (not adjusted for pressure). [unitless]
model : str, default 'isotropic'
Irradiance model. Can be one of ``'isotropic'``, ``'klucher'``,
``'haydavies'``, ``'reindl'``, ``'king'``, ``'perez'``,
``'perez-driesse'``.
model_perez : str, default 'allsitescomposite1990'
Used only if ``model='perez'``. See :py:func:`~pvlib.irradiance.perez`.
Returns
-------
poa_sky_diffuse : numeric
Sky diffuse irradiance in the plane of array. [W/m2]
Raises
------
ValueError
If model is one of ``'haydavies'``, ``'reindl'``, or ``'perez'`` and
``dni_extra`` is not specified.
Notes
-----
Models ``'haydavies'``, ``'reindl'``, ``'perez'`` and ``'perez-driesse'``
require ``'dni_extra'``. Values can be calculated using
:py:func:`~pvlib.irradiance.get_extra_radiation`.
The ``'perez'`` and ``'perez-driesse'`` models require relative airmass
(``airmass``) as input. If ``airmass`` is not provided, it is calculated
using the defaults in :py:func:`~pvlib.atmosphere.get_relative_airmass`.
"""
model = model.lower()
if dni_extra is None and model in {'haydavies', 'reindl',
'perez', 'perez-driesse'}:
raise ValueError(f'dni_extra is required for model {model}')
if model == 'isotropic':
sky = isotropic(surface_tilt, dhi)
elif model == 'klucher':
sky = klucher(surface_tilt, surface_azimuth, dhi, ghi,
solar_zenith, solar_azimuth)
elif model == 'haydavies':
sky = haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
solar_zenith, solar_azimuth)
elif model == 'reindl':
sky = reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
solar_zenith, solar_azimuth)
elif model == 'king':
sky = king(surface_tilt, dhi, ghi, solar_zenith)
elif model == 'perez':
if airmass is None:
airmass = atmosphere.get_relative_airmass(solar_zenith)
sky = perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
solar_zenith, solar_azimuth, airmass,
model=model_perez)
elif model == 'perez-driesse':
# perez_driesse will calculate its own airmass if needed
sky = perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
solar_zenith, solar_azimuth, airmass)
else:
raise ValueError(f'invalid model selection {model}')
return sky
[docs]def poa_components(aoi, dni, poa_sky_diffuse, poa_ground_diffuse):
r'''
Determine in-plane irradiance components.
Combines DNI with sky diffuse and ground-reflected irradiance to calculate
total, direct and diffuse irradiance components in the plane of array.
Parameters
----------
aoi : numeric
Angle of incidence of solar rays with respect to the module
surface, from :func:`aoi`.
dni : numeric
Direct normal irradiance (W/m^2), as measured from a TMY file or
calculated with a clearsky model.
poa_sky_diffuse : numeric
Diffuse irradiance (W/m^2) in the plane of the modules, as
calculated by a diffuse irradiance translation function
poa_ground_diffuse : numeric
Ground reflected irradiance (W/m^2) in the plane of the modules,
as calculated by an albedo model (eg. :func:`grounddiffuse`)
Returns
-------
irrads : OrderedDict or DataFrame
Contains the following keys:
* ``poa_global`` : Total in-plane irradiance (W/m^2)
* ``poa_direct`` : Total in-plane beam irradiance (W/m^2)
* ``poa_diffuse`` : Total in-plane diffuse irradiance (W/m^2)
* ``poa_sky_diffuse`` : In-plane diffuse irradiance from sky (W/m^2)
* ``poa_ground_diffuse`` : In-plane diffuse irradiance from ground
(W/m^2)
Notes
------
Negative beam irradiation due to aoi :math:`> 90^{\circ}` or AOI
:math:`< 0^{\circ}` is set to zero.
'''
poa_direct = np.maximum(dni * np.cos(np.radians(aoi)), 0)
poa_diffuse = poa_sky_diffuse + poa_ground_diffuse
poa_global = poa_direct + poa_diffuse
irrads = OrderedDict()
irrads['poa_global'] = poa_global
irrads['poa_direct'] = poa_direct
irrads['poa_diffuse'] = poa_diffuse
irrads['poa_sky_diffuse'] = poa_sky_diffuse
irrads['poa_ground_diffuse'] = poa_ground_diffuse
if isinstance(poa_direct, pd.Series):
irrads = pd.DataFrame(irrads)
return irrads
[docs]def get_ground_diffuse(surface_tilt, ghi, albedo=.25, surface_type=None):
r'''
Estimate diffuse irradiance on a tilted surface from ground reflections.
Ground diffuse irradiance is calculated as
.. math::
G_{ground} = GHI \times \rho \times \frac{1 - \cos\beta}{2}
where :math:`\rho` is ``albedo`` and :math:`\beta` is ``surface_tilt``.
Parameters
----------
surface_tilt : numeric
Surface tilt angles in decimal degrees. Tilt must be >=0 and
<=180. The tilt angle is defined as degrees from horizontal
(e.g. surface facing up = 0, surface facing horizon = 90).
ghi : numeric
Global horizontal irradiance. :math:`W/m^2`
albedo : numeric, default 0.25
Ground reflectance, typically 0.1-0.4 for surfaces on Earth
(land), may increase over snow, ice, etc. May also be known as
the reflection coefficient. Must be >=0 and <=1. Will be
overridden if ``surface_type`` is supplied.
surface_type : string, optional
If supplied, overrides ``albedo``. ``surface_type`` can be one of
'urban', 'grass', 'fresh grass', 'snow', 'fresh snow', 'asphalt',
'concrete', 'aluminum', 'copper', 'fresh steel', 'dirty steel',
'sea'.
Returns
-------
grounddiffuse : numeric
Ground reflected irradiance. :math:`W/m^2`
Notes
-----
Table of albedo values by ``surface_type`` are from [2]_, [3]_, [4]_;
see :py:data:`~pvlib.irradiance.SURFACE_ALBEDOS`.
References
----------
.. [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute
solar irradiance on inclined surfaces for building energy simulation"
2007, Solar Energy vol. 81. pp. 254-267.
.. [2] https://www.pvsyst.com/help/albedo.htm Accessed January, 2024.
.. [3] http://en.wikipedia.org/wiki/Albedo Accessed January, 2024.
.. [4] Payne, R. E. "Albedo of the Sea Surface". J. Atmos. Sci., 29,
pp. 959–970, 1972.
:doi:`10.1175/1520-0469(1972)029<0959:AOTSS>2.0.CO;2`
'''
if surface_type is not None:
albedo = SURFACE_ALBEDOS[surface_type]
diffuse_irrad = ghi * albedo * (1 - np.cos(np.radians(surface_tilt))) * 0.5
try:
diffuse_irrad.name = 'diffuse_ground'
except AttributeError:
pass
return diffuse_irrad
[docs]def isotropic(surface_tilt, dhi):
r'''
Determine diffuse irradiance from the sky on a tilted surface using
the isotropic sky model.
.. math::
I_{d} = DHI \frac{1 + \cos\beta}{2}
Hottel and Woertz's model treats the sky as a uniform source of
diffuse irradiance. Thus, the diffuse irradiance from the sky (ground
reflected irradiance is not included in this algorithm) on a tilted
surface can be found from the diffuse horizontal irradiance and the
tilt angle of the surface. A discussion of the origin of the
isotropic model can be found in [2]_.
Parameters
----------
surface_tilt : numeric
Surface tilt angle in decimal degrees. Tilt must be >=0 and
<=180. The tilt angle is defined as degrees from horizontal
(e.g. surface facing up = 0, surface facing horizon = 90)
dhi : numeric
Diffuse horizontal irradiance in W/m^2. DHI must be >=0.
Returns
-------
diffuse : numeric
The sky diffuse component of the solar radiation.
References
----------
.. [1] Loutzenhiser P.G. et al. "Empirical validation of models to
compute solar irradiance on inclined surfaces for building energy
simulation" 2007, Solar Energy vol. 81. pp. 254-267
:doi:`10.1016/j.solener.2006.03.009`
.. [2] Kamphuis, N.R. et al. "Perspectives on the origin, derivation,
meaning, and significance of the isotropic sky model" 2020, Solar
Energy vol. 201. pp. 8-12
:doi:`10.1016/j.solener.2020.02.067`
'''
sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5
return sky_diffuse
[docs]def klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith,
solar_azimuth):
r'''
Determine diffuse irradiance from the sky on a tilted surface
using Klucher's 1979 model
.. math::
I_{d} = DHI \frac{1 + \cos\beta}{2} (1 + F' \sin^3(\beta/2))
(1 + F' \cos^2\theta\sin^3\theta_z)
where
.. math::
F' = 1 - (I_{d0} / GHI)^2
Klucher's 1979 model determines the diffuse irradiance from the sky
(ground reflected irradiance is not included in this algorithm) on a
tilted surface using the surface tilt angle, surface azimuth angle,
diffuse horizontal irradiance, direct normal irradiance, global
horizontal irradiance, extraterrestrial irradiance, sun zenith
angle, and sun azimuth angle.
Parameters
----------
surface_tilt : numeric
Surface tilt angles in decimal degrees. surface_tilt must be >=0
and <=180. The tilt angle is defined as degrees from horizontal
(e.g. surface facing up = 0, surface facing horizon = 90)
surface_azimuth : numeric
Surface azimuth angles in decimal degrees. surface_azimuth must
be >=0 and <=360. The Azimuth convention is defined as degrees
east of north (e.g. North = 0, South=180 East = 90, West = 270).
dhi : numeric
Diffuse horizontal irradiance in W/m^2. DHI must be >=0.
ghi : numeric
Global irradiance in W/m^2. DNI must be >=0.
solar_zenith : numeric
Apparent (refraction-corrected) zenith angles in decimal
degrees. solar_zenith must be >=0 and <=180.
solar_azimuth : numeric
Sun azimuth angles in decimal degrees. solar_azimuth must be >=0
and <=360. The Azimuth convention is defined as degrees east of
north (e.g. North = 0, East = 90, West = 270).
Returns
-------
diffuse : numeric
The sky diffuse component of the solar radiation.
References
----------
.. [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute
solar irradiance on inclined surfaces for building energy simulation"
2007, Solar Energy vol. 81. pp. 254-267
.. [2] Klucher, T.M., 1979. Evaluation of models to predict insolation on
tilted surfaces. Solar Energy 23 (2), 111-114.
'''
# zenith angle with respect to panel normal.
cos_tt = aoi_projection(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth)
cos_tt = np.maximum(cos_tt, 0) # GH 526
# silence warning from 0 / 0
with np.errstate(invalid='ignore'):
F = 1 - ((dhi / ghi) ** 2)
try:
# fails with single point input
F.fillna(0, inplace=True)
except AttributeError:
F = np.where(np.isnan(F), 0, F)
term1 = 0.5 * (1 + tools.cosd(surface_tilt))
term2 = 1 + F * (tools.sind(0.5 * surface_tilt) ** 3)
term3 = 1 + F * (cos_tt ** 2) * (tools.sind(solar_zenith) ** 3)
sky_diffuse = dhi * term1 * term2 * term3
return sky_diffuse
[docs]def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
solar_zenith=None, solar_azimuth=None, projection_ratio=None,
return_components=False):
r'''
Determine diffuse irradiance from the sky on a tilted surface using
Hay & Davies' 1980 model
.. math::
I_{d} = DHI ( A R_b + (1 - A) (\frac{1 + \cos\beta}{2}) )
Hay and Davies' 1980 model determines the diffuse irradiance from
the sky (ground reflected irradiance is not included in this
algorithm) on a tilted surface using the surface tilt angle, surface
azimuth angle, diffuse horizontal irradiance, direct normal
irradiance, extraterrestrial irradiance, sun zenith angle, and sun
azimuth angle.
Parameters
----------
surface_tilt : numeric
Surface tilt angles in decimal degrees. The tilt angle is
defined as degrees from horizontal (e.g. surface facing up = 0,
surface facing horizon = 90)
surface_azimuth : numeric
Surface azimuth angles in decimal degrees. The azimuth
convention is defined as degrees east of north (e.g. North=0,
South=180, East=90, West=270).
dhi : numeric
Diffuse horizontal irradiance in W/m^2.
dni : numeric
Direct normal irradiance in W/m^2.
dni_extra : numeric
Extraterrestrial normal irradiance in W/m^2.
solar_zenith : numeric, optional
Solar apparent (refraction-corrected) zenith angles in decimal
degrees. Must supply ``solar_zenith`` and ``solar_azimuth`` or
supply ``projection_ratio``.
solar_azimuth : numeric, optional
Solar azimuth angles in decimal degrees. Must supply
``solar_zenith`` and ``solar_azimuth`` or supply
``projection_ratio``.
projection_ratio : numeric, optional
Ratio of angle of incidence projection to solar zenith angle
projection. Must supply ``solar_zenith`` and ``solar_azimuth``
or supply ``projection_ratio``.
return_components : bool, default False
Flag used to decide whether to return the calculated diffuse components
or not.
Returns
--------
numeric, OrderedDict, or DataFrame
Return type controlled by `return_components` argument.
If ``return_components=False``, `sky_diffuse` is returned.
If ``return_components=True``, `diffuse_components` is returned.
sky_diffuse : numeric
The sky diffuse component of the solar radiation on a tilted
surface.
diffuse_components : OrderedDict (array input) or DataFrame (Series input)
Keys/columns are:
* sky_diffuse: Total sky diffuse
* isotropic
* circumsolar
* horizon
Notes
------
When supplying ``projection_ratio``, consider constraining its values
when zenith angle approaches 90 degrees or angle of incidence
projection is negative. See code for details.
References
-----------
.. [1] Loutzenhiser P.G. et. al. "Empirical validation of models to
compute solar irradiance on inclined surfaces for building energy
simulation" 2007, Solar Energy vol. 81. pp. 254-267
.. [2] Hay, J.E., Davies, J.A., 1980. Calculations of the solar
radiation incident on an inclined surface. In: Hay, J.E., Won, T.K.
(Eds.), Proc. of First Canadian Solar Radiation Data Workshop, 59.
Ministry of Supply and Services, Canada.
'''
# if necessary, calculate ratio of titled and horizontal beam irradiance
if projection_ratio is None:
cos_tt = aoi_projection(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth)
cos_tt = np.maximum(cos_tt, 0) # GH 526
cos_solar_zenith = tools.cosd(solar_zenith)
Rb = cos_tt / np.maximum(cos_solar_zenith, 0.01745) # GH 432
else:
Rb = projection_ratio
# Anisotropy Index
AI = dni / dni_extra
# these are the () and [] sub-terms of the second term of eqn 7
term1 = 1 - AI
term2 = 0.5 * (1 + tools.cosd(surface_tilt))
poa_isotropic = np.maximum(dhi * term1 * term2, 0)
poa_circumsolar = np.maximum(dhi * (AI * Rb), 0)
sky_diffuse = poa_isotropic + poa_circumsolar
if return_components:
diffuse_components = OrderedDict()
diffuse_components['sky_diffuse'] = sky_diffuse
# Calculate the individual components
diffuse_components['isotropic'] = poa_isotropic
diffuse_components['circumsolar'] = poa_circumsolar
diffuse_components['horizon'] = np.where(
np.isnan(diffuse_components['isotropic']), np.nan, 0.)
if isinstance(sky_diffuse, pd.Series):
diffuse_components = pd.DataFrame(diffuse_components)
return diffuse_components
else:
return sky_diffuse
[docs]def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
solar_zenith, solar_azimuth):
r'''
Determine diffuse irradiance from the sky on a tilted surface using
Reindl's 1990 model
.. math::
I_{d} = DHI (A R_b + (1 - A) (\frac{1 + \cos\beta}{2})
(1 + \sqrt{\frac{I_{hb}}{I_h}} \sin^3(\beta/2)) )
Reindl's 1990 model determines the diffuse irradiance from the sky
(ground reflected irradiance is not included in this algorithm) on a
tilted surface using the surface tilt angle, surface azimuth angle,
diffuse horizontal irradiance, direct normal irradiance, global
horizontal irradiance, extraterrestrial irradiance, sun zenith
angle, and sun azimuth angle.
Parameters
----------
surface_tilt : numeric
Surface tilt angles in decimal degrees. The tilt angle is
defined as degrees from horizontal (e.g. surface facing up = 0,
surface facing horizon = 90)
surface_azimuth : numeric
Surface azimuth angles in decimal degrees. The azimuth
convention is defined as degrees east of north (e.g. North = 0,
South=180 East = 90, West = 270).
dhi : numeric
diffuse horizontal irradiance in W/m^2.
dni : numeric
direct normal irradiance in W/m^2.
ghi: numeric
Global irradiance in W/m^2.
dni_extra : numeric
Extraterrestrial normal irradiance in W/m^2.
solar_zenith : numeric
Apparent (refraction-corrected) zenith angles in decimal degrees.
solar_azimuth : numeric
Sun azimuth angles in decimal degrees. The azimuth convention is
defined as degrees east of north (e.g. North = 0, East = 90,
West = 270).
Returns
-------
poa_sky_diffuse : numeric
The sky diffuse component of the solar radiation.
Notes
-----
The poa_sky_diffuse calculation is generated from the Loutzenhiser et al.
(2007) paper, equation 8. Note that I have removed the beam and ground
reflectance portion of the equation and this generates ONLY the diffuse
radiation from the sky and circumsolar, so the form of the equation
varies slightly from equation 8.
References
----------
.. [1] Loutzenhiser P.G. et. al. "Empirical validation of models to
compute solar irradiance on inclined surfaces for building energy
simulation" 2007, Solar Energy vol. 81. pp. 254-267
.. [2] Reindl, D.T., Beckmann, W.A., Duffie, J.A., 1990a. Diffuse
fraction correlations. Solar Energy 45(1), 1-7.
.. [3] Reindl, D.T., Beckmann, W.A., Duffie, J.A., 1990b. Evaluation of
hourly tilted surface radiation models. Solar Energy 45(1), 9-17.
'''
cos_tt = aoi_projection(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth)
cos_tt = np.maximum(cos_tt, 0) # GH 526
# do not apply cos(zen) limit here (needed for HB below)
cos_solar_zenith = tools.cosd(solar_zenith)
# ratio of titled and horizontal beam irradiance
Rb = cos_tt / np.maximum(cos_solar_zenith, 0.01745) # GH 432
# Anisotropy Index
AI = dni / dni_extra
# DNI projected onto horizontal
HB = dni * cos_solar_zenith
HB = np.maximum(HB, 0)
# these are the () and [] sub-terms of the second term of eqn 8
term1 = 1 - AI
term2 = 0.5 * (1 + tools.cosd(surface_tilt))
with np.errstate(invalid='ignore', divide='ignore'):
hb_to_ghi = np.where(ghi == 0, 0, np.divide(HB, ghi))
term3 = 1 + np.sqrt(hb_to_ghi) * (tools.sind(0.5 * surface_tilt)**3)
sky_diffuse = dhi * (AI * Rb + term1 * term2 * term3)
sky_diffuse = np.maximum(sky_diffuse, 0)
return sky_diffuse
[docs]def king(surface_tilt, dhi, ghi, solar_zenith):
'''
Determine diffuse irradiance from the sky on a tilted surface using
the King model.
King's model determines the diffuse irradiance from the sky (ground
reflected irradiance is not included in this algorithm) on a tilted
surface using the surface tilt angle, diffuse horizontal irradiance,
global horizontal irradiance, and sun zenith angle. Note that this
model is not well documented and has not been published in any
fashion (as of January 2012).
Parameters
----------
surface_tilt : numeric
Surface tilt angles in decimal degrees. The tilt angle is
defined as degrees from horizontal (e.g. surface facing up = 0,
surface facing horizon = 90)
dhi : numeric
Diffuse horizontal irradiance in W/m^2.
ghi : numeric
Global horizontal irradiance in W/m^2.
solar_zenith : numeric
Apparent (refraction-corrected) zenith angles in decimal degrees.
Returns
--------
poa_sky_diffuse : numeric
The diffuse component of the solar radiation.
'''
sky_diffuse = (dhi * (1 + tools.cosd(surface_tilt)) / 2 + ghi *
(0.012 * solar_zenith - 0.04) *
(1 - tools.cosd(surface_tilt)) / 2)
sky_diffuse = np.maximum(sky_diffuse, 0)
return sky_diffuse
[docs]def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
solar_zenith, solar_azimuth, airmass,
model='allsitescomposite1990', return_components=False):
'''
Determine diffuse irradiance from the sky on a tilted surface using
one of the Perez models.
Perez models determine the diffuse irradiance from the sky (ground
reflected irradiance is not included in this algorithm) on a tilted
surface using the surface tilt angle, surface azimuth angle, diffuse
horizontal irradiance, direct normal irradiance, extraterrestrial
irradiance, sun zenith angle, sun azimuth angle, and relative (not
pressure-corrected) airmass. Optionally a selector may be used to
use any of Perez's model coefficient sets.
Parameters
----------
surface_tilt : numeric
Surface tilt angles in decimal degrees. surface_tilt must be >=0
and <=180. The tilt angle is defined as degrees from horizontal
(e.g. surface facing up = 0, surface facing horizon = 90)
surface_azimuth : numeric
Surface azimuth angles in decimal degrees. surface_azimuth must
be >=0 and <=360. The azimuth convention is defined as degrees
east of north (e.g. North = 0, South=180 East = 90, West = 270).
dhi : numeric
Diffuse horizontal irradiance in W/m^2. DHI must be >=0.
dni : numeric
Direct normal irradiance in W/m^2. DNI must be >=0.
dni_extra : numeric
Extraterrestrial normal irradiance in W/m^2.
solar_zenith : numeric
apparent (refraction-corrected) zenith angles in decimal
degrees. solar_zenith must be >=0 and <=180.
solar_azimuth : numeric
Sun azimuth angles in decimal degrees. solar_azimuth must be >=0
and <=360. The azimuth convention is defined as degrees east of
north (e.g. North = 0, East = 90, West = 270).
airmass : numeric
Relative (not pressure-corrected) airmass values. If AM is a
DataFrame it must be of the same size as all other DataFrame
inputs. AM must be >=0 (careful using the 1/sec(z) model of AM
generation)
model : string, default 'allsitescomposite1990'
A string which selects the desired set of Perez coefficients. If
model is not provided as an input, the default, '1990' will be
used. All possible model selections are:
* '1990'
* 'allsitescomposite1990' (same as '1990')
* 'allsitescomposite1988'
* 'sandiacomposite1988'
* 'usacomposite1988'
* 'france1988'
* 'phoenix1988'
* 'elmonte1988'
* 'osage1988'
* 'albuquerque1988'
* 'capecanaveral1988'
* 'albany1988'
return_components : bool, default False
Flag used to decide whether to return the calculated diffuse components
or not.
Returns
--------
numeric, OrderedDict, or DataFrame
Return type controlled by `return_components` argument.
If ``return_components=False``, `sky_diffuse` is returned.
If ``return_components=True``, `diffuse_components` is returned.
sky_diffuse : numeric
The sky diffuse component of the solar radiation on a tilted
surface.
diffuse_components : OrderedDict (array input) or DataFrame (Series input)
Keys/columns are:
* sky_diffuse: Total sky diffuse
* isotropic
* circumsolar
* horizon
References
----------
.. [1] Loutzenhiser P.G. et. al. "Empirical validation of models to
compute solar irradiance on inclined surfaces for building energy
simulation" 2007, Solar Energy vol. 81. pp. 254-267
.. [2] Perez, R., Seals, R., Ineichen, P., Stewart, R., Menicucci, D.,
1987. A new simplified version of the Perez diffuse irradiance model
for tilted surfaces. Solar Energy 39(3), 221-232.
.. [3] Perez, R., Ineichen, P., Seals, R., Michalsky, J., Stewart, R.,
1990. Modeling daylight availability and irradiance components from
direct and global irradiance. Solar Energy 44 (5), 271-289.
.. [4] Perez, R. et. al 1988. "The Development and Verification of the
Perez Diffuse Radiation Model". SAND88-7030
'''
kappa = 1.041 # for solar_zenith in radians
z = np.radians(solar_zenith) # convert to radians
# delta is the sky's "brightness"
delta = dhi * airmass / dni_extra
# epsilon is the sky's "clearness"
with np.errstate(invalid='ignore'):
eps = ((dhi + dni) / dhi + kappa * (z ** 3)) / (1 + kappa * (z ** 3))
# numpy indexing below will not work with a Series
if isinstance(eps, pd.Series):
eps = eps.values
# Perez et al define clearness bins according to the following
# rules. 1 = overcast ... 8 = clear (these names really only make
# sense for small zenith angles, but...) these values will
# eventually be used as indicies for coeffecient look ups
ebin = np.digitize(eps, (0., 1.065, 1.23, 1.5, 1.95, 2.8, 4.5, 6.2))
ebin = np.array(ebin) # GH 642
ebin[np.isnan(eps)] = 0
# correct for 0 indexing in coeffecient lookup
# later, ebin = -1 will yield nan coefficients
ebin -= 1
# The various possible sets of Perez coefficients are contained
# in a subfunction to clean up the code.
F1c, F2c = _get_perez_coefficients(model)
# results in invalid eps (ebin = -1) being mapped to nans
nans = np.array([np.nan, np.nan, np.nan])
F1c = np.vstack((F1c, nans))
F2c = np.vstack((F2c, nans))
F1 = (F1c[ebin, 0] + F1c[ebin, 1] * delta + F1c[ebin, 2] * z)
F1 = np.maximum(F1, 0)
F2 = (F2c[ebin, 0] + F2c[ebin, 1] * delta + F2c[ebin, 2] * z)
A = aoi_projection(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth)
A = np.maximum(A, 0)
B = tools.cosd(solar_zenith)
B = np.maximum(B, tools.cosd(85))
# Calculate Diffuse POA from sky dome
term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))
term2 = F1 * A / B
term3 = F2 * tools.sind(surface_tilt)
sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0)
# we've preserved the input type until now, so don't ruin it!
if isinstance(sky_diffuse, pd.Series):
sky_diffuse[np.isnan(airmass)] = 0
else:
sky_diffuse = np.where(np.isnan(airmass), 0, sky_diffuse)
if return_components:
diffuse_components = OrderedDict()
diffuse_components['sky_diffuse'] = sky_diffuse
# Calculate the different components
diffuse_components['isotropic'] = dhi * term1
diffuse_components['circumsolar'] = dhi * term2
diffuse_components['horizon'] = dhi * term3
# Set values of components to 0 when sky_diffuse is 0
mask = sky_diffuse == 0
if isinstance(sky_diffuse, pd.Series):
diffuse_components = pd.DataFrame(diffuse_components)
diffuse_components.loc[mask] = 0
else:
diffuse_components = {k: np.where(mask, 0, v) for k, v in
diffuse_components.items()}
return diffuse_components
else:
return sky_diffuse
def _calc_delta(dhi, dni_extra, solar_zenith, airmass=None):
'''
Compute the delta parameter, which represents sky dome "brightness"
in the Perez and Perez-Driesse models.
Helper function for perez_driesse transposition.
'''
if airmass is None:
# use the same airmass model as in the original perez work
airmass = atmosphere.get_relative_airmass(solar_zenith,
'kastenyoung1989')
max_airmass = atmosphere.get_relative_airmass(90, 'kastenyoung1989')
airmass = np.where(solar_zenith >= 90, max_airmass, airmass)
return dhi / (dni_extra / airmass)
def _calc_zeta(dhi, dni, zenith):
'''
Compute the zeta parameter, which represents sky dome "clearness"
in the Perez-Driesse model.
Helper function for perez_driesse transposition.
'''
dhi = np.asarray(dhi)
dni = np.asarray(dni)
# first calculate what zeta would be without the kappa correction
# using eq. 5 and eq. 13
with np.errstate(invalid='ignore'):
zeta = dni / (dhi + dni)
zeta = np.where(dhi == 0, 0.0, zeta)
# then apply the kappa correction in a manner analogous to eq. 7
kappa = 1.041
kterm = kappa * np.radians(zenith) ** 3
zeta = zeta / (1 - kterm * (zeta - 1))
return zeta
def _f(i, j, zeta):
'''
Evaluate the quadratic splines corresponding to the
allsitescomposite1990 Perez model look-up table.
Helper function for perez_driesse transposition.
'''
knots = np.array(
[0.000, 0.000, 0.000,
0.061, 0.187, 0.333, 0.487, 0.643, 0.778, 0.839,
1.000, 1.000, 1.000])
coefs = np.array(
[[-0.053, +0.529, -0.028, -0.071, +0.061, -0.019],
[-0.008, +0.588, -0.062, -0.060, +0.072, -0.022],
[+0.131, +0.770, -0.167, -0.026, +0.106, -0.032],
[+0.328, +0.471, -0.216, +0.069, -0.105, -0.028],
[+0.557, +0.241, -0.300, +0.086, -0.085, -0.012],
[+0.861, -0.323, -0.355, +0.240, -0.467, -0.008],
[ 1.212, -1.239, -0.444, +0.305, -0.797, +0.047],
[ 1.099, -1.847, -0.365, +0.275, -1.132, +0.124],
[+0.544, +0.157, -0.213, +0.118, -1.455, +0.292],
[+0.544, +0.157, -0.213, +0.118, -1.455, +0.292],
[+0.000, +0.000, +0.000, +0.000, +0.000, +0.000],
[+0.000, +0.000, +0.000, +0.000, +0.000, +0.000],
[+0.000, +0.000, +0.000, +0.000, +0.000, +0.000]])
coefs = coefs.T.reshape((2, 3, 13))
tck = (knots, coefs[i-1, j-1], 2)
return splev(zeta, tck)
[docs]def perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
solar_zenith, solar_azimuth, airmass=None,
return_components=False):
'''
Determine diffuse irradiance from the sky on a tilted surface using
the continuous Perez-Driesse model.
The Perez-Driesse model [1]_ is a reformulation of the 1990 Perez
model [2]_ that provides continuity of the function and of its first
derivatives. This is achieved by replacing the look-up table of
coefficients with quadratic splines.
Parameters
----------
surface_tilt : numeric
Surface tilt angles in decimal degrees. surface_tilt must be >=0
and <=180. The tilt angle is defined as degrees from horizontal
(e.g. surface facing up = 0, surface facing horizon = 90)
surface_azimuth : numeric
Surface azimuth angles in decimal degrees. surface_azimuth must
be >=0 and <=360. The azimuth convention is defined as degrees
east of north (e.g. North = 0, South=180 East = 90, West = 270).
dhi : numeric
Diffuse horizontal irradiance in W/m^2. dhi must be >=0.
dni : numeric
Direct normal irradiance in W/m^2. dni must be >=0.
dni_extra : numeric
Extraterrestrial normal irradiance in W/m^2.
solar_zenith : numeric
apparent (refraction-corrected) zenith angles in decimal
degrees. solar_zenith must be >=0 and <=180.
solar_azimuth : numeric
Sun azimuth angles in decimal degrees. solar_azimuth must be >=0
and <=360. The azimuth convention is defined as degrees east of
north (e.g. North = 0, East = 90, West = 270).
airmass : numeric, optional
Relative (not pressure-corrected) airmass values. If ``airmass`` is a
DataFrame it must be of the same size as all other DataFrame
inputs. The kastenyoung1989 airmass calculation is used internally
and is also recommended when pre-calculating airmass because
it was used in the original model development.
return_components: bool (optional, default=False)
Flag used to decide whether to return the calculated diffuse components
or not.
Returns
--------
numeric, OrderedDict, or DataFrame
Return type controlled by `return_components` argument.
If ``return_components=False``, `sky_diffuse` is returned.
If ``return_components=True``, `diffuse_components` is returned.
sky_diffuse : numeric
The sky diffuse component of the solar radiation on a tilted
surface.
diffuse_components : OrderedDict (array input) or DataFrame (Series input)
Keys/columns are:
* sky_diffuse: Total sky diffuse
* isotropic
* circumsolar
* horizon
Notes
-----
The Perez-Driesse model can be considered a plug-in replacement for the
1990 Perez model using the ``'allsitescomposite1990'`` coefficient set.
Deviations between the two are very small, as demonstrated in [1]_.
Other coefficient sets are not supported because the 1990 set is
based on the largest and most diverse set of empirical data.
References
----------
.. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the
Perez diffuse sky model for forward and reverse transposition.
Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093`
.. [2] Perez, R., Ineichen, P., Seals, R., Michalsky, J., Stewart, R.,
1990. Modeling daylight availability and irradiance components from
direct and global irradiance. Solar Energy 44 (5), 271-289.
See also
--------
perez
isotropic
haydavies
klucher
reindl
king
'''
# Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct., 2023
delta = _calc_delta(dhi, dni_extra, solar_zenith, airmass)
zeta = _calc_zeta(dhi, dni, solar_zenith)
z = np.radians(solar_zenith)
F1 = _f(1, 1, zeta) + _f(1, 2, zeta) * delta + _f(1, 3, zeta) * z
F2 = _f(2, 1, zeta) + _f(2, 2, zeta) * delta + _f(2, 3, zeta) * z
# note the newly recommended upper limit on F1
F1 = np.clip(F1, 0, 0.9)
# lines after this point are identical to the original perez function
# with some checks removed
A = aoi_projection(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth)
A = np.maximum(A, 0)
B = tools.cosd(solar_zenith)
B = np.maximum(B, tools.cosd(85))
# Calculate Diffuse POA from sky dome
term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))
term2 = F1 * A / B
term3 = F2 * tools.sind(surface_tilt)
sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0)
if return_components:
diffuse_components = OrderedDict()
diffuse_components['sky_diffuse'] = sky_diffuse
# Calculate the different components
diffuse_components['isotropic'] = dhi * term1
diffuse_components['circumsolar'] = dhi * term2
diffuse_components['horizon'] = dhi * term3
if isinstance(sky_diffuse, pd.Series):
diffuse_components = pd.DataFrame(diffuse_components)
return diffuse_components
else:
return sky_diffuse
def _poa_from_ghi(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth,
ghi,
dni_extra, airmass, albedo):
'''
Transposition function that includes decomposition of GHI using the
continuous Erbs-Driesse model.
Helper function for ghi_from_poa_driesse_2023.
'''
# Contributed by Anton Driesse (@adriesse), PV Performance Labs. Nov., 2023
erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra)
dni = erbsout['dni']
dhi = erbsout['dhi']
irrads = get_total_irradiance(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth,
dni, ghi, dhi,
dni_extra, airmass, albedo,
model='perez-driesse')
return irrads['poa_global']
def _ghi_from_poa(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth,
poa_global,
dni_extra, airmass, albedo,
xtol=0.01):
'''
Reverse transposition function that uses the scalar bisection from scipy.
Helper function for ghi_from_poa_driesse_2023.
'''
# Contributed by Anton Driesse (@adriesse), PV Performance Labs. Nov., 2023
# propagate nans and zeros quickly
if np.isnan(poa_global):
return np.nan, False, 0
if poa_global <= 0:
return 0.0, True, 0
# function whose root needs to be found
def poa_error(ghi):
poa_hat = _poa_from_ghi(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth,
ghi,
dni_extra, airmass, albedo)
return poa_hat - poa_global
# calculate an upper bound for ghi using clearness index 1.25
ghi_clear = dni_extra * tools.cosd(solar_zenith)
ghi_high = np.maximum(10, 1.25 * ghi_clear)
try:
result = bisect(poa_error,
a=0,
b=ghi_high,
xtol=xtol,
maxiter=25,
full_output=True,
disp=False,
)
except ValueError:
# this occurs when poa_error has the same sign at both end points
ghi = np.nan
conv = False
niter = -1
else:
ghi = result[0]
conv = result[1].converged
niter = result[1].iterations
return ghi, conv, niter
[docs]def ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth,
poa_global,
dni_extra=None, airmass=None, albedo=0.25,
xtol=0.01,
full_output=False):
'''
Estimate global horizontal irradiance (GHI) from global plane-of-array
(POA) irradiance. This reverse transposition algorithm uses a bisection
search together with the continuous Perez-Driesse transposition and
continuous Erbs-Driesse decomposition models, as described in [1]_.
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal. [degree]
surface_azimuth : numeric
Panel azimuth from north. [degree]
solar_zenith : numeric
Solar zenith angle. [degree]
solar_azimuth : numeric
Solar azimuth angle. [degree]
poa_global : numeric
Plane-of-array global irradiance, aka global tilted irradiance. [W/m^2]
dni_extra : None or numeric, default None
Extraterrestrial direct normal irradiance. [W/m^2]
airmass : None or numeric, default None
Relative airmass (not adjusted for pressure). [unitless]
albedo : numeric, default 0.25
Ground surface albedo. [unitless]
xtol : numeric, default 0.01
Convergence criterion. The estimated GHI will be within xtol of the
true value. Must be positive. [W/m^2]
full_output : boolean, default False
If full_output is False, only ghi is returned, otherwise the return
value is (ghi, converged, niter). (see Returns section for details).
Returns
-------
ghi : numeric
Estimated GHI. [W/m^2]
converged : boolean, optional
Present if full_output=True. Indicates which elements converged
successfully.
niter : integer, optional
Present if full_output=True. Indicates how many bisection iterations
were done.
Notes
-----
Since :py:func:`scipy.optimize.bisect` is not vectorized, high-resolution
time series can be quite slow to process.
References
----------
.. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the
Perez diffuse sky model for forward and reverse transposition.
Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093`
See also
--------
perez_driesse
erbs_driesse
gti_dirint
'''
# Contributed by Anton Driesse (@adriesse), PV Performance Labs. Nov., 2023
if xtol <= 0:
raise ValueError(f"xtol too small ({xtol:g} <= 0)")
ghi_from_poa_array = np.vectorize(_ghi_from_poa)
ghi, conv, niter = ghi_from_poa_array(surface_tilt, surface_azimuth,
solar_zenith, solar_azimuth,
poa_global,
dni_extra, airmass, albedo,
xtol=xtol)
if isinstance(poa_global, pd.Series):
ghi = pd.Series(ghi, poa_global.index)
conv = pd.Series(conv, poa_global.index)
niter = pd.Series(niter, poa_global.index)
if full_output:
return ghi, conv, niter
else:
return ghi
[docs]def clearsky_index(ghi, clearsky_ghi, max_clearsky_index=2.0):
"""
Calculate the clearsky index.
The clearsky index is the ratio of global to clearsky global irradiance.
Negative and non-finite clearsky index values will be truncated to zero.
Parameters
----------
ghi : numeric
Global horizontal irradiance in W/m^2.
clearsky_ghi : numeric
Modeled clearsky GHI
max_clearsky_index : numeric, default 2.0
Maximum value of the clearsky index. The default, 2.0, allows
for over-irradiance events typically seen in sub-hourly data.
Returns
-------
clearsky_index : numeric
Clearsky index
"""
clearsky_index = ghi / clearsky_ghi
# set +inf, -inf, and nans to zero
clearsky_index = np.where(~np.isfinite(clearsky_index), 0,
clearsky_index)
# but preserve nans in the input arrays
input_is_nan = ~np.isfinite(ghi) | ~np.isfinite(clearsky_ghi)
clearsky_index = np.where(input_is_nan, np.nan, clearsky_index)
clearsky_index = np.maximum(clearsky_index, 0)
clearsky_index = np.minimum(clearsky_index, max_clearsky_index)
# preserve input type
if isinstance(ghi, pd.Series):
clearsky_index = pd.Series(clearsky_index, index=ghi.index)
return clearsky_index
[docs]def clearness_index(ghi, solar_zenith, extra_radiation, min_cos_zenith=0.065,
max_clearness_index=2.0):
"""
Calculate the clearness index.
The clearness index is the ratio of global to extraterrestrial
irradiance on a horizontal plane [1]_.
Parameters
----------
ghi : numeric
Global horizontal irradiance in W/m^2.
solar_zenith : numeric
True (not refraction-corrected) solar zenith angle in decimal
degrees.
extra_radiation : numeric
Irradiance incident at the top of the atmosphere
min_cos_zenith : numeric, default 0.065
Minimum value of cos(zenith) to allow when calculating global
clearness index `kt`. Equivalent to zenith = 86.273 degrees.
max_clearness_index : numeric, default 2.0
Maximum value of the clearness index. The default, 2.0, allows
for over-irradiance events typically seen in sub-hourly data.
NREL's SRRL Fortran code used 0.82 for hourly data.
Returns
-------
kt : numeric
Clearness index
References
----------
.. [1] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
Global Horizontal to Direct Normal Insolation", Technical
Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research
Institute, 1987.
"""
cos_zenith = tools.cosd(solar_zenith)
I0h = extra_radiation * np.maximum(cos_zenith, min_cos_zenith)
# consider adding
# with np.errstate(invalid='ignore', divide='ignore'):
# to kt calculation, but perhaps it's good to allow these
# warnings to the users that override min_cos_zenith
kt = ghi / I0h
kt = np.maximum(kt, 0)
kt = np.minimum(kt, max_clearness_index)
return kt
[docs]def clearness_index_zenith_independent(clearness_index, airmass,
max_clearness_index=2.0):
"""
Calculate the zenith angle independent clearness index.
See [1]_ for details.
Parameters
----------
clearness_index : numeric
Ratio of global to extraterrestrial irradiance on a horizontal
plane
airmass : numeric
Airmass
max_clearness_index : numeric, default 2.0
Maximum value of the clearness index. The default, 2.0, allows
for over-irradiance events typically seen in sub-hourly data.
NREL's SRRL Fortran code used 0.82 for hourly data.
Returns
-------
kt_prime : numeric
Zenith independent clearness index
References
----------
.. [1] Perez, R., P. Ineichen, E. Maxwell, R. Seals and A. Zelenka,
(1992). "Dynamic Global-to-Direct Irradiance Conversion Models".
ASHRAE Transactions-Research Series, pp. 354-369
"""
# Perez eqn 1
kt_prime = clearness_index / _kt_kt_prime_factor(airmass)
kt_prime = np.maximum(kt_prime, 0)
kt_prime = np.minimum(kt_prime, max_clearness_index)
return kt_prime
def _kt_kt_prime_factor(airmass):
"""
Calculate the conversion factor between kt and kt prime.
Function is useful because DIRINT and GTI-DIRINT both use this.
"""
# consider adding
# airmass = np.maximum(airmass, 12) # GH 450
return 1.031 * np.exp(-1.4 / (0.9 + 9.4 / airmass)) + 0.1
[docs]def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
min_cos_zenith=0.065, max_zenith=87, max_airmass=12):
"""
Estimate Direct Normal Irradiance from Global Horizontal Irradiance
using the DISC model.
The DISC algorithm converts global horizontal irradiance to direct
normal irradiance through empirical relationships between the global
and direct clearness indices.
The pvlib implementation limits the clearness index to 1.
The original report describing the DISC model [1]_ uses the
relative airmass rather than the absolute (pressure-corrected)
airmass. However, the NREL implementation of the DISC model [2]_
uses absolute airmass. PVLib Matlab also uses the absolute airmass.
pvlib python defaults to absolute airmass, but the relative airmass
can be used by supplying `pressure=None`.
Parameters
----------
ghi : numeric
Global horizontal irradiance in W/m^2.
solar_zenith : numeric
True (not refraction-corrected) solar zenith angles in decimal
degrees.
datetime_or_doy : int, float, array, pd.DatetimeIndex
Day of year or array of days of year e.g.
pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
pressure : numeric or None, default 101325
Site pressure in Pascal. Uses absolute (pressure-corrected) airmass
by default. Set to ``None`` to use relative airmass.
min_cos_zenith : numeric, default 0.065
Minimum value of cos(zenith) to allow when calculating global
clearness index `kt`. Equivalent to zenith = 86.273 degrees.
max_zenith : numeric, default 87
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
Contains the following keys:
* ``dni``: The modeled direct normal irradiance
in W/m^2 provided by the
Direct Insolation Simulation Code (DISC) model.
* ``kt``: Ratio of global to extraterrestrial
irradiance on a horizontal plane.
* ``airmass``: Airmass
References
----------
.. [1] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
Global Horizontal to Direct Normal Insolation", Technical
Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research
Institute, 1987.
.. [2] Maxwell, E. "DISC Model", Excel Worksheet.
https://www.nrel.gov/grid/solar-resource/disc.html
See Also
--------
dirint
"""
# 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')
kt = clearness_index(ghi, solar_zenith, I0, min_cos_zenith=min_cos_zenith,
max_clearness_index=1)
am = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
if pressure is not None:
am = atmosphere.get_absolute_airmass(am, pressure)
Kn, am = _disc_kn(kt, am, max_airmass=max_airmass)
dni = Kn * I0
bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0)
dni = np.where(bad_values, 0, dni)
output = OrderedDict()
output['dni'] = dni
output['kt'] = kt
output['airmass'] = am
if isinstance(datetime_or_doy, pd.DatetimeIndex):
output = pd.DataFrame(output, index=datetime_or_doy)
return output
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
is_cloudy = (kt <= 0.6)
# Use Horner's method to compute polynomials efficiently
a = np.where(
is_cloudy,
0.512 + kt*(-1.56 + kt*(2.286 - 2.222*kt)),
-5.743 + kt*(21.77 + kt*(-27.49 + 11.56*kt)))
b = np.where(
is_cloudy,
0.37 + 0.962*kt,
41.4 + kt*(-118.5 + kt*(66.05 + 31.9*kt)))
c = np.where(
is_cloudy,
-0.28 + kt*(0.932 - 2.048*kt),
-47.01 + kt*(184.2 + kt*(-222.0 + 73.81*kt)))
delta_kn = a + b * np.exp(c*am)
Knc = 0.866 + am*(-0.122 + am*(0.0121 + am*(-0.000653 + 1.4e-05*am)))
Kn = Knc - delta_kn
return Kn, am
[docs]def dirint(ghi, solar_zenith, times, pressure=101325., use_delta_kt_prime=True,
temp_dew=None, min_cos_zenith=0.065, max_zenith=87):
"""
Determine DNI from GHI using the DIRINT modification of the DISC
model.
Implements the modified DISC model known as "DIRINT" introduced in
[1]_. DIRINT predicts direct normal irradiance (DNI) from measured
global horizontal irradiance (GHI). DIRINT improves upon the DISC
model by using time-series GHI data and dew point temperature
information. The effectiveness of the DIRINT model improves with
each piece of information provided.
The pvlib implementation limits the clearness index to 1.
Parameters
----------
ghi : array-like
Global horizontal irradiance in W/m^2.
solar_zenith : array-like
True (not refraction-corrected) solar_zenith angles in decimal
degrees.
times : DatetimeIndex
pressure : float or array-like, default 101325.0
The site pressure in Pascal. Pressure may be measured or an
average pressure may be calculated from site altitude.
use_delta_kt_prime : bool, default True
If True, indicates that the stability index delta_kt_prime is
included in the model. The stability index adjusts the estimated
DNI in response to dynamics in the time series of GHI. It is
recommended that delta_kt_prime is not used if the time between
GHI points is 1.5 hours or greater. If use_delta_kt_prime=True,
input data must be Series.
temp_dew : float, or array-like, optional
Surface dew point temperatures, in degrees C. Values of temp_dew
may be numeric or NaN. Any single time period point with a
temp_dew=NaN does not have dew point improvements applied. If
temp_dew is not provided, then dew point improvements are not
applied.
min_cos_zenith : numeric, default 0.065
Minimum value of cos(zenith) to allow when calculating global
clearness index `kt`. Equivalent to zenith = 86.273 degrees.
max_zenith : numeric, default 87
Maximum value of zenith to allow in DNI calculation. DNI will be
set to 0 for times with zenith values greater than `max_zenith`.
Returns
-------
dni : array-like
The modeled direct normal irradiance in W/m^2 provided by the
DIRINT model.
Notes
-----
DIRINT model requires time series data (ie. one of the inputs must
be a vector of length > 2).
References
----------
.. [1] Perez, R., P. Ineichen, E. Maxwell, R. Seals and A. Zelenka,
(1992). "Dynamic Global-to-Direct Irradiance Conversion Models".
ASHRAE Transactions-Research Series, pp. 354-369
.. [2] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
Global Horizontal to Direct Normal Insolation", Technical Report No.
SERI/TR-215-3087, Golden, CO: Solar Energy Research Institute, 1987.
"""
disc_out = disc(ghi, solar_zenith, times, pressure=pressure,
min_cos_zenith=min_cos_zenith, max_zenith=max_zenith)
airmass = disc_out['airmass']
kt = disc_out['kt']
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)
dirint_coeffs = _dirint_coeffs(times, kt_prime, solar_zenith, w,
delta_kt_prime)
# Perez eqn 5
dni = disc_out['dni'] * dirint_coeffs
return dni
def _dirint_from_dni_ktprime(dni, kt_prime, solar_zenith, use_delta_kt_prime,
temp_dew):
"""
Calculate DIRINT DNI from supplied DISC DNI and Kt'.
Supports :py:func:`gti_dirint`
"""
times = dni.index
delta_kt_prime = _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime,
times)
w = _temp_dew_dirint(temp_dew, times)
dirint_coeffs = _dirint_coeffs(times, kt_prime, solar_zenith, w,
delta_kt_prime)
dni_dirint = dni * dirint_coeffs
return dni_dirint
def _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime, times):
"""
Calculate delta_kt_prime (Perez eqn 2 and eqn 3), or return a default value
for use with :py:func:`_dirint_bins`.
"""
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))
else:
# do not change unless also modifying _dirint_bins
delta_kt_prime = pd.Series(-1, index=times)
return delta_kt_prime
def _temp_dew_dirint(temp_dew, times):
"""
Calculate precipitable water from surface dew point temp (Perez eqn 4),
or return a default value for use with :py:func:`_dirint_bins`.
"""
if temp_dew is not None:
# Perez eqn 4
w = pd.Series(np.exp(0.07 * temp_dew - 0.075), index=times)
else:
# do not change unless also modifying _dirint_bins
w = pd.Series(-1, index=times)
return w
def _dirint_coeffs(times, 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
delta_kt_prime : stability index
Returns
-------
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)
# get the coefficients
coeffs = _get_dirint_coeffs()
# 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)
return dirint_coeffs
def _dirint_bins(times, 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
delta_kt_prime : stability index
Returns
-------
tuple of kt_prime_bin, zenith_bin, w_bin, delta_kt_prime_bin
"""
# @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[(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
kt_prime_bin[(kt_prime >= 0.56) & (kt_prime < 0.7)] = 4
kt_prime_bin[(kt_prime >= 0.7) & (kt_prime < 0.8)] = 5
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[(zenith >= 0) & (zenith < 25)] = 1
zenith_bin[(zenith >= 25) & (zenith < 40)] = 2
zenith_bin[(zenith >= 40) & (zenith < 55)] = 3
zenith_bin[(zenith >= 55) & (zenith < 70)] = 4
zenith_bin[(zenith >= 70) & (zenith < 80)] = 5
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[(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[(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
delta_kt_prime_bin[(delta_kt_prime >= 0.035) & (delta_kt_prime < 0.07)] = 3
delta_kt_prime_bin[(delta_kt_prime >= 0.07) & (delta_kt_prime < 0.15)] = 4
delta_kt_prime_bin[(delta_kt_prime >= 0.15) & (delta_kt_prime < 0.3)] = 5
delta_kt_prime_bin[(delta_kt_prime >= 0.3) & (delta_kt_prime <= 1)] = 6
delta_kt_prime_bin[delta_kt_prime == -1] = 7
return kt_prime_bin, zenith_bin, w_bin, delta_kt_prime_bin
[docs]def dirindex(ghi, ghi_clearsky, dni_clearsky, zenith, times, pressure=101325.,
use_delta_kt_prime=True, temp_dew=None, min_cos_zenith=0.065,
max_zenith=87):
"""
Determine DNI from GHI using the DIRINDEX model.
The DIRINDEX model [1]_ modifies the DIRINT model implemented in
:py:func:`pvlib.irradiance.dirint` by taking into account information
from a clear sky model. It is recommended that ``ghi_clearsky`` be
calculated using the Ineichen clear sky model
:py:func:`pvlib.clearsky.ineichen` with ``perez_enhancement=True``.
The pvlib implementation limits the clearness index to 1.
Parameters
----------
ghi : array-like
Global horizontal irradiance in W/m^2.
ghi_clearsky : array-like
Global horizontal irradiance from clear sky model, in W/m^2.
dni_clearsky : array-like
Direct normal irradiance from clear sky model, in W/m^2.
zenith : array-like
True (not refraction-corrected) zenith angles in decimal
degrees. If Z is a vector it must be of the same size as all
other vector inputs. Z must be >=0 and <=180.
times : DatetimeIndex
pressure : float or array-like, default 101325.0
The site pressure in Pascal. Pressure may be measured or an
average pressure may be calculated from site altitude.
use_delta_kt_prime : bool, default True
If True, indicates that the stability index delta_kt_prime is
included in the model. The stability index adjusts the estimated
DNI in response to dynamics in the time series of GHI. It is
recommended that delta_kt_prime is not used if the time between
GHI points is 1.5 hours or greater. If use_delta_kt_prime=True,
input data must be Series.
temp_dew : float, or array-like, optional
Surface dew point temperatures, in degrees C. Values of temp_dew
may be numeric or NaN. Any single time period point with a
temp_dew=NaN does not have dew point improvements applied. If
temp_dew is not provided, then dew point improvements are not
applied.
min_cos_zenith : numeric, default 0.065
Minimum value of cos(zenith) to allow when calculating global
clearness index `kt`. Equivalent to zenith = 86.273 degrees.
max_zenith : numeric, default 87
Maximum value of zenith to allow in DNI calculation. DNI will be
set to 0 for times with zenith values greater than `max_zenith`.
Returns
-------
dni : array-like
The modeled direct normal irradiance in W/m^2.
Notes
-----
DIRINDEX model requires time series data (ie. one of the inputs must
be a vector of length > 2).
References
----------
.. [1] Perez, R., Ineichen, P., Moore, K., Kmiecik, M., Chain, C., George,
R., & Vignola, F. (2002). A new operational model for satellite-derived
irradiances: description and validation. Solar Energy, 73(5), 307-317.
"""
dni_dirint = dirint(ghi, zenith, times, pressure=pressure,
use_delta_kt_prime=use_delta_kt_prime,
temp_dew=temp_dew, min_cos_zenith=min_cos_zenith,
max_zenith=max_zenith)
dni_dirint_clearsky = dirint(ghi_clearsky, zenith, times,
pressure=pressure,
use_delta_kt_prime=use_delta_kt_prime,
temp_dew=temp_dew,
min_cos_zenith=min_cos_zenith,
max_zenith=max_zenith)
dni_dirindex = dni_clearsky * dni_dirint / dni_dirint_clearsky
dni_dirindex[dni_dirindex < 0] = 0.
return dni_dirindex
[docs]def gti_dirint(poa_global, aoi, solar_zenith, solar_azimuth, times,
surface_tilt, surface_azimuth, pressure=101325.,
use_delta_kt_prime=True, temp_dew=None, albedo=.25,
model='perez', model_perez='allsitescomposite1990',
calculate_gt_90=True, max_iterations=30):
"""
Determine GHI, DNI, DHI from POA global using the GTI DIRINT model.
The GTI DIRINT model is described in [1]_.
.. warning::
Model performance is poor for AOI greater than approximately
80 degrees `and` plane of array irradiance greater than
approximately 200 W/m^2.
Parameters
----------
poa_global : array-like
Plane of array global irradiance in W/m^2.
aoi : array-like
Angle of incidence of solar rays with respect to the module
surface normal.
solar_zenith : array-like
True (not refraction-corrected) solar zenith angles in decimal
degrees.
solar_azimuth : array-like
Solar azimuth angles in decimal degrees.
times : DatetimeIndex
Time indices for the input array-like data.
surface_tilt : numeric
Surface tilt angles in decimal degrees. Tilt must be >=0 and
<=180. The tilt angle is defined as degrees from horizontal
(e.g. surface facing up = 0, surface facing horizon = 90).
surface_azimuth : numeric
Surface azimuth angles in decimal degrees. surface_azimuth must
be >=0 and <=360. The Azimuth convention is defined as degrees
east of north (e.g. North = 0, South=180 East = 90, West = 270).
pressure : numeric, default 101325.0
The site pressure in Pascal. Pressure may be measured or an
average pressure may be calculated from site altitude.
use_delta_kt_prime : bool, default True
If True, indicates that the stability index delta_kt_prime is
included in the model. The stability index adjusts the estimated
DNI in response to dynamics in the time series of GHI. It is
recommended that delta_kt_prime is not used if the time between
GHI points is 1.5 hours or greater. If use_delta_kt_prime=True,
input data must be Series.
temp_dew : float, or array-like, optional
Surface dew point temperatures, in degrees C. Values of temp_dew
may be numeric or NaN. Any single time period point with a
temp_dew=NaN does not have dew point improvements applied. If
temp_dew is not provided, then dew point improvements are not
applied.
albedo : numeric, default 0.25
Ground surface albedo. [unitless]
model : String, default 'perez'
Irradiance model. See :py:func:`get_sky_diffuse` for allowed values.
model_perez : String, default 'allsitescomposite1990'
Used only if model='perez'. See :py:func:`perez`.
calculate_gt_90 : bool, default True
Controls if the algorithm evaluates inputs with AOI >= 90 degrees.
If False, returns nan for AOI >= 90 degrees. Significant speed ups
can be achieved by setting this parameter to False.
max_iterations : int, default 30
Maximum number of iterations for the aoi < 90 deg algorithm.
Returns
-------
data : DataFrame
Contains the following keys/columns:
* ``ghi``: the modeled global horizontal irradiance in W/m^2.
* ``dni``: the modeled direct normal irradiance in W/m^2.
* ``dhi``: the modeled diffuse horizontal irradiance in
W/m^2.
References
----------
.. [1] B. Marion, A model for deriving the direct normal and
diffuse horizontal irradiance from the global tilted
irradiance, Solar Energy 122, 1037-1046.
:doi:`10.1016/j.solener.2015.10.024`
"""
aoi_lt_90 = aoi < 90
# for AOI less than 90 degrees
ghi, dni, dhi, kt_prime = _gti_dirint_lt_90(
poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth, times,
surface_tilt, surface_azimuth, pressure=pressure,
use_delta_kt_prime=use_delta_kt_prime, temp_dew=temp_dew,
albedo=albedo, model=model, model_perez=model_perez,
max_iterations=max_iterations)
# for AOI greater than or equal to 90 degrees
if calculate_gt_90:
ghi_gte_90, dni_gte_90, dhi_gte_90 = _gti_dirint_gte_90(
poa_global, aoi, solar_zenith, solar_azimuth,
surface_tilt, times, kt_prime,
pressure=pressure, temp_dew=temp_dew, albedo=albedo)
else:
ghi_gte_90, dni_gte_90, dhi_gte_90 = np.nan, np.nan, np.nan
# put the AOI < 90 and AOI >= 90 conditions together
output = OrderedDict()
output['ghi'] = ghi.where(aoi_lt_90, ghi_gte_90)
output['dni'] = dni.where(aoi_lt_90, dni_gte_90)
output['dhi'] = dhi.where(aoi_lt_90, dhi_gte_90)
output = pd.DataFrame(output, index=times)
return output
def _gti_dirint_lt_90(poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth,
times, surface_tilt, surface_azimuth, pressure=101325.,
use_delta_kt_prime=True, temp_dew=None, albedo=.25,
model='perez', model_perez='allsitescomposite1990',
max_iterations=30):
"""
GTI-DIRINT model for AOI < 90 degrees. See Marion 2015 Section 2.1.
See gti_dirint signature for parameter details.
"""
I0 = get_extra_radiation(times, 1370, 'spencer')
cos_zenith = tools.cosd(solar_zenith)
# I0h as in Marion 2015 eqns 1, 3
I0h = I0 * np.maximum(0.065, cos_zenith)
airmass = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
airmass = atmosphere.get_absolute_airmass(airmass, pressure)
# these coeffs and diff variables and the loop below
# implement figure 1 of Marion 2015
# make coeffs that is at least 30 elements long so that all
# coeffs can be assigned as specified in Marion 2015.
# slice below will limit iterations if necessary
coeffs = np.empty(max(30, max_iterations))
coeffs[0:3] = 1
coeffs[3:10] = 0.5
coeffs[10:20] = 0.25
coeffs[20:] = 0.125
coeffs = coeffs[:max_iterations] # covers case where max_iterations < 30
# initialize diff
diff = pd.Series(9999, index=times)
best_diff = diff
# initialize poa_global_i
poa_global_i = poa_global
for iteration, coeff in enumerate(coeffs):
# test if difference between modeled GTI and
# measured GTI (poa_global) is less than 1 W/m^2
# only test for aoi less than 90 deg
best_diff_lte_1 = best_diff <= 1
best_diff_lte_1_lt_90 = best_diff_lte_1[aoi_lt_90]
if best_diff_lte_1_lt_90.all():
# all aoi < 90 points have a difference <= 1, so break loop
break
# 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)[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,
use_delta_kt_prime, temp_dew)
# calculate DHI using Marion eqn 3 (identify 1st term on RHS as GHI)
# I0h has a minimum zenith projection, but multiplier of DNI does not
ghi = kt * I0h # Kt * I0 * max(0.065, cos(zen))
dhi = ghi - dni * cos_zenith # no cos(zen) restriction here
# following SSC code
dni = np.maximum(dni, 0)
ghi = np.maximum(ghi, 0)
dhi = np.maximum(dhi, 0)
# use DNI and DHI to model GTI
# GTI-DIRINT uses perez transposition model, but we allow for
# any model here
all_irrad = get_total_irradiance(
surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
dni, ghi, dhi, dni_extra=I0, airmass=airmass,
albedo=albedo, model=model, model_perez=model_perez)
gti_model = all_irrad['poa_global']
# calculate new diff
diff = gti_model - poa_global
# determine if the new diff is smaller in magnitude
# than the old diff
diff_abs = diff.abs()
smallest_diff = diff_abs < best_diff
# save the best differences
best_diff = diff_abs.where(smallest_diff, best_diff)
# on first iteration, the best values are the only values
if iteration == 0:
best_ghi = ghi
best_dni = dni
best_dhi = dhi
best_kt_prime = kt_prime
else:
# save new DNI, DHI, DHI if they provide the best consistency
# otherwise use the older values.
best_ghi = ghi.where(smallest_diff, best_ghi)
best_dni = dni.where(smallest_diff, best_dni)
best_dhi = dhi.where(smallest_diff, best_dhi)
best_kt_prime = kt_prime.where(smallest_diff, best_kt_prime)
# calculate adjusted inputs for next iteration. Marion eqn 4
poa_global_i = np.maximum(1.0, poa_global_i - coeff * diff)
else:
# we are here because we ran out of coeffs to loop over and
# therefore we have exceeded max_iterations
import warnings
failed_points = best_diff[aoi_lt_90][~best_diff_lte_1_lt_90]
warnings.warn(
('%s points failed to converge after %s iterations. best_diff:\n%s'
% (len(failed_points), max_iterations, failed_points)),
RuntimeWarning)
# return the best data, whether or not the solution converged
return best_ghi, best_dni, best_dhi, best_kt_prime
def _gti_dirint_gte_90(poa_global, aoi, solar_zenith, solar_azimuth,
surface_tilt, times, kt_prime,
pressure=101325., temp_dew=None, albedo=.25):
"""
GTI-DIRINT model for AOI >= 90 degrees. See Marion 2015 Section 2.2.
See gti_dirint signature for parameter details.
"""
kt_prime_gte_90 = _gti_dirint_gte_90_kt_prime(aoi, solar_zenith,
solar_azimuth, times,
kt_prime)
I0 = get_extra_radiation(times, 1370, 'spencer')
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)[0] * I0, 0)
dni_gte_90 = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith,
False, temp_dew)
dni_gte_90_proj = dni_gte_90 * tools.cosd(solar_zenith)
cos_surface_tilt = tools.cosd(surface_tilt)
# isotropic sky plus ground diffuse
dhi_gte_90 = (
(2 * poa_global - dni_gte_90_proj * albedo * (1 - cos_surface_tilt)) /
(1 + cos_surface_tilt + albedo * (1 - cos_surface_tilt)))
ghi_gte_90 = dni_gte_90_proj + dhi_gte_90
return ghi_gte_90, dni_gte_90, dhi_gte_90
def _gti_dirint_gte_90_kt_prime(aoi, solar_zenith, solar_azimuth, times,
kt_prime):
"""
Determine kt' values to be used in GTI-DIRINT AOI >= 90 deg case.
See Marion 2015 Section 2.2.
For AOI >= 90 deg: average of the kt_prime values for 65 < AOI < 80
in each day's morning and afternoon. Morning and afternoon are treated
separately.
For AOI < 90 deg: NaN.
See gti_dirint signature for parameter details.
Returns
-------
kt_prime_gte_90 : Series
Index is `times`.
"""
# kt_prime values from DIRINT calculation for AOI < 90 case
# set the kt_prime from sunrise to AOI=90 to be equal to
# the kt_prime for 65 < AOI < 80 during the morning.
# similar for the afternoon. repeat for every day.
aoi_gte_90 = aoi >= 90
aoi_65_80 = (aoi > 65) & (aoi < 80)
zenith_lt_90 = solar_zenith < 90
morning = solar_azimuth < 180
afternoon = solar_azimuth > 180
aoi_65_80_morning = aoi_65_80 & morning
aoi_65_80_afternoon = aoi_65_80 & afternoon
zenith_lt_90_aoi_gte_90_morning = zenith_lt_90 & aoi_gte_90 & morning
zenith_lt_90_aoi_gte_90_afternoon = zenith_lt_90 & aoi_gte_90 & afternoon
kt_prime_gte_90 = []
for date, data in kt_prime.groupby(times.date):
kt_prime_am_avg = data[aoi_65_80_morning].mean()
kt_prime_pm_avg = data[aoi_65_80_afternoon].mean()
kt_prime_by_date = pd.Series(np.nan, index=data.index)
kt_prime_by_date[zenith_lt_90_aoi_gte_90_morning] = kt_prime_am_avg
kt_prime_by_date[zenith_lt_90_aoi_gte_90_afternoon] = kt_prime_pm_avg
kt_prime_gte_90.append(kt_prime_by_date)
kt_prime_gte_90 = pd.concat(kt_prime_gte_90)
return kt_prime_gte_90
[docs]def erbs(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
r"""
Estimate DNI and DHI from GHI using the Erbs model.
The Erbs model [1]_ estimates the diffuse fraction DF from global
horizontal irradiance through an empirical relationship between DF
and the ratio of GHI to extraterrestrial irradiance, Kt. The
function uses the diffuse fraction to compute DHI as
.. math::
DHI = DF \times GHI
DNI is then estimated as
.. math::
DNI = (GHI - DHI)/\cos(Z)
where Z is the zenith angle.
Parameters
----------
ghi: numeric
Global horizontal irradiance in W/m^2.
zenith: numeric
True (not refraction-corrected) zenith angles in decimal degrees.
datetime_or_doy : int, float, array, pd.DatetimeIndex
Day of year or array of days of year e.g.
pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
min_cos_zenith : numeric, default 0.065
Minimum value of cos(zenith) to allow when calculating global
clearness index `kt`. Equivalent to zenith = 86.273 degrees.
max_zenith : numeric, default 87
Maximum value of zenith to allow in DNI calculation. DNI will be
set to 0 for times with zenith values greater than `max_zenith`.
Returns
-------
data : OrderedDict or DataFrame
Contains the following keys/columns:
* ``dni``: the modeled direct normal irradiance in W/m^2.
* ``dhi``: the modeled diffuse horizontal irradiance in
W/m^2.
* ``kt``: Ratio of global to extraterrestrial irradiance
on a horizontal plane.
References
----------
.. [1] D. G. Erbs, S. A. Klein and J. A. Duffie, Estimation of the
diffuse radiation fraction for hourly, daily and monthly-average
global radiation, Solar Energy 28(4), pp 293-302, 1982. Eq. 1
See also
--------
dirint
disc
orgill_hollands
boland
"""
dni_extra = get_extra_radiation(datetime_or_doy)
kt = clearness_index(ghi, zenith, dni_extra, min_cos_zenith=min_cos_zenith,
max_clearness_index=1)
# For Kt <= 0.22, set the diffuse fraction
df = 1 - 0.09*kt
# For Kt > 0.22 and Kt <= 0.8, set the diffuse fraction
df = np.where((kt > 0.22) & (kt <= 0.8),
0.9511 - 0.1604*kt + 4.388*kt**2 -
16.638*kt**3 + 12.336*kt**4,
df)
# For Kt > 0.8, set the diffuse fraction
df = np.where(kt > 0.8, 0.165, df)
dhi = df * ghi
dni = (ghi - dhi) / tools.cosd(zenith)
bad_values = (zenith > max_zenith) | (ghi < 0) | (dni < 0)
dni = np.where(bad_values, 0, dni)
# ensure that closure relationship remains valid
dhi = np.where(bad_values, ghi, dhi)
data = OrderedDict()
data['dni'] = dni
data['dhi'] = dhi
data['kt'] = kt
if isinstance(datetime_or_doy, pd.DatetimeIndex):
data = pd.DataFrame(data, index=datetime_or_doy)
return data
[docs]def erbs_driesse(ghi, zenith, datetime_or_doy=None, dni_extra=None,
min_cos_zenith=0.065, max_zenith=87):
r"""
Estimate DNI and DHI from GHI using the continuous Erbs-Driesse model.
The Erbs-Driesse model [1]_ is a reformulation of the original Erbs
model [2]_ that provides continuity of the function and its first
derivative at the two transition points.
.. math::
DHI = DF \times GHI
DNI is then estimated as
.. math::
DNI = (GHI - DHI)/\cos(Z)
where Z is the zenith angle.
Parameters
----------
ghi: numeric
Global horizontal irradiance in W/m^2.
zenith: numeric
True (not refraction-corrected) zenith angles in decimal degrees.
datetime_or_doy : int, float, array or pd.DatetimeIndex, optional
Day of year or array of days of year e.g.
pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
Either datetime_or_doy or dni_extra must be provided.
dni_extra : numeric, optional
Extraterrestrial normal irradiance.
dni_extra can be provided if available to avoid recalculating it
inside this function. In this case datetime_or_doy is not required.
min_cos_zenith : numeric, default 0.065
Minimum value of cos(zenith) to allow when calculating global
clearness index `kt`. Equivalent to zenith = 86.273 degrees.
max_zenith : numeric, default 87
Maximum value of zenith to allow in DNI calculation. DNI will be
set to 0 for times with zenith values greater than `max_zenith`.
Returns
-------
data : OrderedDict or DataFrame
Contains the following keys/columns:
* ``dni``: the modeled direct normal irradiance in W/m^2.
* ``dhi``: the modeled diffuse horizontal irradiance in
W/m^2.
* ``kt``: Ratio of global to extraterrestrial irradiance
on a horizontal plane.
Raises
------
ValueError
If neither datetime_or_doy nor dni_extra is provided.
Notes
-----
The diffuse fraction DHI/GHI of the Erbs-Driesse model deviates from the
original Erbs model by less than 0.0005.
References
----------
.. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the
Perez diffuse sky model for forward and reverse transposition.
Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093`
.. [2] D. G. Erbs, S. A. Klein and J. A. Duffie, Estimation of the
diffuse radiation fraction for hourly, daily and monthly-average
global radiation, Solar Energy 28(4), pp 293-302, 1982. Eq. 1
See also
--------
erbs
dirint
disc
orgill_hollands
boland
"""
# Contributed by Anton Driesse (@adriesse), PV Performance Labs. Aug., 2023
# central polynomial coefficients with float64 precision
p = [+12.26911439571261000,
-16.47050842469730700,
+04.24692671521831700,
-00.11390583806313881,
+00.94629663357100100]
if datetime_or_doy is None and dni_extra is None:
raise ValueError('Either datetime_or_doy or dni_extra '
'must be provided.')
if dni_extra is None:
dni_extra = get_extra_radiation(datetime_or_doy)
# negative ghi should not reach this point, but just in case
ghi = np.maximum(0, ghi)
kt = clearness_index(ghi, zenith, dni_extra, min_cos_zenith=min_cos_zenith,
max_clearness_index=1)
# For all Kt, set the default diffuse fraction
df = 1 - 0.09 * kt
# For Kt > 0.216, update the diffuse fraction
df = np.where(kt > 0.216, np.polyval(p, kt), df)
# For Kt > 0.792, update the diffuse fraction again
df = np.where(kt > 0.792, 0.165, df)
dhi = df * ghi
dni = (ghi - dhi) / tools.cosd(zenith)
bad_values = (zenith > max_zenith) | (ghi < 0) | (dni < 0)
dni = np.where(bad_values, 0, dni)
# ensure that closure relationship remains valid
dhi = np.where(bad_values, ghi, dhi)
data = OrderedDict()
data['dni'] = dni
data['dhi'] = dhi
data['kt'] = kt
if isinstance(datetime_or_doy, pd.DatetimeIndex):
data = pd.DataFrame(data, index=datetime_or_doy)
elif isinstance(ghi, pd.Series):
data = pd.DataFrame(data, index=ghi.index)
return data
[docs]def orgill_hollands(ghi, zenith, datetime_or_doy, dni_extra=None,
min_cos_zenith=0.065, max_zenith=87):
"""Estimate DNI and DHI from GHI using the Orgill and Hollands model.
The Orgill and Hollands model [1]_ estimates the diffuse fraction DF from
global horizontal irradiance through an empirical relationship between
hourly DF observations (in Toronto, Canada) and the ratio of GHI to
extraterrestrial irradiance, Kt.
Parameters
----------
ghi: numeric
Global horizontal irradiance in W/m^2.
zenith: numeric
True (not refraction-corrected) zenith angles in decimal degrees.
datetime_or_doy : int, float, array, pd.DatetimeIndex
Day of year or array of days of year e.g.
pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
dni_extra : numeric, optional
Extraterrestrial direct normal irradiance. [W/m2]
min_cos_zenith : numeric, default 0.065
Minimum value of cos(zenith) to allow when calculating global
clearness index `kt`. Equivalent to zenith = 86.273 degrees.
max_zenith : numeric, default 87
Maximum value of zenith to allow in DNI calculation. DNI will be
set to 0 for times with zenith values greater than `max_zenith`.
Returns
-------
data : OrderedDict or DataFrame
Contains the following keys/columns:
* ``dni``: the modeled direct normal irradiance in W/m^2.
* ``dhi``: the modeled diffuse horizontal irradiance in
W/m^2.
* ``kt``: Ratio of global to extraterrestrial irradiance
on a horizontal plane.
References
----------
.. [1] Orgill, J.F., Hollands, K.G.T., Correlation equation for hourly
diffuse radiation on a horizontal surface, Solar Energy 19(4), pp 357–359,
1977. Eqs. 3(a), 3(b) and 3(c)
:doi:`10.1016/0038-092X(77)90006-8`
See Also
--------
dirint
disc
erbs
boland
"""
if dni_extra is None:
dni_extra = get_extra_radiation(datetime_or_doy)
kt = clearness_index(ghi, zenith, dni_extra, min_cos_zenith=min_cos_zenith,
max_clearness_index=1)
# For Kt < 0.35, set the diffuse fraction
df = 1 - 0.249*kt
# For Kt >= 0.35 and Kt <= 0.75, set the diffuse fraction
df = np.where((kt >= 0.35) & (kt <= 0.75),
1.557 - 1.84*kt, df)
# For Kt > 0.75, set the diffuse fraction
df = np.where(kt > 0.75, 0.177, df)
dhi = df * ghi
dni = (ghi - dhi) / tools.cosd(zenith)
bad_values = (zenith > max_zenith) | (ghi < 0) | (dni < 0)
dni = np.where(bad_values, 0, dni)
# ensure that closure relationship remains valid
dhi = np.where(bad_values, ghi, dhi)
data = OrderedDict()
data['dni'] = dni
data['dhi'] = dhi
data['kt'] = kt
if isinstance(datetime_or_doy, pd.DatetimeIndex):
data = pd.DataFrame(data, index=datetime_or_doy)
return data
[docs]def boland(ghi, solar_zenith, datetime_or_doy, a_coeff=8.645, b_coeff=0.613,
min_cos_zenith=0.065, max_zenith=87):
r"""
Estimate DNI and DHI from GHI using the Boland clearness index model.
The Boland model [1]_, [2]_ estimates the diffuse fraction, DF, from global
horizontal irradiance, GHI, through an empirical relationship between DF
and the clearness index, :math:`k_t`, the ratio of GHI to horizontal
extraterrestrial irradiance.
.. math::
\mathit{DF} = \frac{1}{1 + \exp\left(a \left(k_t - b\right)\right)}
Parameters
----------
ghi: numeric
Global horizontal irradiance. [W/m^2]
solar_zenith: numeric
True (not refraction-corrected) zenith angles in decimal degrees.
datetime_or_doy : numeric, pandas.DatetimeIndex
Day of year or array of days of year e.g.
pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
a_coeff : float, default 8.645
Logistic curve fit coefficient.
b_coeff : float, default 0.613
Logistic curve fit coefficient.
min_cos_zenith : numeric, default 0.065
Minimum value of cos(zenith) to allow when calculating global
clearness index :math:`k_t`. Equivalent to zenith = 86.273 degrees.
max_zenith : numeric, default 87
Maximum value of zenith to allow in DNI calculation. DNI will be
set to 0 for times with zenith values greater than `max_zenith`.
Returns
-------
data : OrderedDict or DataFrame
Contains the following keys/columns:
* ``dni``: the modeled direct normal irradiance in W/m^2.
* ``dhi``: the modeled diffuse horizontal irradiance in
W/m^2.
* ``kt``: Ratio of global to extraterrestrial irradiance
on a horizontal plane.
References
----------
.. [1] J. Boland, B. Ridley (2008) Models of Diffuse Solar Fraction. In:
Badescu V. (eds) Modeling Solar Radiation at the Earth’s Surface.
Springer, Berlin, Heidelberg. :doi:`10.1007/978-3-540-77455-6_8`
.. [2] John Boland, Lynne Scott, and Mark Luther, Modelling the diffuse
fraction of global solar radiation on a horizontal surface,
Environmetrics 12(2), pp 103-116, 2001,
:doi:`10.1002/1099-095X(200103)12:2%3C103::AID-ENV447%3E3.0.CO;2-2`
See also
--------
dirint
disc
erbs
orgill_hollands
Notes
-----
Boland diffuse fraction differs from other decomposition algorithms by use
of a logistic function to fit the entire range of clearness index,
:math:`k_t`. Parameters ``a_coeff`` and ``b_coeff`` are reported in [2]_
for different time intervals:
* 15-minute: ``a = 8.645`` and ``b = 0.613``
* 1-hour: ``a = 7.997`` and ``b = 0.586``
"""
dni_extra = get_extra_radiation(datetime_or_doy)
kt = clearness_index(
ghi, solar_zenith, dni_extra, min_cos_zenith=min_cos_zenith,
max_clearness_index=1)
# Boland equation
df = 1.0 / (1.0 + np.exp(a_coeff * (kt - b_coeff)))
# NOTE: [2] has different coefficients, for different time intervals
# 15-min: df = 1 / (1 + exp(8.645 * (kt - 0.613)))
# 1-hour: df = 1 / (1 + exp(7.997 * (kt - 0.586)))
dhi = df * ghi
dni = (ghi - dhi) / tools.cosd(solar_zenith)
bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0)
dni = np.where(bad_values, 0, dni)
# ensure that closure relationship remains valid
dhi = np.where(bad_values, ghi, dhi)
data = OrderedDict()
data['dni'] = dni
data['dhi'] = dhi
data['kt'] = kt
if isinstance(datetime_or_doy, pd.DatetimeIndex):
data = pd.DataFrame(data, index=datetime_or_doy)
return data
[docs]def campbell_norman(zenith, transmittance, pressure=101325.0,
dni_extra=1367.0):
'''
Determine DNI, DHI, GHI from extraterrestrial flux, transmittance,
and atmospheric pressure.
Parameters
----------
zenith: pd.Series
True (not refraction-corrected) zenith angles in decimal
degrees. If Z is a vector it must be of the same size as all
other vector inputs. Z must be >=0 and <=180.
transmittance: float
Atmospheric transmittance between 0 and 1.
pressure: float, default 101325.0
Air pressure
dni_extra: float, default 1367.0
Direct irradiance incident at the top of the atmosphere.
Returns
-------
irradiance: DataFrame
Modeled direct normal irradiance, direct horizontal irradiance,
and global horizontal irradiance in W/m^2
References
----------
.. [1] Campbell, G. S., J. M. Norman (1998) An Introduction to
Environmental Biophysics. 2nd Ed. New York: Springer.
'''
tau = transmittance
airmass = atmosphere.get_relative_airmass(zenith, model='simple')
airmass = atmosphere.get_absolute_airmass(airmass, pressure=pressure)
dni = dni_extra*tau**airmass
cos_zen = tools.cosd(zenith)
dhi = 0.3 * (1.0 - tau**airmass) * dni_extra * cos_zen
ghi = dhi + dni * cos_zen
irrads = OrderedDict()
irrads['ghi'] = ghi
irrads['dni'] = dni
irrads['dhi'] = dhi
if isinstance(ghi, pd.Series):
irrads = pd.DataFrame(irrads)
return irrads
def _liujordan(zenith, transmittance, airmass, dni_extra=1367.0):
'''
Determine DNI, DHI, GHI from extraterrestrial flux, transmittance,
and optical air mass number.
Liu and Jordan, 1960, developed a simplified direct radiation model.
DHI is from an empirical equation for diffuse radiation from Liu and
Jordan, 1960.
Parameters
----------
zenith: pd.Series
True (not refraction-corrected) zenith angles in decimal
degrees. If Z is a vector it must be of the same size as all
other vector inputs. Z must be >=0 and <=180.
transmittance: float
Atmospheric transmittance between 0 and 1.
pressure: float, default 101325.0
Air pressure
dni_extra: float, default 1367.0
Direct irradiance incident at the top of the atmosphere.
Returns
-------
irradiance: DataFrame
Modeled direct normal irradiance, direct horizontal irradiance,
and global horizontal irradiance in W/m^2
References
----------
.. [1] Campbell, G. S., J. M. Norman (1998) An Introduction to
Environmental Biophysics. 2nd Ed. New York: Springer.
.. [2] Liu, B. Y., R. C. Jordan, (1960). "The interrelationship and
characteristic distribution of direct, diffuse, and total solar
radiation". Solar Energy 4:1-19
'''
tau = transmittance
dni = dni_extra*tau**airmass
dhi = 0.3 * (1.0 - tau**airmass) * dni_extra * np.cos(np.radians(zenith))
ghi = dhi + dni * np.cos(np.radians(zenith))
irrads = OrderedDict()
irrads['ghi'] = ghi
irrads['dni'] = dni
irrads['dhi'] = dhi
if isinstance(ghi, pd.Series):
irrads = pd.DataFrame(irrads)
return irrads
def _get_perez_coefficients(perezmodel):
'''
Find coefficients for the Perez model
Parameters
----------
perezmodel : string, default 'allsitescomposite1990'
a character string which selects the desired set of Perez
coefficients. If model is not provided as an input, the default,
'1990' will be used.
All possible model selections are:
* '1990'
* 'allsitescomposite1990' (same as '1990')
* 'allsitescomposite1988'
* 'sandiacomposite1988'
* 'usacomposite1988'
* 'france1988'
* 'phoenix1988'
* 'elmonte1988'
* 'osage1988'
* 'albuquerque1988'
* 'capecanaveral1988'
* 'albany1988'
Returns
--------
F1coeffs, F2coeffs : (array, array)
F1 and F2 coefficients for the Perez model
References
----------
.. [1] Loutzenhiser P.G. et. al. "Empirical validation of models to
compute solar irradiance on inclined surfaces for building energy
simulation" 2007, Solar Energy vol. 81. pp. 254-267
.. [2] Perez, R., Seals, R., Ineichen, P., Stewart, R., Menicucci, D.,
1987. A new simplified version of the Perez diffuse irradiance model
for tilted surfaces. Solar Energy 39(3), 221-232.
.. [3] Perez, R., Ineichen, P., Seals, R., Michalsky, J., Stewart, R.,
1990. Modeling daylight availability and irradiance components from
direct and global irradiance. Solar Energy 44 (5), 271-289.
.. [4] Perez, R. et. al 1988. "The Development and Verification of the
Perez Diffuse Radiation Model". SAND88-7030
'''
coeffdict = {
'allsitescomposite1990': [
[-0.0080, 0.5880, -0.0620, -0.0600, 0.0720, -0.0220],
[0.1300, 0.6830, -0.1510, -0.0190, 0.0660, -0.0290],
[0.3300, 0.4870, -0.2210, 0.0550, -0.0640, -0.0260],
[0.5680, 0.1870, -0.2950, 0.1090, -0.1520, -0.0140],
[0.8730, -0.3920, -0.3620, 0.2260, -0.4620, 0.0010],
[1.1320, -1.2370, -0.4120, 0.2880, -0.8230, 0.0560],
[1.0600, -1.6000, -0.3590, 0.2640, -1.1270, 0.1310],
[0.6780, -0.3270, -0.2500, 0.1560, -1.3770, 0.2510]],
'allsitescomposite1988': [
[-0.0180, 0.7050, -0.071, -0.0580, 0.1020, -0.0260],
[0.1910, 0.6450, -0.1710, 0.0120, 0.0090, -0.0270],
[0.4400, 0.3780, -0.2560, 0.0870, -0.1040, -0.0250],
[0.7560, -0.1210, -0.3460, 0.1790, -0.3210, -0.0080],
[0.9960, -0.6450, -0.4050, 0.2600, -0.5900, 0.0170],
[1.0980, -1.2900, -0.3930, 0.2690, -0.8320, 0.0750],
[0.9730, -1.1350, -0.3780, 0.1240, -0.2580, 0.1490],
[0.6890, -0.4120, -0.2730, 0.1990, -1.6750, 0.2370]],
'sandiacomposite1988': [
[-0.1960, 1.0840, -0.0060, -0.1140, 0.1800, -0.0190],
[0.2360, 0.5190, -0.1800, -0.0110, 0.0200, -0.0380],
[0.4540, 0.3210, -0.2550, 0.0720, -0.0980, -0.0460],
[0.8660, -0.3810, -0.3750, 0.2030, -0.4030, -0.0490],
[1.0260, -0.7110, -0.4260, 0.2730, -0.6020, -0.0610],
[0.9780, -0.9860, -0.3500, 0.2800, -0.9150, -0.0240],
[0.7480, -0.9130, -0.2360, 0.1730, -1.0450, 0.0650],
[0.3180, -0.7570, 0.1030, 0.0620, -1.6980, 0.2360]],
'usacomposite1988': [
[-0.0340, 0.6710, -0.0590, -0.0590, 0.0860, -0.0280],
[0.2550, 0.4740, -0.1910, 0.0180, -0.0140, -0.0330],
[0.4270, 0.3490, -0.2450, 0.0930, -0.1210, -0.0390],
[0.7560, -0.2130, -0.3280, 0.1750, -0.3040, -0.0270],
[1.0200, -0.8570, -0.3850, 0.2800, -0.6380, -0.0190],
[1.0500, -1.3440, -0.3480, 0.2800, -0.8930, 0.0370],
[0.9740, -1.5070, -0.3700, 0.1540, -0.5680, 0.1090],
[0.7440, -1.8170, -0.2560, 0.2460, -2.6180, 0.2300]],
'france1988': [
[0.0130, 0.7640, -0.1000, -0.0580, 0.1270, -0.0230],
[0.0950, 0.9200, -0.1520, 0, 0.0510, -0.0200],
[0.4640, 0.4210, -0.2800, 0.0640, -0.0510, -0.0020],
[0.7590, -0.0090, -0.3730, 0.2010, -0.3820, 0.0100],
[0.9760, -0.4000, -0.4360, 0.2710, -0.6380, 0.0510],
[1.1760, -1.2540, -0.4620, 0.2950, -0.9750, 0.1290],
[1.1060, -1.5630, -0.3980, 0.3010, -1.4420, 0.2120],
[0.9340, -1.5010, -0.2710, 0.4200, -2.9170, 0.2490]],
'phoenix1988': [
[-0.0030, 0.7280, -0.0970, -0.0750, 0.1420, -0.0430],
[0.2790, 0.3540, -0.1760, 0.0300, -0.0550, -0.0540],
[0.4690, 0.1680, -0.2460, 0.0480, -0.0420, -0.0570],
[0.8560, -0.5190, -0.3400, 0.1760, -0.3800, -0.0310],
[0.9410, -0.6250, -0.3910, 0.1880, -0.3600, -0.0490],
[1.0560, -1.1340, -0.4100, 0.2810, -0.7940, -0.0650],
[0.9010, -2.1390, -0.2690, 0.1180, -0.6650, 0.0460],
[0.1070, 0.4810, 0.1430, -0.1110, -0.1370, 0.2340]],
'elmonte1988': [
[0.0270, 0.7010, -0.1190, -0.0580, 0.1070, -0.0600],
[0.1810, 0.6710, -0.1780, -0.0790, 0.1940, -0.0350],
[0.4760, 0.4070, -0.2880, 0.0540, -0.0320, -0.0550],
[0.8750, -0.2180, -0.4030, 0.1870, -0.3090, -0.0610],
[1.1660, -1.0140, -0.4540, 0.2110, -0.4100, -0.0440],
[1.1430, -2.0640, -0.2910, 0.0970, -0.3190, 0.0530],
[1.0940, -2.6320, -0.2590, 0.0290, -0.4220, 0.1470],
[0.1550, 1.7230, 0.1630, -0.1310, -0.0190, 0.2770]],
'osage1988': [
[-0.3530, 1.4740, 0.0570, -0.1750, 0.3120, 0.0090],
[0.3630, 0.2180, -0.2120, 0.0190, -0.0340, -0.0590],
[-0.0310, 1.2620, -0.0840, -0.0820, 0.2310, -0.0170],
[0.6910, 0.0390, -0.2950, 0.0910, -0.1310, -0.0350],
[1.1820, -1.3500, -0.3210, 0.4080, -0.9850, -0.0880],
[0.7640, 0.0190, -0.2030, 0.2170, -0.2940, -0.1030],
[0.2190, 1.4120, 0.2440, 0.4710, -2.9880, 0.0340],
[3.5780, 22.2310, -10.7450, 2.4260, 4.8920, -5.6870]],
'albuquerque1988': [
[0.0340, 0.5010, -0.0940, -0.0630, 0.1060, -0.0440],
[0.2290, 0.4670, -0.1560, -0.0050, -0.0190, -0.0230],
[0.4860, 0.2410, -0.2530, 0.0530, -0.0640, -0.0220],
[0.8740, -0.3930, -0.3970, 0.1810, -0.3270, -0.0370],
[1.1930, -1.2960, -0.5010, 0.2810, -0.6560, -0.0450],
[1.0560, -1.7580, -0.3740, 0.2260, -0.7590, 0.0340],
[0.9010, -4.7830, -0.1090, 0.0630, -0.9700, 0.1960],
[0.8510, -7.0550, -0.0530, 0.0600, -2.8330, 0.3300]],
'capecanaveral1988': [
[0.0750, 0.5330, -0.1240, -0.0670, 0.0420, -0.0200],
[0.2950, 0.4970, -0.2180, -0.0080, 0.0030, -0.0290],
[0.5140, 0.0810, -0.2610, 0.0750, -0.1600, -0.0290],
[0.7470, -0.3290, -0.3250, 0.1810, -0.4160, -0.0300],
[0.9010, -0.8830, -0.2970, 0.1780, -0.4890, 0.0080],
[0.5910, -0.0440, -0.1160, 0.2350, -0.9990, 0.0980],
[0.5370, -2.4020, 0.3200, 0.1690, -1.9710, 0.3100],
[-0.8050, 4.5460, 1.0720, -0.2580, -0.9500, 0.7530]],
'albany1988': [
[0.0120, 0.5540, -0.0760, -0.0520, 0.0840, -0.0290],
[0.2670, 0.4370, -0.1940, 0.0160, 0.0220, -0.0360],
[0.4200, 0.3360, -0.2370, 0.0740, -0.0520, -0.0320],
[0.6380, -0.0010, -0.2810, 0.1380, -0.1890, -0.0120],
[1.0190, -1.0270, -0.3420, 0.2710, -0.6280, 0.0140],
[1.1490, -1.9400, -0.3310, 0.3220, -1.0970, 0.0800],
[1.4340, -3.9940, -0.4920, 0.4530, -2.3760, 0.1170],
[1.0070, -2.2920, -0.4820, 0.3900, -3.3680, 0.2290]], }
array = np.array(coeffdict[perezmodel])
F1coeffs = array[:, 0:3]
F2coeffs = array[:, 3:7]
return F1coeffs, F2coeffs
def _get_dirint_coeffs():
"""
A place to stash the dirint coefficients.
Returns
-------
np.array with shape ``(6, 6, 7, 5)``.
Ordering is ``[kt_prime_bin, zenith_bin, delta_kt_prime_bin, w_bin]``
"""
# To allow for maximum copy/paste from the MATLAB 1-indexed code,
# we create and assign values to an oversized array.
# Then, we return the [1:, 1:, :, :] slice.
coeffs = np.zeros((7, 7, 7, 5))
coeffs[1, 1, :, :] = [
[0.385230, 0.385230, 0.385230, 0.462880, 0.317440],
[0.338390, 0.338390, 0.221270, 0.316730, 0.503650],
[0.235680, 0.235680, 0.241280, 0.157830, 0.269440],
[0.830130, 0.830130, 0.171970, 0.841070, 0.457370],
[0.548010, 0.548010, 0.478000, 0.966880, 1.036370],
[0.548010, 0.548010, 1.000000, 3.012370, 1.976540],
[0.582690, 0.582690, 0.229720, 0.892710, 0.569950]]
coeffs[1, 2, :, :] = [
[0.131280, 0.131280, 0.385460, 0.511070, 0.127940],
[0.223710, 0.223710, 0.193560, 0.304560, 0.193940],
[0.229970, 0.229970, 0.275020, 0.312730, 0.244610],
[0.090100, 0.184580, 0.260500, 0.687480, 0.579440],
[0.131530, 0.131530, 0.370190, 1.380350, 1.052270],
[1.116250, 1.116250, 0.928030, 3.525490, 2.316920],
[0.090100, 0.237000, 0.300040, 0.812470, 0.664970]]
coeffs[1, 3, :, :] = [
[0.587510, 0.130000, 0.400000, 0.537210, 0.832490],
[0.306210, 0.129830, 0.204460, 0.500000, 0.681640],
[0.224020, 0.260620, 0.334080, 0.501040, 0.350470],
[0.421540, 0.753970, 0.750660, 3.706840, 0.983790],
[0.706680, 0.373530, 1.245670, 0.864860, 1.992630],
[4.864400, 0.117390, 0.265180, 0.359180, 3.310820],
[0.392080, 0.493290, 0.651560, 1.932780, 0.898730]]
coeffs[1, 4, :, :] = [
[0.126970, 0.126970, 0.126970, 0.126970, 0.126970],
[0.810820, 0.810820, 0.810820, 0.810820, 0.810820],
[3.241680, 2.500000, 2.291440, 2.291440, 2.291440],
[4.000000, 3.000000, 2.000000, 0.975430, 1.965570],
[12.494170, 12.494170, 8.000000, 5.083520, 8.792390],
[21.744240, 21.744240, 21.744240, 21.744240, 21.744240],
[3.241680, 12.494170, 1.620760, 1.375250, 2.331620]]
coeffs[1, 5, :, :] = [
[0.126970, 0.126970, 0.126970, 0.126970, 0.126970],
[0.810820, 0.810820, 0.810820, 0.810820, 0.810820],
[3.241680, 2.500000, 2.291440, 2.291440, 2.291440],
[4.000000, 3.000000, 2.000000, 0.975430, 1.965570],
[12.494170, 12.494170, 8.000000, 5.083520, 8.792390],
[21.744240, 21.744240, 21.744240, 21.744240, 21.744240],
[3.241680, 12.494170, 1.620760, 1.375250, 2.331620]]
coeffs[1, 6, :, :] = [
[0.126970, 0.126970, 0.126970, 0.126970, 0.126970],
[0.810820, 0.810820, 0.810820, 0.810820, 0.810820],
[3.241680, 2.500000, 2.291440, 2.291440, 2.291440],
[4.000000, 3.000000, 2.000000, 0.975430, 1.965570],
[12.494170, 12.494170, 8.000000, 5.083520, 8.792390],
[21.744240, 21.744240, 21.744240, 21.744240, 21.744240],
[3.241680, 12.494170, 1.620760, 1.375250, 2.331620]]
coeffs[2, 1, :, :] = [
[0.337440, 0.337440, 0.969110, 1.097190, 1.116080],
[0.337440, 0.337440, 0.969110, 1.116030, 0.623900],
[0.337440, 0.337440, 1.530590, 1.024420, 0.908480],
[0.584040, 0.584040, 0.847250, 0.914940, 1.289300],
[0.337440, 0.337440, 0.310240, 1.435020, 1.852830],
[0.337440, 0.337440, 1.015010, 1.097190, 2.117230],
[0.337440, 0.337440, 0.969110, 1.145730, 1.476400]]
coeffs[2, 2, :, :] = [
[0.300000, 0.300000, 0.700000, 1.100000, 0.796940],
[0.219870, 0.219870, 0.526530, 0.809610, 0.649300],
[0.386650, 0.386650, 0.119320, 0.576120, 0.685460],
[0.746730, 0.399830, 0.470970, 0.986530, 0.785370],
[0.575420, 0.936700, 1.649200, 1.495840, 1.335590],
[1.319670, 4.002570, 1.276390, 2.644550, 2.518670],
[0.665190, 0.678910, 1.012360, 1.199940, 0.986580]]
coeffs[2, 3, :, :] = [
[0.378870, 0.974060, 0.500000, 0.491880, 0.665290],
[0.105210, 0.263470, 0.407040, 0.553460, 0.582590],
[0.312900, 0.345240, 1.144180, 0.854790, 0.612280],
[0.119070, 0.365120, 0.560520, 0.793720, 0.802600],
[0.781610, 0.837390, 1.270420, 1.537980, 1.292950],
[1.152290, 1.152290, 1.492080, 1.245370, 2.177100],
[0.424660, 0.529550, 0.966910, 1.033460, 0.958730]]
coeffs[2, 4, :, :] = [
[0.310590, 0.714410, 0.252450, 0.500000, 0.607600],
[0.975190, 0.363420, 0.500000, 0.400000, 0.502800],
[0.175580, 0.196250, 0.476360, 1.072470, 0.490510],
[0.719280, 0.698620, 0.657770, 1.190840, 0.681110],
[0.426240, 1.464840, 0.678550, 1.157730, 0.978430],
[2.501120, 1.789130, 1.387090, 2.394180, 2.394180],
[0.491640, 0.677610, 0.685610, 1.082400, 0.735410]]
coeffs[2, 5, :, :] = [
[0.597000, 0.500000, 0.300000, 0.310050, 0.413510],
[0.314790, 0.336310, 0.400000, 0.400000, 0.442460],
[0.166510, 0.460440, 0.552570, 1.000000, 0.461610],
[0.401020, 0.559110, 0.403630, 1.016710, 0.671490],
[0.400360, 0.750830, 0.842640, 1.802600, 1.023830],
[3.315300, 1.510380, 2.443650, 1.638820, 2.133990],
[0.530790, 0.745850, 0.693050, 1.458040, 0.804500]]
coeffs[2, 6, :, :] = [
[0.597000, 0.500000, 0.300000, 0.310050, 0.800920],
[0.314790, 0.336310, 0.400000, 0.400000, 0.237040],
[0.166510, 0.460440, 0.552570, 1.000000, 0.581990],
[0.401020, 0.559110, 0.403630, 1.016710, 0.898570],
[0.400360, 0.750830, 0.842640, 1.802600, 3.400390],
[3.315300, 1.510380, 2.443650, 1.638820, 2.508780],
[0.204340, 1.157740, 2.003080, 2.622080, 1.409380]]
coeffs[3, 1, :, :] = [
[1.242210, 1.242210, 1.242210, 1.242210, 1.242210],
[0.056980, 0.056980, 0.656990, 0.656990, 0.925160],
[0.089090, 0.089090, 1.040430, 1.232480, 1.205300],
[1.053850, 1.053850, 1.399690, 1.084640, 1.233340],
[1.151540, 1.151540, 1.118290, 1.531640, 1.411840],
[1.494980, 1.494980, 1.700000, 1.800810, 1.671600],
[1.018450, 1.018450, 1.153600, 1.321890, 1.294670]]
coeffs[3, 2, :, :] = [
[0.700000, 0.700000, 1.023460, 0.700000, 0.945830],
[0.886300, 0.886300, 1.333620, 0.800000, 1.066620],
[0.902180, 0.902180, 0.954330, 1.126690, 1.097310],
[1.095300, 1.075060, 1.176490, 1.139470, 1.096110],
[1.201660, 1.201660, 1.438200, 1.256280, 1.198060],
[1.525850, 1.525850, 1.869160, 1.985410, 1.911590],
[1.288220, 1.082810, 1.286370, 1.166170, 1.119330]]
coeffs[3, 3, :, :] = [
[0.600000, 1.029910, 0.859890, 0.550000, 0.813600],
[0.604450, 1.029910, 0.859890, 0.656700, 0.928840],
[0.455850, 0.750580, 0.804930, 0.823000, 0.911000],
[0.526580, 0.932310, 0.908620, 0.983520, 0.988090],
[1.036110, 1.100690, 0.848380, 1.035270, 1.042380],
[1.048440, 1.652720, 0.900000, 2.350410, 1.082950],
[0.817410, 0.976160, 0.861300, 0.974780, 1.004580]]
coeffs[3, 4, :, :] = [
[0.782110, 0.564280, 0.600000, 0.600000, 0.665740],
[0.894480, 0.680730, 0.541990, 0.800000, 0.669140],
[0.487460, 0.818950, 0.841830, 0.872540, 0.709040],
[0.709310, 0.872780, 0.908480, 0.953290, 0.844350],
[0.863920, 0.947770, 0.876220, 1.078750, 0.936910],
[1.280350, 0.866720, 0.769790, 1.078750, 0.975130],
[0.725420, 0.869970, 0.868810, 0.951190, 0.829220]]
coeffs[3, 5, :, :] = [
[0.791750, 0.654040, 0.483170, 0.409000, 0.597180],
[0.566140, 0.948990, 0.971820, 0.653570, 0.718550],
[0.648710, 0.637730, 0.870510, 0.860600, 0.694300],
[0.637630, 0.767610, 0.925670, 0.990310, 0.847670],
[0.736380, 0.946060, 1.117590, 1.029340, 0.947020],
[1.180970, 0.850000, 1.050000, 0.950000, 0.888580],
[0.700560, 0.801440, 0.961970, 0.906140, 0.823880]]
coeffs[3, 6, :, :] = [
[0.500000, 0.500000, 0.586770, 0.470550, 0.629790],
[0.500000, 0.500000, 1.056220, 1.260140, 0.658140],
[0.500000, 0.500000, 0.631830, 0.842620, 0.582780],
[0.554710, 0.734730, 0.985820, 0.915640, 0.898260],
[0.712510, 1.205990, 0.909510, 1.078260, 0.885610],
[1.899260, 1.559710, 1.000000, 1.150000, 1.120390],
[0.653880, 0.793120, 0.903320, 0.944070, 0.796130]]
coeffs[4, 1, :, :] = [
[1.000000, 1.000000, 1.050000, 1.170380, 1.178090],
[0.960580, 0.960580, 1.059530, 1.179030, 1.131690],
[0.871470, 0.871470, 0.995860, 1.141910, 1.114600],
[1.201590, 1.201590, 0.993610, 1.109380, 1.126320],
[1.065010, 1.065010, 0.828660, 0.939970, 1.017930],
[1.065010, 1.065010, 0.623690, 1.119620, 1.132260],
[1.071570, 1.071570, 0.958070, 1.114130, 1.127110]]
coeffs[4, 2, :, :] = [
[0.950000, 0.973390, 0.852520, 1.092200, 1.096590],
[0.804120, 0.913870, 0.980990, 1.094580, 1.042420],
[0.737540, 0.935970, 0.999940, 1.056490, 1.050060],
[1.032980, 1.034540, 0.968460, 1.032080, 1.015780],
[0.900000, 0.977210, 0.945960, 1.008840, 0.969960],
[0.600000, 0.750000, 0.750000, 0.844710, 0.899100],
[0.926800, 0.965030, 0.968520, 1.044910, 1.032310]]
coeffs[4, 3, :, :] = [
[0.850000, 1.029710, 0.961100, 1.055670, 1.009700],
[0.818530, 0.960010, 0.996450, 1.081970, 1.036470],
[0.765380, 0.953500, 0.948260, 1.052110, 1.000140],
[0.775610, 0.909610, 0.927800, 0.987800, 0.952100],
[1.000990, 0.881880, 0.875950, 0.949100, 0.893690],
[0.902370, 0.875960, 0.807990, 0.942410, 0.917920],
[0.856580, 0.928270, 0.946820, 1.032260, 0.972990]]
coeffs[4, 4, :, :] = [
[0.750000, 0.857930, 0.983800, 1.056540, 0.980240],
[0.750000, 0.987010, 1.013730, 1.133780, 1.038250],
[0.800000, 0.947380, 1.012380, 1.091270, 0.999840],
[0.800000, 0.914550, 0.908570, 0.999190, 0.915230],
[0.778540, 0.800590, 0.799070, 0.902180, 0.851560],
[0.680190, 0.317410, 0.507680, 0.388910, 0.646710],
[0.794920, 0.912780, 0.960830, 1.057110, 0.947950]]
coeffs[4, 5, :, :] = [
[0.750000, 0.833890, 0.867530, 1.059890, 0.932840],
[0.979700, 0.971470, 0.995510, 1.068490, 1.030150],
[0.858850, 0.987920, 1.043220, 1.108700, 1.044900],
[0.802400, 0.955110, 0.911660, 1.045070, 0.944470],
[0.884890, 0.766210, 0.885390, 0.859070, 0.818190],
[0.615680, 0.700000, 0.850000, 0.624620, 0.669300],
[0.835570, 0.946150, 0.977090, 1.049350, 0.979970]]
coeffs[4, 6, :, :] = [
[0.689220, 0.809600, 0.900000, 0.789500, 0.853990],
[0.854660, 0.852840, 0.938200, 0.923110, 0.955010],
[0.938600, 0.932980, 1.010390, 1.043950, 1.041640],
[0.843620, 0.981300, 0.951590, 0.946100, 0.966330],
[0.694740, 0.814690, 0.572650, 0.400000, 0.726830],
[0.211370, 0.671780, 0.416340, 0.297290, 0.498050],
[0.843540, 0.882330, 0.911760, 0.898420, 0.960210]]
coeffs[5, 1, :, :] = [
[1.054880, 1.075210, 1.068460, 1.153370, 1.069220],
[1.000000, 1.062220, 1.013470, 1.088170, 1.046200],
[0.885090, 0.993530, 0.942590, 1.054990, 1.012740],
[0.920000, 0.950000, 0.978720, 1.020280, 0.984440],
[0.850000, 0.908500, 0.839940, 0.985570, 0.962180],
[0.800000, 0.800000, 0.810080, 0.950000, 0.961550],
[1.038590, 1.063200, 1.034440, 1.112780, 1.037800]]
coeffs[5, 2, :, :] = [
[1.017610, 1.028360, 1.058960, 1.133180, 1.045620],
[0.920000, 0.998970, 1.033590, 1.089030, 1.022060],
[0.912370, 0.949930, 0.979770, 1.020420, 0.981770],
[0.847160, 0.935300, 0.930540, 0.955050, 0.946560],
[0.880260, 0.867110, 0.874130, 0.972650, 0.883420],
[0.627150, 0.627150, 0.700000, 0.774070, 0.845130],
[0.973700, 1.006240, 1.026190, 1.071960, 1.017240]]
coeffs[5, 3, :, :] = [
[1.028710, 1.017570, 1.025900, 1.081790, 1.024240],
[0.924980, 0.985500, 1.014100, 1.092210, 0.999610],
[0.828570, 0.934920, 0.994950, 1.024590, 0.949710],
[0.900810, 0.901330, 0.928830, 0.979570, 0.913100],
[0.761030, 0.845150, 0.805360, 0.936790, 0.853460],
[0.626400, 0.546750, 0.730500, 0.850000, 0.689050],
[0.957630, 0.985480, 0.991790, 1.050220, 0.987900]]
coeffs[5, 4, :, :] = [
[0.992730, 0.993880, 1.017150, 1.059120, 1.017450],
[0.975610, 0.987160, 1.026820, 1.075440, 1.007250],
[0.871090, 0.933190, 0.974690, 0.979840, 0.952730],
[0.828750, 0.868090, 0.834920, 0.905510, 0.871530],
[0.781540, 0.782470, 0.767910, 0.764140, 0.795890],
[0.743460, 0.693390, 0.514870, 0.630150, 0.715660],
[0.934760, 0.957870, 0.959640, 0.972510, 0.981640]]
coeffs[5, 5, :, :] = [
[0.965840, 0.941240, 0.987100, 1.022540, 1.011160],
[0.988630, 0.994770, 0.976590, 0.950000, 1.034840],
[0.958200, 1.018080, 0.974480, 0.920000, 0.989870],
[0.811720, 0.869090, 0.812020, 0.850000, 0.821050],
[0.682030, 0.679480, 0.632450, 0.746580, 0.738550],
[0.668290, 0.445860, 0.500000, 0.678920, 0.696510],
[0.926940, 0.953350, 0.959050, 0.876210, 0.991490]]
coeffs[5, 6, :, :] = [
[0.948940, 0.997760, 0.850000, 0.826520, 0.998470],
[1.017860, 0.970000, 0.850000, 0.700000, 0.988560],
[1.000000, 0.950000, 0.850000, 0.606240, 0.947260],
[1.000000, 0.746140, 0.751740, 0.598390, 0.725230],
[0.922210, 0.500000, 0.376800, 0.517110, 0.548630],
[0.500000, 0.450000, 0.429970, 0.404490, 0.539940],
[0.960430, 0.881630, 0.775640, 0.596350, 0.937680]]
coeffs[6, 1, :, :] = [
[1.030000, 1.040000, 1.000000, 1.000000, 1.049510],
[1.050000, 0.990000, 0.990000, 0.950000, 0.996530],
[1.050000, 0.990000, 0.990000, 0.820000, 0.971940],
[1.050000, 0.790000, 0.880000, 0.820000, 0.951840],
[1.000000, 0.530000, 0.440000, 0.710000, 0.928730],
[0.540000, 0.470000, 0.500000, 0.550000, 0.773950],
[1.038270, 0.920180, 0.910930, 0.821140, 1.034560]]
coeffs[6, 2, :, :] = [
[1.041020, 0.997520, 0.961600, 1.000000, 1.035780],
[0.948030, 0.980000, 0.900000, 0.950360, 0.977460],
[0.950000, 0.977250, 0.869270, 0.800000, 0.951680],
[0.951870, 0.850000, 0.748770, 0.700000, 0.883850],
[0.900000, 0.823190, 0.727450, 0.600000, 0.839870],
[0.850000, 0.805020, 0.692310, 0.500000, 0.788410],
[1.010090, 0.895270, 0.773030, 0.816280, 1.011680]]
coeffs[6, 3, :, :] = [
[1.022450, 1.004600, 0.983650, 1.000000, 1.032940],
[0.943960, 0.999240, 0.983920, 0.905990, 0.978150],
[0.936240, 0.946480, 0.850000, 0.850000, 0.930320],
[0.816420, 0.885000, 0.644950, 0.817650, 0.865310],
[0.742960, 0.765690, 0.561520, 0.700000, 0.827140],
[0.643870, 0.596710, 0.474460, 0.600000, 0.651200],
[0.971740, 0.940560, 0.714880, 0.864380, 1.001650]]
coeffs[6, 4, :, :] = [
[0.995260, 0.977010, 1.000000, 1.000000, 1.035250],
[0.939810, 0.975250, 0.939980, 0.950000, 0.982550],
[0.876870, 0.879440, 0.850000, 0.900000, 0.917810],
[0.873480, 0.873450, 0.751470, 0.850000, 0.863040],
[0.761470, 0.702360, 0.638770, 0.750000, 0.783120],
[0.734080, 0.650000, 0.600000, 0.650000, 0.715660],
[0.942160, 0.919100, 0.770340, 0.731170, 0.995180]]
coeffs[6, 5, :, :] = [
[0.952560, 0.916780, 0.920000, 0.900000, 1.005880],
[0.928620, 0.994420, 0.900000, 0.900000, 0.983720],
[0.913070, 0.850000, 0.850000, 0.800000, 0.924280],
[0.868090, 0.807170, 0.823550, 0.600000, 0.844520],
[0.769570, 0.719870, 0.650000, 0.550000, 0.733500],
[0.580250, 0.650000, 0.600000, 0.500000, 0.628850],
[0.904770, 0.852650, 0.708370, 0.493730, 0.949030]]
coeffs[6, 6, :, :] = [
[0.911970, 0.800000, 0.800000, 0.800000, 0.956320],
[0.912620, 0.682610, 0.750000, 0.700000, 0.950110],
[0.653450, 0.659330, 0.700000, 0.600000, 0.856110],
[0.648440, 0.600000, 0.641120, 0.500000, 0.695780],
[0.570000, 0.550000, 0.598800, 0.400000, 0.560150],
[0.475230, 0.500000, 0.518640, 0.339970, 0.520230],
[0.743440, 0.592190, 0.603060, 0.316930, 0.794390]]
return coeffs[1:, 1:, :, :]
[docs]def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1.1,
zenith_threshold_for_zero_dni=88.0,
zenith_threshold_for_clearsky_limit=80.0):
"""
Determine DNI from GHI and DHI.
When calculating the DNI from GHI and DHI the calculated DNI may be
unreasonably high or negative for zenith angles close to 90 degrees
(sunrise/sunset transitions). This function identifies unreasonable DNI
values and sets them to NaN. If the clearsky DNI is given unreasonably high
values are cut off.
Parameters
----------
ghi : Series
Global horizontal irradiance.
dhi : Series
Diffuse horizontal irradiance.
zenith : Series
True (not refraction-corrected) zenith angles in decimal
degrees. Angles must be >=0 and <=180.
clearsky_dni : Series, optional
Clearsky direct normal irradiance.
clearsky_tolerance : float, default 1.1
If 'clearsky_dni' is given this parameter can be used to allow a
tolerance by how much the calculated DNI value can be greater than
the clearsky value before it is identified as an unreasonable value.
zenith_threshold_for_zero_dni : float, default 88.0
Non-zero DNI values for zenith angles greater than or equal to
'zenith_threshold_for_zero_dni' will be set to NaN.
zenith_threshold_for_clearsky_limit : float, default 80.0
DNI values for zenith angles greater than or equal to
'zenith_threshold_for_clearsky_limit' and smaller the
'zenith_threshold_for_zero_dni' that are greater than the clearsky DNI
(times allowed tolerance) will be corrected. Only applies if
'clearsky_dni' is not None.
Returns
-------
dni : Series
The modeled direct normal irradiance.
"""
# calculate DNI
dni = (ghi - dhi) / tools.cosd(zenith)
# cutoff negative values
dni[dni < 0] = float('nan')
# set non-zero DNI values for zenith angles >=
# zenith_threshold_for_zero_dni to NaN
dni[(zenith >= zenith_threshold_for_zero_dni) & (dni != 0)] = float('nan')
# correct DNI values for zenith angles greater or equal to the
# zenith_threshold_for_clearsky_limit and smaller than the
# upper_cutoff_zenith that are greater than the clearsky DNI (times
# clearsky_tolerance)
if clearsky_dni is not None:
max_dni = clearsky_dni * clearsky_tolerance
dni[(zenith >= zenith_threshold_for_clearsky_limit) &
(zenith < zenith_threshold_for_zero_dni) &
(dni > max_dni)] = max_dni
return dni
[docs]def complete_irradiance(solar_zenith,
ghi=None,
dhi=None,
dni=None,
dni_clear=None):
r"""
Use the component sum equations to calculate the missing series, using
the other available time series. One of the three parameters (ghi, dhi,
dni) is passed as None, and the other associated series passed are used to
calculate the missing series value.
The "component sum" or "closure" equation relates the three
primary irradiance components as follows:
.. math::
GHI = DHI + DNI \cos(\theta_z)
Parameters
----------
solar_zenith : Series
Zenith angles in decimal degrees, with datetime index.
Angles must be >=0 and <=180. Must have the same datetime index
as ghi, dhi, and dni series, when available.
ghi : Series, optional
Pandas series of dni data, with datetime index. Must have the same
datetime index as dni, dhi, and zenith series, when available.
dhi : Series, optional
Pandas series of dni data, with datetime index. Must have the same
datetime index as ghi, dni, and zenith series, when available.
dni : Series, optional
Pandas series of dni data, with datetime index. Must have the same
datetime index as ghi, dhi, and zenith series, when available.
dni_clear : Series, optional
Pandas series of clearsky dni data. Must have the same datetime index
as ghi, dhi, dni, and zenith series, when available. See
:py:func:`dni` for details.
Returns
-------
component_sum_df : Dataframe
Pandas series of 'ghi', 'dhi', and 'dni' columns with datetime index
"""
if ghi is not None and dhi is not None and dni is None:
dni = pvlib.irradiance.dni(ghi, dhi, solar_zenith,
clearsky_dni=dni_clear,
clearsky_tolerance=1.1)
elif dni is not None and dhi is not None and ghi is None:
ghi = (dhi + dni * tools.cosd(solar_zenith))
elif dni is not None and ghi is not None and dhi is None:
dhi = (ghi - dni * tools.cosd(solar_zenith))
else:
raise ValueError(
"Please check that exactly one of ghi, dhi and dni parameters "
"is set to None"
)
# Merge the outputs into a master dataframe containing 'ghi', 'dhi',
# and 'dni' columns
component_sum_df = pd.DataFrame({'ghi': ghi,
'dhi': dhi,
'dni': dni})
return component_sum_df
[docs]def louche(ghi, solar_zenith, datetime_or_doy, max_zenith=90):
"""
Determine DNI and DHI from GHI using the Louche model.
Parameters
----------
ghi : numeric
Global horizontal irradiance. [W/m^2]
solar_zenith : numeric
True (not refraction-corrected) zenith angles in decimal
degrees. Angles must be >=0 and <=90.
datetime_or_doy : numeric, pandas.DatetimeIndex
Day of year or array of days of year e.g.
pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
Returns
-------
data: OrderedDict or DataFrame
Contains the following keys/columns:
* ``dni``: the modeled direct normal irradiance in W/m^2.
* ``dhi``: the modeled diffuse horizontal irradiance in
W/m^2.
* ``kt``: Ratio of global to extraterrestrial irradiance
on a horizontal plane.
References
-------
.. [1] Louche A, Notton G, Poggi P, Simonnot G. Correlations for direct
normal and global horizontal irradiation on a French Mediterranean site.
Solar Energy 1991;46:261-6. :doi:`10.1016/0038-092X(91)90072-5`
"""
I0 = get_extra_radiation(datetime_or_doy)
Kt = clearness_index(ghi, solar_zenith, I0)
kb = -10.627*Kt**5 + 15.307*Kt**4 - 5.205 * \
Kt**3 + 0.994*Kt**2 - 0.059*Kt + 0.002
dni = kb*I0
dhi = ghi - dni*tools.cosd(solar_zenith)
bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0)
dni = np.where(bad_values, 0, dni)
# ensure that closure relationship remains valid
dhi = np.where(bad_values, ghi, dhi)
data = OrderedDict()
data['dni'] = dni
data['dhi'] = dhi
data['kt'] = Kt
if isinstance(datetime_or_doy, pd.DatetimeIndex):
data = pd.DataFrame(data, index=datetime_or_doy)
return data