Add diffuse components support: get_sky_diffuse, poa_components, and get_total_irradiance

docs/sphinx/source/whatsnew/v0.15.3.rst

@@ -18,6 +18,8 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
+* Add support for diffuse irradiance components to :py:func:`pvlib.irradiance.isotropic`
+  when ``return_components=True``. (:issue:`2750`, :pull:`2787`)
 * Ensure all timezones are available in all OSs. (:issue:`2795`, :pull:`2809`)
 
 

pvlib/irradiance.py

@@ -281,7 +281,8 @@ def get_total_irradiance(surface_tilt, surface_azimuth,
                          dni, ghi, dhi, dni_extra=None, airmass=None,
                          albedo=0.25, surface_type=None,
                          model='isotropic',
-                         model_perez='allsitescomposite1990'):
+                         model_perez='allsitescomposite1990',
+                         diffuse_components=False):
     r"""
     Determine total in-plane irradiance and its beam, sky diffuse and ground
     reflected components, using the specified sky diffuse irradiance model.
@@ -332,10 +333,17 @@ def get_total_irradiance(surface_tilt, surface_azimuth,
         ``'perez-driesse'``.
     model_perez : str, default 'allsitescomposite1990'
         Used only if ``model='perez'``. See :py:func:`~pvlib.irradiance.perez`.
+    diffuse_components : bool, default False
+        If `True`, returns values for the different diffuse irradiance
+        components available from the selected model
+        (e.g., isotropic, circumsolar, horizon brightening).
+        If `False`, only the total diffuse irradiance is returned.
+        This option is not available for the ``'klucher'`` and
+        ``'king'`` models.
 
     Returns
     -------
-    total_irrad : OrderedDict or DataFrame
+    total_irrad : dict or DataFrame
         Contains keys/columns ``'poa_global', 'poa_direct', 'poa_diffuse',
         'poa_sky_diffuse', 'poa_ground_diffuse'``. [Wm⁻²]
 
@@ -353,7 +361,7 @@ def get_total_irradiance(surface_tilt, surface_azimuth,
     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)
+        model_perez=model_perez, return_components=diffuse_components)
 
     poa_ground_diffuse = get_ground_diffuse(surface_tilt, ghi, albedo,
                                             surface_type)
@@ -366,7 +374,8 @@ 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'):
+                    model_perez='allsitescomposite1990',
+                    return_components=False):
     r"""
     Determine in-plane sky diffuse irradiance component
     using the specified sky diffuse irradiance model.
@@ -408,11 +417,22 @@ def get_sky_diffuse(surface_tilt, surface_azimuth,
         ``'perez-driesse'``.
     model_perez : str, default 'allsitescomposite1990'
         Used only if ``model='perez'``. See :py:func:`~pvlib.irradiance.perez`.
+    return_components : bool, default False
+        If `True`, returns values for the different diffuse irradiance
+        components available from the selected model
+        (e.g., isotropic, circumsolar, horizon brightening).
+        If `False`, only the total diffuse irradiance is returned.
+        This option is not available for the ``'klucher'`` and
+        ``'king'`` models.
 
     Returns
     -------
-    poa_sky_diffuse : numeric
-        Sky diffuse irradiance in the plane of array. [Wm⁻²]
+    numeric, Dict, or DataFrame
+        Return type controlled by ``return_components`` argument.
+        If `False`, total sky diffuse irradiance in the plane of array
+        is returned (numeric). [Wm⁻²]
+        If `True`, the different diffuse components are returned
+        (Dict or DataFrame). [Wm⁻²]
 
     Raises
     ------
@@ -438,33 +458,40 @@ def get_sky_diffuse(surface_tilt, surface_azimuth,
 
     model = model.lower()
 
+    if return_components and model in {'klucher', 'king'}:
+        raise ValueError('return_components is not supported for'
+                         f' model {model}')
+
     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)
+        sky = isotropic(surface_tilt, dhi, return_components=return_components)
     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)
+                        solar_zenith, solar_azimuth,
+                        return_components=return_components)
     elif model == 'reindl':
         sky = reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
-                     solar_zenith, solar_azimuth)
+                     solar_zenith, solar_azimuth,
+                     return_components=return_components)
     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)
+                    model=model_perez, return_components=return_components)
     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)
+                            solar_zenith, solar_azimuth, airmass,
+                            return_components=return_components)
     else:
         raise ValueError(f'invalid model selection {model}')
 
@@ -488,7 +515,7 @@ def poa_components(aoi, dni, poa_sky_diffuse, poa_ground_diffuse):
         Direct normal irradiance, as measured from a TMY file or
         calculated with a clearsky model. See :term:`dni`. [Wm⁻²]
 
-    poa_sky_diffuse : numeric
+    poa_sky_diffuse : numeric, Dict or DataFrame
         Diffuse irradiance in the plane of the modules, as
         calculated by a diffuse irradiance translation function. [Wm⁻²]
 
