Implement dirindex decomposition model (#283)

* Starting implementation of DIRINDEX model

* making call to DIRINT model
* make sure in test that data is returned

* Added calculation of DIRINDEX dni values with clearsky model inputs

* providing clear sky ghi and dni as arguments so that the user can use
 different clear sky models
* only returning positive values from the calculation
* an option could be to provide dni_dirint and dni_dirint_clearsky as
 arguments as well. Not sure if that's preferable
* added docstrings

* Adding parenthesis for Python 3

* Added v0.4.3.txt file and references + updated api.rst

* Updated dirindex test to compare obtained to expected values

* tests for consistency of the algorithm
* doesn't exactly test for validity of implementation, need reference
 values?
* adding comparison to dirint values but not sure how relevant it is

* Making test more robust + decorator for min build on travis

* writing explicit values for input arguments, instead of calling
functions
* adding ``@needs_numpy_1_10`` decorator for continuous integration

Author: anomam

docs/sphinx/source/api.rst

@@ -137,6 +137,7 @@ DNI estimation models
 
    irradiance.disc
    irradiance.dirint
+   irradiance.dirindex
    irradiance.erbs
    irradiance.liujordan
 

docs/sphinx/source/whatsnew.rst

@@ -6,6 +6,7 @@ What's New
 
 These are new features and improvements of note in each release.
 
+.. include:: whatsnew/v0.4.3.txt
 .. include:: whatsnew/v0.4.2.txt
 .. include:: whatsnew/v0.4.1.txt
 .. include:: whatsnew/v0.4.0.txt

docs/sphinx/source/whatsnew/v0.4.3.txt

@@ -0,0 +1,15 @@
+.. _whatsnew_0430:
+
+v0.4.3 (December 13, 2016)
+--------------------------
+
+Enhancements
+~~~~~~~~~~~~
+
+* Adding implementation of Perez's DIRINDEX model based on existing DIRINT
+  model implementation. (:issue:`282`)
+
+Code Contributors
+~~~~~~~~~~~~~~~~~
+
+* Marc Anoma

pvlib/irradiance.py

@@ -1223,7 +1223,7 @@ def dirint(ghi, zenith, times, pressure=101325., use_delta_kt_prime=True,
     Notes
     -----
     DIRINT model requires time series data (ie. one of the inputs must
-    be a vector of length > 2.
+    be a vector of length > 2).
 
     References
     ----------
@@ -1319,6 +1319,86 @@ def dirint(ghi, zenith, times, pressure=101325., use_delta_kt_prime=True,
     return dni
 
 
+def dirindex(ghi, ghi_clearsky, dni_clearsky, zenith, times, pressure=101325.,
+             use_delta_kt_prime=True, temp_dew=None):
+    """
+    Determine DNI from GHI using the DIRINDEX model, which is a modification of
+    the DIRINT model with information from a clear sky model.
+
+    DIRINDEX [1] improves upon the DIRINT model by taking into account turbidity
+    when used with the Ineichen clear sky model results.
+
+    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
+        The site pressure in Pascal. Pressure may be measured or an
+        average pressure may be calculated from site altitude.
+
+    use_delta_kt_prime : bool
+        Indicates if the user would like to utilize the time-series
+        nature of the GHI measurements. A value of ``False`` will not
+        use the time-series improvements, any other numeric value will
+        use time-series improvements. It is recommended that time-series
+        data only be used if the time between measured data points is
+        less than 1.5 hours. If none of the input arguments are vectors,
+        then time-series improvements are not used (because it's not a
+        time-series). If True, input data must be Series.
+
+    temp_dew : None, float, or array-like
+        Surface dew point temperatures, in degrees C. Values of temp_dew
+        may be numeric or NaN. Any single time period point with a
+        DewPtTemp=NaN does not have dew point improvements applied. If
+        DewPtTemp is not provided, then dew point improvements are not
+        applied.
+
+    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)
+
+    dni_dirint_clearsky = dirint(ghi_clearsky, zenith, times, pressure=pressure,
+                                 use_delta_kt_prime=use_delta_kt_prime,
+                                 temp_dew=temp_dew)
+
+    dni_dirindex = dni_clearsky * dni_dirint / dni_dirint_clearsky
+
+    dni_dirindex[dni_dirindex < 0] = 0.
+
+    return dni_dirindex
+
+
 def erbs(ghi, zenith, doy):
     r"""
     Estimate DNI and DHI from GHI using the Erbs model.

pvlib/test/test_irradiance.py

@@ -42,6 +42,7 @@
 dt_np64 = np.datetime64(dt_datetime)
 value = 1383.636203
 
+
 @pytest.mark.parametrize('input, expected', [
     (doy, value),
     (np.float64(doy), value),
@@ -194,7 +195,6 @@ def test_perez_arrays():
     assert_allclose(out, expected, atol=1e-2)
 
 
-
 def test_liujordan():
     expected = pd.DataFrame(np.
         array([[863.859736967, 653.123094076, 220.65905025]]),
@@ -286,6 +286,7 @@ def test_dirint():
     dirint_data = irradiance.dirint(clearsky_data['ghi'], ephem_data['zenith'],
                                     ephem_data.index, pressure=pressure)
 
+
 def test_dirint_value():
     times = pd.DatetimeIndex(['2014-06-24T12-0700','2014-06-24T18-0700'])
     ghi = pd.Series([1038.62, 254.53], index=times)
@@ -318,6 +319,7 @@ def test_dirint_tdew():
     assert_almost_equal(dirint_data.values,
                         np.array([892.9,  636.5]), 1)
 
+
 def test_dirint_no_delta_kt():
     times = pd.DatetimeIndex(['2014-06-24T12-0700','2014-06-24T18-0700'])
     ghi = pd.Series([1038.62, 254.53], index=times)
@@ -328,6 +330,7 @@ def test_dirint_no_delta_kt():
     assert_almost_equal(dirint_data.values,
                         np.array([861.9,  670.4]), 1)
 
+
 def test_dirint_coeffs():
     coeffs = irradiance._get_dirint_coeffs()
     assert coeffs[0,0,0,0] == 0.385230
@@ -365,3 +368,38 @@ def test_erbs_all_scalar():
 
     for k, v in out.items():
         assert_allclose(v, expected[k], 5)
+
+@needs_numpy_1_10
+def test_dirindex():
+    clearsky_data = tus.get_clearsky(times, model='ineichen',
+                                     linke_turbidity=3)
+    ghi = pd.Series([0, 0, 1038.62, 254.53], index=times)
+    ghi_clearsky = pd.Series(
+        np.array([0., 79.73860422, 1042.48031487, 257.20751138]),
+        index=times
+    )
+    dni_clearsky = pd.Series(
+        np.array([0., 316.1949056, 939.95469881, 646.22886049]),
+        index=times
+    )
+    zenith = pd.Series(
+        np.array([124.0390863, 82.85457044, 10.56413562, 72.41687122]),
+        index=times
+    )
+    pressure = 93193.
+    tdew = 10.
+    out = irradiance.dirindex(ghi, ghi_clearsky, dni_clearsky,
+                              zenith, times, pressure=pressure,
+                              temp_dew=tdew)
+    dirint_close_values = irradiance.dirint(ghi, zenith, times,
+                                            pressure=pressure,
+                                            use_delta_kt_prime=True,
+                                            temp_dew=tdew).values
+    expected_out = np.array([np.nan, 0., 748.31562753, 630.72592644])
+
+    tolerance = 1e-8
+    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)