Add optional return_components parameter to irradiance.haydavies

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

@@ -10,6 +10,7 @@ Deprecations
 Enhancements
 ~~~~~~~~~~~~
 * Multiple code style issues fixed that were reported by LGTM analysis. (:issue:`1275`, :pull:`1559`)
+* Option to return individual components of Hay Davies sky transposition model (:issue:`1553`, :pull:`1568`)
 
 Bug fixes
 ~~~~~~~~~
@@ -34,4 +35,5 @@ Requirements
 
 Contributors
 ~~~~~~~~~~~~
-* Christian Orner (:ghuser:`chrisorner`)
+* Christian Orner (:ghuser:`chrisorner`)
+* Saurabh Aneja (:ghuser:`spaneja`)

pvlib/irradiance.py

@@ -739,7 +739,8 @@ def klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith,
 
 
 def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
-              solar_zenith=None, solar_azimuth=None, projection_ratio=None):
+              solar_zenith=None, solar_azimuth=None, projection_ratio=None,
+              return_components=False):
     r'''
     Determine diffuse irradiance from the sky on a tilted surface using
     Hay & Davies' 1980 model
@@ -790,6 +791,10 @@ def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
         projection. Must supply ``solar_zenith`` and ``solar_azimuth``
         or supply ``projection_ratio``.
 
+    return_components: bool (optional, default=False)
+        Flag used to decide whether to return the calculated diffuse components
+        or not.
+
     Returns
     --------
     sky_diffuse : numeric
@@ -833,7 +838,33 @@ def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     sky_diffuse = dhi * (AI * Rb + term1 * term2)
     sky_diffuse = np.maximum(sky_diffuse, 0)
 
-    return sky_diffuse
+    # we've preserved the input type until now, so don't ruin it!
+    # if isinstance(sky_diffuse, pd.Series):
+    #     sky_diffuse[np.isnan(surface_tilt)] = 0
+    # else:
+    #     sky_diffuse = np.where(np.isnan(surface_tilt), 0, sky_diffuse)
+
+    if return_components:
+        diffuse_components = OrderedDict()
+        diffuse_components['sky_diffuse'] = sky_diffuse
+
+        # Calculate the individual components
+        diffuse_components['isotropic'] = np.maximum(dhi * term1 * term2, 0)
+        diffuse_components['circumsolar'] = np.maximum(dhi * (AI * Rb), 0)
+        diffuse_components['horizon'] = np.where(
+            np.isnan(diffuse_components['isotropic']), np.nan, 0.)
+
+        # Set values of components to 0 when sky_diffuse is 0
+        mask = sky_diffuse == 0
+        if isinstance(sky_diffuse, pd.Series):
+            diffuse_components = pd.DataFrame(diffuse_components)
+            diffuse_components.loc[mask] = 0
+        else:
+            diffuse_components = {k: np.where(mask, 0, v) for k, v in
+                                  diffuse_components.items()}
+        return diffuse_components
+    else:
+        return sky_diffuse
 
 
 def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,

pvlib/tests/test_irradiance.py

@@ -207,6 +207,36 @@ def test_haydavies(irrad_data, ephem_data, dni_et):
     assert_allclose(result, [0, 27.1775, 102.9949, 33.1909], atol=1e-4)
 
 
+def test_haydavies_components(irrad_data, ephem_data, dni_et):
+    result = irradiance.haydavies(
+        40, 180, irrad_data['dhi'], irrad_data['dni'], dni_et,
+        ephem_data['apparent_zenith'], ephem_data['azimuth'],
+        return_components=True)
+    result_val = irradiance.haydavies(
+        40, 180, irrad_data['dhi'].values, irrad_data['dni'].values, dni_et,
+        ephem_data['apparent_zenith'].values, ephem_data['azimuth'].values,
+        return_components=True)
+    expected = pd.DataFrame(np.array(
+        [[0, 27.1775, 102.9949, 33.1909],
+         [0, 27.1775, 30.1818, 27.9837],
+         [0, 0, 72.8130, 5.2071],
+         [0, 0, 0, 0]]).T,
+        columns=['sky_diffuse', 'isotropic', 'circumsolar', 'horizon'],
+        index=irrad_data.index
+    )
+    # values
+    assert_allclose(result.sky_diffuse, [0, 27.1775, 102.9949, 33.1909],
+                    atol=1e-4)
+    assert_allclose(result.isotropic, [0, 27.1775, 30.1818, 27.9837],
+                    atol=1e-4)
+    assert_allclose(result.circumsolar, [0, 0, 72.8130, 5.2071], atol=1e-4)
+    assert_allclose(result.horizon, [0, 0, 0, 0], atol=1e-4)
+    assert_frame_equal(result, expected, check_less_precise=4)
+    assert isinstance(result, pd.DataFrame)
+    assert_allclose(result_val['sky_diffuse'], [0, 27.1775, 102.9949, 33.1909],
+                    atol=1e-4)
+
+
 def test_reindl(irrad_data, ephem_data, dni_et):
     result = irradiance.reindl(
         40, 180, irrad_data['dhi'], irrad_data['dni'], irrad_data['ghi'],