@@ -499,7 +526,7 @@ def poa_components(aoi, dni, poa_sky_diffuse, poa_ground_diffuse):
 
     Returns
     -------
-    irrads : OrderedDict or DataFrame
+    irrads : Dict or DataFrame
         Contains the following keys:
 
         * ``poa_global`` : Total in-plane irradiance. [Wm⁻²]
@@ -509,21 +536,37 @@ def poa_components(aoi, dni, poa_sky_diffuse, poa_ground_diffuse):
         * ``poa_ground_diffuse`` : In-plane diffuse irradiance from ground.
           [Wm⁻²]
 
+        If ``poa_sky_diffuse`` is a Dict or DataFrame, ``irrads`` will
+        contain additional keys for each of the diffuse components returned by
+        the selected diffuse irradiance model.
+
     Notes
     ------
     Negative beam irradiation due to AOI > 90° or AOI < 0° is set to zero.
     '''
 
+    if isinstance(poa_sky_diffuse, dict):
+        sky_components = poa_sky_diffuse.copy()
+        total_poa_sky_diffuse = sky_components.pop('poa_sky_diffuse')
+    elif isinstance(poa_sky_diffuse, pd.DataFrame):
+        sky_components = poa_sky_diffuse.to_dict(orient='series')
+        total_poa_sky_diffuse = sky_components.pop('poa_sky_diffuse')
+    else:
+        sky_components = {}
+        total_poa_sky_diffuse = poa_sky_diffuse
+
     poa_direct = np.maximum(dni * np.cos(np.radians(aoi)), 0)
-    poa_diffuse = poa_sky_diffuse + poa_ground_diffuse
+    poa_diffuse = total_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
+    irrads = {
+        'poa_global': poa_global,
+        'poa_direct': poa_direct,
+        'poa_diffuse': poa_diffuse,
+        'poa_sky_diffuse': total_poa_sky_diffuse,
+        'poa_ground_diffuse': poa_ground_diffuse,
+        **sky_components
+    }
 
     if isinstance(poa_direct, pd.Series):
         irrads = pd.DataFrame(irrads)
@@ -599,7 +642,7 @@ def get_ground_diffuse(surface_tilt, ghi, albedo=.25, surface_type=None):
     return diffuse_irrad
 
 
-def isotropic(surface_tilt, dhi):
+def isotropic(surface_tilt, dhi, return_components=False):
     r'''
     Determine diffuse irradiance from the sky on a tilted surface using
     the isotropic sky model.
@@ -623,11 +666,27 @@ def isotropic(surface_tilt, dhi):
     dhi : numeric
         Diffuse horizontal irradiance, must be >=0. See :term:`dhi`.
 
+    return_components : bool, default `False`
+        If `False`, ``sky_diffuse`` is returned.
+        If `True`, ``diffuse_components`` is returned.
+        For this model, return_components does not add more information,
+        but it is included for consistency with the other sky diffuse models.
+
     Returns
     -------
-    diffuse : numeric
+    numeric, Dict, or DataFrame
+        Return type controlled by ``return_components`` argument.
+        If `False`, ``sky_diffuse`` is returned.
+        If `True`, ``diffuse_components`` is returned.
+
+    sky_diffuse : numeric
         The sky diffuse component of the solar radiation. [Wm⁻²]
 
+    diffuse_components : Dict (array input) or DataFrame (Series input)
+        Keys/columns are:
+            * poa_sky_diffuse: Total sky diffuse
+            * poa_isotropic
+
     References
     ----------
     .. [1] Loutzenhiser P.G. et al. "Empirical validation of models to
