diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 0ac8aa4fde..ccfeb80981 100644
--- a/docs/sphinx/source/api.rst
+++ b/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
diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst
index d4f132baf6..251df4ae92 100644
--- a/docs/sphinx/source/whatsnew/v0.6.1.rst
+++ b/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`)
\ No newline at end of file
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 38c9f35781..c79ce618a9 100644
--- a/pvlib/irradiance.py
+++ b/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):
     """
diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py
index 508c3a9dd2..285c5a6bd8 100644
--- a/pvlib/test/test_irradiance.py
+++ b/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])