Repository: pvlib/pvanalytics (python)

Issue:
Values outside QC test domain are incorrectly flagged as in error
I've found that the QC checks in pvanalytics flag values outside the testing domain as erroneously, which I strongly believe is against the purpose of the testing domain.

As an example, the [BSRN comparison and diffuse ratio tests](https://epic.awi.de/30083/1/BSRN_recommended_QC_tests_V2.pdf) are stated to only be performed for GHI>50 W/m^2. The reason for this threshold is that the uncertainty at lower irradiance is too high to perform the test reliably. **This is NOT saying that values lower than 50 W/m^2 should be flagged.**

If values lower than 50 W/m^2 should be flagged, it certainly does not make sense for the BSRN lower limits on the individual components to be -4 W/m^2. For example, the below data point which matches the closure equation perfectly (i.e., GHI = DHI+DNI*cos(Z)) is flagged as erroneously because it's below 50 W/m^2.

        qcrad_consistency_flag = check_irradiance_consistency_qcrad(
            solar_zenith=pd.Series(80),
            ghi=40,
            dhi=40,
            dni=0)

I have a series of QC checks that I would like to add, but this issue is a blocker for now.

This issue is related to https://github.com/pvlib/pvanalytics/issues/188#issuecomment-1481728592 and https://github.com/pvlib/pvanalytics/pull/167#issuecomment-1373484368

Pull request:
Add outside_domain parameter to irradiance ratio QC tests
<!-- Thank you for your contribution! Please don't hesitate to ask for help if you are unsure of how to accomplish any of the items. 
You are free to strikethrough any checklist items that do not apply or to add additional items that are not on this list. -->

- [x] Closes #191
- [x] Added tests to cover all new or modified code.
- [x] Clearly documented all new API functions with [PEP257](https://www.python.org/dev/peps/pep-0257/) and [numpydoc](https://numpydoc.readthedocs.io/en/latest/format.html) compliant docstrings.
- ~[ ] Added new API functions to `docs/api.rst`.~
- [x] Non-API functions clearly documented with docstrings or comments as necessary.
- [x] Adds description and name entries in the appropriate "what's new" file 
      in [`docs/whatsnew`](https://github.com/pvlib/pvanalytics/tree/main/docs/whatsnew) 
      for all changes. Includes link to the GitHub Issue with `` :issue:`num` `` 
      or this Pull Request with `` :pull:`num` ``. Includes contributor name 
      and/or GitHub username (link with `` :ghuser:`user` ``).
- [x] Pull request is nearly complete and ready for detailed review.
- [x] Maintainer: Appropriate GitHub Labels and Milestone are assigned to the Pull Request and linked Issue.

<!-- Please provide a brief description of the problem and the proposed solution or new feature (if not already fully described in a linked issue): -->
This PR adds an optional `out_of_bounds` parameter to the irradiance QC tests which allows the user to specify whether to Fail (False) or pass (True) a test when the conditions are outside of the valid test domain.



Patch:
diff --git a/docs/whatsnew/0.2.2.rst b/docs/whatsnew/0.2.2.rst
index 82983d6..71eb44c 100644
--- a/docs/whatsnew/0.2.2.rst
+++ b/docs/whatsnew/0.2.2.rst
@@ -6,6 +6,9 @@
 
 Enhancements
 ~~~~~~~~~~~~
+* Added optional keyword `outside_domain` to
+  :py:func:`~pvanalytics.quality.irradiance.check_irradiance_consistency_qcrad`.
+  (:pull:`214`)
 
 
 Bug Fixes
@@ -26,3 +29,4 @@ Testing
 
 Contributors
 ~~~~~~~~~~~~
+* Adam R. Jensen (:ghuser:`AdamRJensen`)
diff --git a/pvanalytics/quality/irradiance.py b/pvanalytics/quality/irradiance.py
index c4b08fd..6967805 100644
--- a/pvanalytics/quality/irradiance.py
+++ b/pvanalytics/quality/irradiance.py
@@ -269,29 +269,32 @@ def _get_bounds(bounds):
 
 
 def _check_irrad_ratio(ratio, ghi, sza, bounds):
-    # unpack bounds dict
+    # unpack bounds
     ghi_lb, ghi_ub, sza_lb, sza_ub, ratio_lb, ratio_ub = _get_bounds(bounds)
-    # for zenith set inclusive_lower to handle edge cases, e.g., zenith=0
-    return (
+
+    within_domain = (
         quality.util.check_limits(
             sza, lower_bound=sza_lb, upper_bound=sza_ub, inclusive_lower=True)
         & quality.util.check_limits(
-            ghi, lower_bound=ghi_lb, upper_bound=ghi_ub)
-        & quality.util.check_limits(
-            ratio, lower_bound=ratio_lb, upper_bound=ratio_ub)
+            ghi, lower_bound=ghi_lb, upper_bound=ghi_ub, inclusive_lower=True)
     )
 
+    flag = within_domain & quality.util.check_limits(
+        ratio, lower_bound=ratio_lb, upper_bound=ratio_ub)
+
+    return flag, within_domain
+
 
 def check_irradiance_consistency_qcrad(solar_zenith, ghi, dhi, dni,
-                                       param=None):
-    """Check consistency of GHI, DHI and DNI using QCRad criteria.
+                                       param=None, outside_domain=False):
+    r"""Check consistency of GHI, DHI and DNI using QCRad criteria.
 
     Uses criteria given in [1]_ to validate the ratio of irradiance
     components.
 
-    .. warning:: Not valid for night time. While you can pass data
-       from night time to this function, be aware that the truth
-       values returned for that data will not be valid.
+    .. warning:: Not valid for night time or low irradiance. When the input
+       data fall outside the test domain, the returned value is set by the
+       ``outside_domain`` parameter.
 
     Parameters
     ----------
@@ -303,12 +306,15 @@ def check_irradiance_consistency_qcrad(solar_zenith, ghi, dhi, dni,
         Diffuse horizontal irradiance in :math:`W/m^2`
     dni : Series
         Direct normal irradiance in :math:`W/m^2`
-    param : dict
+    param : dict, optional
         keys are 'ghi_ratio' and 'dhi_ratio'. For each key, value is a dict
         with keys 'high_zenith' and 'low_zenith'; for each of these keys,
         value is a dict with keys 'zenith_bounds', 'ghi_bounds', and
         'ratio_bounds' and value is an ordered pair [lower, upper]
         of float.
+    outside_domain : scalar, default False
+        Value to return when the tests are not applicable, i.e., when the
+        input data fall outside the test domain.
 
     Returns
     -------
@@ -319,6 +325,15 @@ def check_irradiance_consistency_qcrad(solar_zenith, ghi, dhi, dni,
 
     Notes
     -----
+    The QCRad algorithm checks that the input GHI is consistent with the
+    component sum :math:`DNI \times \cos ( zenith ) + DHI` of input DNI and
+    DHI, and that the ratio :math:`\frac{DHI}{GHI}` is reasonable.
+
+    In these two parts, the ``ghi_bounds`` are applied differently. In the
+    components test, the bounds are applied to the component sum of diffuse and
+    direct irradiance, whereas in the diffuse ratio test the bounds are applied
+    to the measured ``ghi``.
+
     Copyright (c) 2019 SolarArbiter. See the file
     LICENSES/SOLARFORECASTARBITER_LICENSE at the top level directory
     of this distribution and at `<https://github.com/pvlib/
@@ -341,18 +356,29 @@ def check_irradiance_consistency_qcrad(solar_zenith, ghi, dhi, dni,
     dhi_ratio = dhi / ghi
 
     bounds = param['ghi_ratio']
-    consistent_components = (
-        _check_irrad_ratio(ratio=ghi_ratio, ghi=component_sum,
-                           sza=solar_zenith, bounds=bounds['high_zenith'])
-        | _check_irrad_ratio(ratio=ghi_ratio, ghi=component_sum,
-                             sza=solar_zenith, bounds=bounds['low_zenith']))
+    flag_lz, within_domain_lz = _check_irrad_ratio(
+        ratio=ghi_ratio, ghi=component_sum, sza=solar_zenith,
+        bounds=bounds['low_zenith'])
+    flag_hz, within_domain_hz = _check_irrad_ratio(
+        ratio=ghi_ratio, ghi=component_sum, sza=solar_zenith,
+        bounds=bounds['high_zenith'])
+
+    consistent_components = ((flag_lz & within_domain_lz) |
+                             (flag_hz & within_domain_hz))
+    consistent_components[~(within_domain_lz | within_domain_hz)] = \
+        outside_domain
 
     bounds = param['dhi_ratio']
-    diffuse_ratio_limit = (
-        _check_irrad_ratio(ratio=dhi_ratio, ghi=ghi, sza=solar_zenith,
-                           bounds=bounds['high_zenith'])
-        | _check_irrad_ratio(ratio=dhi_ratio, ghi=ghi, sza=solar_zenith,
-                             bounds=bounds['low_zenith']))
+    flag_lz, within_domain_lz = _check_irrad_ratio(
+        ratio=dhi_ratio, ghi=ghi, sza=solar_zenith,
+        bounds=bounds['low_zenith'])
+    flag_hz, within_domain_hz = _check_irrad_ratio(
+        ratio=dhi_ratio, ghi=ghi, sza=solar_zenith,
+        bounds=bounds['high_zenith'])
+    diffuse_ratio_limit = ((flag_lz & within_domain_lz) |
+                           (flag_hz & within_domain_hz))
+    diffuse_ratio_limit[~(within_domain_lz | within_domain_hz)] = \
+        outside_domain
 
     return consistent_components, diffuse_ratio_limit
 


Test patch:
diff --git a/pvanalytics/tests/quality/test_irradiance.py b/pvanalytics/tests/quality/test_irradiance.py
index fc7d708..2b28e1f 100644
--- a/pvanalytics/tests/quality/test_irradiance.py
+++ b/pvanalytics/tests/quality/test_irradiance.py
@@ -42,23 +42,28 @@ def irradiance_qcrad():
     output = pd.DataFrame(
         columns=['ghi', 'dhi', 'dni', 'solar_zenith', 'dni_extra',
                  'ghi_limit_flag', 'dhi_limit_flag', 'dni_limit_flag',
-                 'consistent_components', 'diffuse_ratio_limit'],
-        data=np.array([[-100, 100, 100, 30, 1370, 0, 1, 1, 0, 0],
-                       [100, -100, 100, 30, 1370, 1, 0, 1, 0, 0],
-                       [100, 100, -100, 30, 1370, 1, 1, 0, 0, 1],
-                       [1000, 100, 900, 0, 1370, 1, 1, 1, 1, 1],
-                       [1000, 200, 800, 15, 1370, 1, 1, 1, 1, 1],
-                       [1000, 200, 800, 60, 1370, 0, 1, 1, 0, 1],
-                       [1000, 300, 850, 80, 1370, 0, 0, 1, 0, 1],
-                       [1000, 500, 800, 90, 1370, 0, 0, 1, 0, 1],
-                       [500, 100, 1100, 0, 1370, 1, 1, 1, 0, 1],
-                       [1000, 300, 1200, 0, 1370, 1, 1, 1, 0, 1],
-                       [500, 600, 100, 60, 1370, 1, 1, 1, 0, 0],
-                       [500, 600, 400, 80, 1370, 0, 0, 1, 0, 0],
-                       [500, 500, 300, 80, 1370, 0, 0, 1, 1, 1],
-                       [0, 0, 0, 93, 1370, 1, 1, 1, 0, 0]]))
+                 'consistent_components', 'diffuse_ratio_limit',
+                 'consistent_components_outside_domain',
+                 'diffuse_ratio_limit_outside_domain',
+                 ],
+        data=np.array([[-100, 100, 100, 30, 1370, 0, 1, 1, 0, 0, 0, 1],
+                       [100, -100, 100, 30, 1370, 1, 0, 1, 0, 0, 1, 0],
+                       [100, 100, -100, 30, 1370, 1, 1, 0, 0, 1, 1, 1],
+                       [1000, 100, 900, 0, 1370, 1, 1, 1, 1, 1, 1, 1],
+                       [1000, 200, 800, 15, 1370, 1, 1, 1, 1, 1, 1, 1],
+                       [1000, 200, 800, 60, 1370, 0, 1, 1, 0, 1, 0, 1],
+                       [1000, 300, 850, 80, 1370, 0, 0, 1, 0, 1, 0, 1],
+                       [1000, 500, 800, 90, 1370, 0, 0, 1, 0, 1, 0, 1],
+                       [500, 100, 1100, 0, 1370, 1, 1, 1, 0, 1, 0, 1],
+                       [1000, 300, 1200, 0, 1370, 1, 1, 1, 0, 1, 0, 1],
+                       [500, 600, 100, 60, 1370, 1, 1, 1, 0, 0, 0, 0],
+                       [500, 600, 400, 80, 1370, 0, 0, 1, 0, 0, 0, 0],
+                       [500, 500, 300, 80, 1370, 0, 0, 1, 1, 1, 1, 1],
+                       [0, 0, 0, 93, 1370, 1, 1, 1, 0, 0, 1, 1],
+                       [100, 100, 0, 95, 1370, 0, 0, 1, 0, 0, 1, 1],
+                       ]))
     dtypes = ['float64', 'float64', 'float64', 'float64', 'float64',
-              'bool', 'bool', 'bool', 'bool', 'bool']
+              'bool', 'bool', 'bool', 'bool', 'bool', 'bool', 'bool']
     for (col, typ) in zip(output.columns, dtypes):
         output[col] = output[col].astype(typ)
     return output
@@ -179,6 +184,20 @@ def test_check_irradiance_consistency_qcrad(irradiance_qcrad):
                         check_names=False)
 
 
+def test_check_irradiance_consistency_qcrad_outside_domain(irradiance_qcrad):
+    expected = irradiance_qcrad
+    cons_comp, diffuse = irradiance.check_irradiance_consistency_qcrad(
+        expected['solar_zenith'], expected['ghi'],
+        expected['dhi'], expected['dni'],
+        outside_domain=True)
+    assert_series_equal(cons_comp,
+                        expected['consistent_components_outside_domain'],
+                        check_names=False)
+    assert_series_equal(diffuse,
+                        expected['diffuse_ratio_limit_outside_domain'],
+                        check_names=False)
+
+
 @pytest.fixture
 def times():
     """One hour of times at 10 minute frequency.