Add irradiance.clearsky_index (#571)

* basic code for irradiance.clearsky_index, still needs full docstring and tests

* Add a float conversion to irradiance.clearsky_index to allow integer inputs

* added tests for irradiance.clearsky_index

* fix typo in irradiance.clearsky_index test

* fix typo in irradiance.clearsky_index test

* Update pvlib/test/test_irradiance.py to remove extra spaces

* fixed incorrect documentation for irradiance.clearsky_index

* rename poorly-named variables in irradiance.clearsky_index

* truncate nan/inf values to zero in irradiance.clearsky_index

* update whatsnew and api.rst

* improved formatting to satisfy stickler

* preserve input nans for irradiance.clearsky_index

* remove trailing whitespace

* refactor irradiance.clearsky_index to simplify type-handling logic.  also changed behavior for input inf values to set clearsky index to zero

* revert change to inf/nan handling; input infs yield clearsky_index=nan

Author: kevinsa5

docs/sphinx/source/api.rst

@@ -190,6 +190,7 @@ Clearness index models
 
    irradiance.clearness_index
    irradiance.clearness_index_zenith_independent
+   irradiance.clearsky_index
 
 
 PV Modeling

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

@@ -63,6 +63,8 @@ Enhancements
 * Created :py:func:`pvlib.pvsystem.pvsyst_celltemp` to implement PVsyst's cell temperature model. (:issue:`552`)
 * Created :py:func:`pvlib.bifacial.pvfactors_timeseries` to use open-source `pvfactors` package to calculate back surface irradiance (:issue:`421`)
 * Add `PVSystem` class method :py:func:`~pvlib.pvsystem.PVSystem.pvsyst_celltemp` (:issue:`633`)
+* Add :py:func:`pvlib.irradiance.clearsky_index` to calculate clear-sky index
+  from measured GHI and modeled clear-sky GHI. (:issue:`551`)
 
 
 Bug fixes
@@ -102,3 +104,4 @@ Contributors
 * Jonathan Gaffiot (:ghuser:`jgaffiot`)
 * Marc Anoma (:ghuser:`anomam`)
 * Anton Driesse (:ghuser:`adriesse`)
+* Kevin Anderson (:ghuser:`kevinsa5`)

pvlib/irradiance.py

@@ -1198,6 +1198,48 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
         return sky_diffuse
 
 
+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
+
+
 def clearness_index(ghi, solar_zenith, extra_radiation, min_cos_zenith=0.065,
                     max_clearness_index=2.0):
     """

pvlib/test/test_irradiance.py

@@ -766,6 +766,45 @@ def test_kt_kt_prime_factor(airmass_kt):
     assert_allclose(out, expected, atol=1e-5)
 
 
+def test_clearsky_index():
+    ghi = np.array([-1., 0., 1., 500., 1000., np.nan])
+    ghi_measured, ghi_modeled = np.meshgrid(ghi, ghi)
+    # default max_clearsky_index
+    with np.errstate(invalid='ignore', divide='ignore'):
+        out = irradiance.clearsky_index(ghi_measured, ghi_modeled)
+    expected = np.array(
+        [[1.    , 0.    , 0.    , 0.    , 0.    , np.nan],
+         [0.    , 0.    , 0.    , 0.    , 0.    , np.nan],
+         [0.    , 0.    , 1.    , 2.    , 2.    , np.nan],
+         [0.    , 0.    , 0.002 , 1.    , 2.    , np.nan],
+         [0.    , 0.    , 0.001 , 0.5   , 1.    , np.nan],
+         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]])
+    assert_allclose(out, expected, atol=0.001)
+    # specify max_clearsky_index
+    with np.errstate(invalid='ignore', divide='ignore'):
+        out = irradiance.clearsky_index(ghi_measured, ghi_modeled,
+                                        max_clearsky_index=1.5)
+    expected = np.array(
+        [[1.    , 0.    , 0.    , 0.    , 0.    , np.nan],
+         [0.    , 0.    , 0.    , 0.    , 0.    , np.nan],
+         [0.    , 0.    , 1.    , 1.5   , 1.5   , np.nan],
+         [0.    , 0.    , 0.002 , 1.    , 1.5   , np.nan],
+         [0.    , 0.    , 0.001 , 0.5   , 1.    , np.nan],
+         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]])
+    assert_allclose(out, expected, atol=0.001)
+    # scalars
+    out = irradiance.clearsky_index(10, 1000)
+    expected = 0.01
+    assert_allclose(out, expected, atol=0.001)
+    # series
+    times = pd.DatetimeIndex(start='20180601', periods=2, freq='12H')
+    ghi_measured = pd.Series([100,  500], index=times)
+    ghi_modeled = pd.Series([500, 1000], index=times)
+    out = irradiance.clearsky_index(ghi_measured, ghi_modeled)
+    expected = pd.Series([0.2, 0.5], index=times)
+    assert_series_equal(out, expected)
+
+
 def test_clearness_index():
     ghi = np.array([-1, 0, 1, 1000])
     solar_zenith = np.array([180, 90, 89.999, 0])