Implement pvlib.irradiance.louche decomposition model (#1705)

* added pvl_louche

* stickler

* add tests

* stickler

* stickler

* made suggested changes

* updated what's new, doc API changes & added notes

* correct comments

* Update pvlib/irradiance.py

Co-authored-by: Cliff Hansen <[email protected]>

* complete PR with formatting changes

* stickler

* Update pvlib/irradiance.py

Co-authored-by: Kevin Anderson <[email protected]>

* done necessary changes

* Update docs/sphinx/source/reference/irradiance/decomposition.rst

Co-authored-by: Cliff Hansen <[email protected]>

* removed SSC comment

* Apply suggestions from code review

---------

Co-authored-by: Cliff Hansen <[email protected]>
Co-authored-by: Kevin Anderson <[email protected]>

Author: 3 people

docs/sphinx/source/reference/irradiance/decomposition.rst

@@ -15,3 +15,5 @@ DNI estimation models
    irradiance.boland
    irradiance.campbell_norman
    irradiance.gti_dirint
+   irradiance.louche
+

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

@@ -23,6 +23,7 @@ Deprecations
 
 Enhancements
 ~~~~~~~~~~~~
+* Added a new irradiance decomposition model :py:func:`pvlib.irradiance.louche`. (:pull:`1705`)
 * Add optional encoding parameter to :py:func:`pvlib.iotools.read_tmy3`.
   (:issue:`1732`, :pull:`1737`)
 * Added function to retrieve horizon data from PVGIS 
@@ -79,3 +80,4 @@ Contributors
 * Andy Lam (:ghuser:`@andylam598`)
 * :ghuser:`ooprathamm`
 * Kevin Anderson (:ghuser:`kandersolar`)
+

pvlib/irradiance.py

@@ -3117,3 +3117,64 @@ def complete_irradiance(solar_zenith,
                                      'dhi': dhi,
                                      'dni': dni})
     return component_sum_df
+
+
+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

pvlib/tests/test_irradiance.py

@@ -246,6 +246,7 @@ def test_haydavies_components(irrad_data, ephem_data, dni_et):
     assert_allclose(result['horizon'], expected['horizon'][-1], atol=1e-4)
     assert isinstance(result, dict)
 
+
 def test_reindl(irrad_data, ephem_data, dni_et):
     result = irradiance.reindl(
         40, 180, irrad_data['dhi'], irrad_data['dni'], irrad_data['ghi'],
@@ -903,8 +904,12 @@ def test_dirindex(times):
     assert np.allclose(out, expected_out, rtol=tolerance, atol=0,
                        equal_nan=True)
     tol_dirint = 0.2
-    assert np.allclose(out.values, dirint_close_values, rtol=tol_dirint, atol=0,
-                       equal_nan=True)
+    assert np.allclose(
+        out.values,
+        dirint_close_values,
+        rtol=tol_dirint,
+        atol=0,
+        equal_nan=True)
 
 
 def test_dirindex_min_cos_zenith_max_zenith():
@@ -1203,3 +1208,20 @@ def test_complete_irradiance():
                                        dhi=None,
                                        dni=i.dni,
                                        dni_clear=clearsky.dni)
+
+
+def test_louche():
+
+    index = pd.DatetimeIndex(['20190101']*3 + ['20190620']*1)
+    ghi = pd.Series([0, 50, 1000, 1000], index=index)
+    zenith = pd.Series([91, 85, 10, 10], index=index)
+    expected = pd.DataFrame(np.array(
+        [[0, 0, 0],
+         [130.089669, 38.661938, 0.405724],
+         [828.498650, 184.088106, 0.718133],
+         [887.407348, 126.074364, 0.768214]]),
+        columns=['dni', 'dhi', 'kt'], index=index)
+
+    out = irradiance.louche(ghi, zenith, index)
+
+    assert_frame_equal(out, expected)