Safely calculate dni from dhi and ghi (#250)

* creating basic function to calculate DNI from GHI and DHI

* replace equation with new function for test purposes

* change types to satisfy type checks

* remove warning if new function is used

* adapt input and output of dni method

* adapt docstring to new inputs

* add clearsky method to correct calculated dni

* add cutoff method to correct calculated dni

* add logging warning and debug messages

* increase code quality

* remove non-ascii characters

* Remove user warning

* Replace input variable location by zenith and clearsky_dni

* Correct docstring

* Remove unnecessary calculation

* Correct negative DNI

Negative values of the calculated DNI will be set to zero or if the new
input value 'set_to_nan' is 'True' will be set to 'NaN'.

* Implement clearsky method

* Include tolerance factor for clearsky method and adapt function call in modelchain

* Enable set_to_nan in cutoff method

* Raise value error if chosen method does not exist

* Debugging info how many values needed correction

* Correct parameter types in docstring

* Correct test_complete_irradiance

* Remove debugging info

* Remove method and introduce new parameters for cutoff zenith angles

* Adapt docstring to new parameters

* Remove commented lines

* Adapt complete_irradiance definition of modelchain to new parameters in DNI calculation

* Add test function for DNI calculation

* Merge upstream changes on master

The only conflict was that both me and the upstream master added a test
at the end of a test file so the solution was easy: just take both.

* Add accidentally un-added files

I had a mishap which removed these files from being version controlled.
This should fix it.

* Remove kwargs from dni function

* Remove unnecessary dni_tmp

* Remove blank line at end of docstring

* Add max_dni to do calculation only once

* Correct comment

* Correct docstring for lower_cutoff_zenith

* Remove kwargs from complete_irradiance

* Correct docstrings

* Use same default value for clearsky_tolerance as in ModelChain

* Change names of lower_ and upper_cutoff_zenith and adapt docstring and comments

lower_cutoff_zenith is renamed to zenith_threshold_for_clearsky_limit
upper_cutoff_zenith is renamed to zenith_threshold_for_zero_dni

* add safely calculate dni from dhi and ghi (#250)

Author: uvchik

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

@@ -13,9 +13,16 @@ Bug fixes
 * Added lower accuracy formulas for equation of time, declination, hour angle
   and solar zenith.
 
+Enhancements
+~~~~~~~~~~~~~
+
+* Added irradiance.dni method that determines DNI from GHI and DHI and corrects
+  unreasonable DNI values during sunrise/sunset transitions
+
 
 Contributors
 ~~~~~~~~~~~~
 
 * Will Holmgren
 * Mark Mikofski
+* Birgit Schachler

pvlib/irradiance.py

@@ -2048,3 +2048,75 @@ 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.1,
+        zenith_threshold_for_zero_dni=88.0,
+        zenith_threshold_for_clearsky_limit=80.0):
+    """
+    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.1.
+
+    zenith_threshold_for_zero_dni : float
+        Non-zero DNI values for zenith angles greater than or equal to
+        'zenith_threshold_for_zero_dni' will be set to NaN. Default: 88.
+
+    zenith_threshold_for_clearsky_limit : float
+        DNI values for zenith angles greater than or equal to
+        'zenith_threshold_for_clearsky_limit' and smaller the
+        'zenith_threshold_for_zero_dni' that are greater than the clearsky DNI
+        (times allowed tolerance) will be corrected. Only applies if
+        'clearsky_dni' is not None. Default: 80.
+
+    Returns
+    -------
+    dni : Series
+        The modeled direct normal irradiance.
+    """
+
+    # calculate DNI
+    dni = (ghi - dhi) / tools.cosd(zenith)
+
+    # cutoff negative values
+    dni[dni < 0] = float('nan')
+
+    # set non-zero DNI values for zenith angles >=
+    # zenith_threshold_for_zero_dni to NaN
+    dni[(zenith >= zenith_threshold_for_zero_dni) & (dni != 0)] = float('nan')
+
+    # correct DNI values for zenith angles greater or equal to the
+    # zenith_threshold_for_clearsky_limit and smaller than the
+    # upper_cutoff_zenith that are greater than the clearsky DNI (times
+    # clearsky_tolerance)
+    if clearsky_dni is not None:
+        max_dni = clearsky_dni * clearsky_tolerance
+        dni[(zenith >= zenith_threshold_for_clearsky_limit) &
+            (zenith < zenith_threshold_for_zero_dni) &
+            (dni > max_dni)] = max_dni
+    return dni

pvlib/modelchain.py

@@ -676,18 +676,22 @@ 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=self.location.get_clearsky(times).dni,
+                clearsky_tolerance=1.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
     """