@@ -642,7 +701,18 @@ def isotropic(surface_tilt, dhi):
     '''
     sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5
 
-    return sky_diffuse
+    if return_components:
+        diffuse_components = {
+            'poa_sky_diffuse': sky_diffuse,
+            'poa_isotropic': sky_diffuse
+        }
+
+        if isinstance(sky_diffuse, pd.Series):
+            diffuse_components = pd.DataFrame(diffuse_components)
+
+        return diffuse_components
+    else:
+        return sky_diffuse
 
 
 def klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith,

tests/test_irradiance.py

@@ -168,6 +168,32 @@ def test_isotropic_series(irrad_data):
     assert_allclose(result, [0, 35.728402, 104.601328, 54.777191], atol=1e-4)
 
 
+def test_isotropic_components(irrad_data):
+    keys = ['poa_sky_diffuse', 'poa_isotropic']
+    expected = pd.DataFrame(np.array(
+        [[0, 35.728402, 104.601328, 54.777191],
+         [0, 35.728402, 104.601328, 54.777191]]).T,
+        columns=keys,
+        index=irrad_data.index
+    )
+    # pandas
+    result = irradiance.isotropic(
+        40, irrad_data['dhi'], return_components=True)
+    assert_frame_equal(result, expected, check_less_precise=4)
+    # numpy
+    result = irradiance.isotropic(
+        40, irrad_data['dhi'].values, return_components=True)
+    for key in keys:
+        assert_allclose(result[key], expected[key], atol=1e-4)
+    assert isinstance(result, dict)
+    # scalar
+    result = irradiance.isotropic(
+        40, irrad_data['dhi'].values[-1], return_components=True)
+    for key in keys:
+        assert_allclose(result[key], expected[key].iloc[-1], atol=1e-4)
+    assert isinstance(result, dict)
+
+
 def test_klucher_series_float():
     # klucher inputs
     surface_tilt, surface_azimuth = 40.0, 180.0
@@ -508,6 +534,28 @@ def test_get_total_irradiance(irrad_data, ephem_data, dni_et,
                                           'poa_ground_diffuse']
 
 
+def test_get_total_irradiance_diffuse_components(irrad_data, ephem_data,
+                                                 dni_et, relative_airmass):
+    models = ['perez', 'perez-driesse']
+
+    for model in models:
+        total = irradiance.get_total_irradiance(
+            32, 180,
+            ephem_data['apparent_zenith'], ephem_data['azimuth'],
+            dni=irrad_data['dni'], ghi=irrad_data['ghi'],
+            dhi=irrad_data['dhi'],
+            dni_extra=dni_et, airmass=relative_airmass,
+            model=model,
+            surface_type='urban',
+            diffuse_components=True)
+
+        assert total.columns.tolist() == ['poa_global', 'poa_direct',
+                                          'poa_diffuse', 'poa_sky_diffuse',
+                                          'poa_ground_diffuse',
+                                          'poa_isotropic', 'poa_circumsolar',
+                                          'poa_horizon']
+
+
 @pytest.mark.parametrize('model', ['isotropic', 'klucher',
                                    'haydavies', 'reindl', 'king',
                                    'perez', 'perez-driesse'])
@@ -597,6 +645,59 @@ def test_poa_components(irrad_data, ephem_data, dni_et, relative_airmass):
     assert_frame_equal(out, expected)
 
 
+def test_poa_components_diffuse_components_perez(irrad_data, ephem_data,
+                                                 dni_et, relative_airmass):
+    aoi = irradiance.aoi(40, 180, ephem_data['apparent_zenith'],
+                         ephem_data['azimuth'])
+    gr_sand = irradiance.get_ground_diffuse(40, irrad_data['ghi'],
+                                            surface_type='sand')
+    diff_perez = irradiance.perez(
+        40, 180, irrad_data['dhi'], irrad_data['dni'], dni_et,
+        ephem_data['apparent_zenith'], ephem_data['azimuth'], relative_airmass,
+        return_components=True)
+    out = irradiance.poa_components(
+        aoi, irrad_data['dni'], diff_perez, gr_sand)
+    expected = pd.DataFrame(np.array(
+        [[0.,  -0.,   0.,   0.,
+            0., 0.,   0.,   0.],
+         [35.19456561,   0.,  35.19456561,  31.4635077,
+            3.73105791,   26.841386,   0.000000,   4.622122],
+         [956.18253696, 798.31939281, 157.86314414, 109.08433162,
+          48.77881252,   41.621826,   61.619987,   5.842518],
+         [90.99624896,  33.50143401,  57.49481495,  45.45978964,
+          12.03502531,   31.726961,   4.479664,   9.253165]]),
+        columns=['poa_global', 'poa_direct', 'poa_diffuse', 'poa_sky_diffuse',
+                 'poa_ground_diffuse', 'poa_isotropic', 'poa_circumsolar',
+                 'poa_horizon'],
+        index=irrad_data.index)
+    assert_frame_equal(out, expected)
+
+
+def test_poa_components_diffuse_components_isotropic(irrad_data, ephem_data,
+                                                     dni_et, relative_airmass):
+    aoi = irradiance.aoi(40, 180, ephem_data['apparent_zenith'],
+                         ephem_data['azimuth'])
+    gr_sand = irradiance.get_ground_diffuse(40, irrad_data['ghi'],
+                                            surface_type='sand')
+    diff_isotropic = irradiance.isotropic(
+        40, irrad_data['dhi'], return_components=True)
+    out = irradiance.poa_components(
+        aoi, irrad_data['dni'], diff_isotropic, gr_sand)
+    expected = pd.DataFrame(np.array(
+        [[0.,  -0.,   0.,   0.,
+            0., 0.],
+         [39.459460,   0.000000,   39.459460,   35.728402,
+          3.731058,   35.728402],
+         [951.699533,   798.319393,   153.380140,   104.601328,
+          48.778813,   104.601328],
+         [100.313650,   33.501434,   66.812216,   54.777191,
+          12.035025,   54.777191]]),
+        columns=['poa_global', 'poa_direct', 'poa_diffuse', 'poa_sky_diffuse',
+                 'poa_ground_diffuse', 'poa_isotropic'],
+        index=irrad_data.index)
+    assert_frame_equal(out, expected)
+
+
 @pytest.mark.parametrize('pressure,expected', [
     (93193,  [[830.46567,   0.79742,   0.93505],
               [676.18340,   0.63782,   3.02102]]),