Safely calculate dni from dhi and ghi

pvlib/irradiance.py

@@ -2048,3 +2048,70 @@ def _get_dirint_coeffs():
         [0.743440, 0.592190, 0.603060, 0.316930, 0.794390]]
 
     return coeffs[1:, 1:, :, :]
+
+
+def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1,
+        upper_cutoff_zenith=88.0, lower_cutoff_zenith=80.0, **kwargs):
+    """
+    Determine DNI from GHI and DHI.
+
+    When calculating the DNI from GHI and DHI the calculated DNI may be
+    unreasonably high or negative for zenith angles close to 90 degrees
+    (sunrise/sunset transitions). This function identifies unreasonable DNI
+    values and sets them to NaN. If the clearsky DNI is given unreasonably high
+    values are cut off.
+
+    Parameters
+    ----------
+    ghi : Series
+        Global horizontal irradiance.
+
+    dhi : Series
+        Diffuse horizontal irradiance.
+
+    zenith : Series
+        True (not refraction-corrected) zenith angles in decimal
+        degrees. Angles must be >=0 and <=180.
+
+    clearsky_dni : None or Series
+        Clearsky direct normal irradiance. Default: None.
+
+    clearsky_tolerance : float
+        If 'clearsky_dni' is given this parameter can be used to allow a
+        tolerance by how much the calculated DNI value can be greater than
+        the clearsky value before it is identified as an unreasonable value.
+        Default: 1.
+
+    upper_cutoff_zenith : float
+        Zenith angle above which nonzero DNI values are set to NaN.
+        Default: 88.
+
+    lower_cutoff_zenith : float
+        Zenith angle above which DNI values greater the clearsky DNI (plus
+        allowed tolerance) are corrected. Only applies if 'clearsky_dni' is not
+        None. Default: 80.
+
+    Returns
+    -------
+    dni : Series
+        The modeled direct normal irradiance.
+
+    """
+
+    # calculate DNI
+    dni_tmp = (ghi - dhi) / tools.cosd(zenith)
+    dni = dni_tmp.copy()
+
+    # cutoff negative values
+    dni[dni < 0] = float('nan')
+
+    # set non-zero DNI values for zenith angles >= upper_cutoff_zenith to NaN
+    dni[(zenith >= upper_cutoff_zenith) & (dni != 0)] = float('nan')
+
+    # correct DNI values for zenith angles between lower_cutoff_zenith and
+    # upper_cutoff_zenith that are greater than the clearsky tolerance
+    if clearsky_dni is not None:
+        dni[(zenith >= lower_cutoff_zenith) & (zenith <= upper_cutoff_zenith) &
+            (dni > (clearsky_dni * clearsky_tolerance))] = (clearsky_dni *
+                                                            clearsky_tolerance)
+    return dni

pvlib/modelchain.py

@@ -623,7 +623,7 @@ def effective_irradiance_model(self):
             fd*self.total_irrad['poa_diffuse'])
         return self
 
-    def complete_irradiance(self, times=None, weather=None):
+    def complete_irradiance(self, times=None, weather=None, **kwargs):
         """
         Determine the missing irradiation columns. Only two of the
         following data columns (dni, ghi, dhi) are needed to calculate
@@ -676,18 +676,23 @@ def complete_irradiance(self, times=None, weather=None):
                    "Results can be too high or negative.\n" +
                    "Help to improve this function on github:\n" +
                    "https://github.com/pvlib/pvlib-python \n")
-        warnings.warn(wrn_txt, UserWarning)
+
         if {'ghi', 'dhi'} <= icolumns and 'dni' not in icolumns:
             logging.debug('Estimate dni from ghi and dhi')
-            self.weather.loc[:, 'dni'] = (
-                (self.weather.loc[:, 'ghi'] - self.weather.loc[:, 'dhi']) /
-                tools.cosd(self.solar_position.loc[:, 'zenith']))
+            self.weather.loc[:, 'dni'] = pvlib.irradiance.dni(
+                self.weather.loc[:, 'ghi'], self.weather.loc[:, 'dhi'],
+                self.solar_position.zenith,
+                clearsky_dni=kwargs.get('clearsky_dni',
+                                        self.location.get_clearsky(times).dni),
+                clearsky_tolerance=kwargs.get('clearsky_tolerance', 1))
         elif {'dni', 'dhi'} <= icolumns and 'ghi' not in icolumns:
+            warnings.warn(wrn_txt, UserWarning)
             logging.debug('Estimate ghi from dni and dhi')
             self.weather.loc[:, 'ghi'] = (
                 self.weather.dni * tools.cosd(self.solar_position.zenith) +
                 self.weather.dhi)
         elif {'dni', 'ghi'} <= icolumns and 'dhi' not in icolumns:
+            warnings.warn(wrn_txt, UserWarning)
             logging.debug('Estimate dhi from dni and ghi')
             self.weather.loc[:, 'dhi'] = (
                 self.weather.ghi - self.weather.dni *

pvlib/test/test_irradiance.py

@@ -369,6 +369,7 @@ 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',
@@ -403,3 +404,21 @@ def test_dirindex():
     tol_dirint = 0.2
     assert np.allclose(out.values, dirint_close_values, rtol=tol_dirint, atol=0,
                        equal_nan=True)
+
+
+def test_dni():
+    ghi = pd.Series([90, 100, 100, 100, 100])
+    dhi = pd.Series([100, 90, 50, 50, 50])
+    zenith = pd.Series([80, 100, 85, 70, 85])
+    clearsky_dni = pd.Series([50, 50, 200, 50, 300])
+
+    dni = irradiance.dni(ghi, dhi, zenith,
+                         clearsky_dni=clearsky_dni, clearsky_tolerance=2)
+    assert_series_equal(dni,
+                        pd.Series([float('nan'), float('nan'), 400,
+                                   146.190220008, 573.685662283]))
+
+    dni = irradiance.dni(ghi, dhi, zenith)
+    assert_series_equal(dni,
+                        pd.Series([float('nan'), float('nan'), 573.685662283,
+                                   146.190220008, 573.685662283]))

pvlib/test/test_modelchain.py

@@ -479,22 +479,22 @@ def test_complete_irradiance_clean_run(system, location):
 def test_complete_irradiance(system, location):
     """Check calculations"""
     mc = ModelChain(system, location)
-    times = pd.date_range('2010-07-05 9:00:00', periods=2, freq='H')
-    i = pd.DataFrame({'dni': [30.354455, 77.22822],
-                      'dhi': [372.103976116, 497.087579068],
-                      'ghi': [356.543700, 465.44400]}, index=times)
+    times = pd.date_range('2010-07-05 7:00:00-0700', periods=2, freq='H')
+    i = pd.DataFrame({'dni': [49.756966, 62.153947],
+                      'ghi': [372.103976116, 497.087579068],
+                      'dhi': [356.543700, 465.44400]}, index=times)
 
     mc.complete_irradiance(times, weather=i[['ghi', 'dni']])
     assert_series_equal(mc.weather['dhi'],
-                        pd.Series([372.103976116, 497.087579068],
+                        pd.Series([356.543700, 465.44400],
                                   index=times, name='dhi'))
 
     mc.complete_irradiance(times, weather=i[['dhi', 'dni']])
     assert_series_equal(mc.weather['ghi'],
-                        pd.Series([356.543700, 465.44400],
+                        pd.Series([372.103976116, 497.087579068],
                                   index=times, name='ghi'))
 
     mc.complete_irradiance(times, weather=i[['dhi', 'ghi']])
     assert_series_equal(mc.weather['dni'],
-                        pd.Series([30.354455, 77.22822],
+                        pd.Series([49.756966, 62.153947],
                                   index=times, name='dni'))

pvlib/tools.py

@@ -15,12 +15,12 @@ def cosd(angle):
 
     Parameters
     ----------
-    angle : float
+    angle : float or array-like
         Angle in degrees
 
     Returns
     -------
-    result : float
+    result : float or array-like
         Cosine of the angle
     """