Update BSRN closure/consistency quality flag function

docs/source/documentation.rst

@@ -18,6 +18,7 @@ Code documentation
    quality.bsrn_limits
    quality.bsrn_limits_flag
    quality.diffuse_fraction_flag
+   quality.closure_flag
    iotools.read_t16
    iotools.write_t16
    processing.resample_to_freq

src/solarpy/quality/__init__.py

@@ -1,3 +1,4 @@
 from solarpy.quality.limits import bsrn_limits  # noqa: F401
 from solarpy.quality.limits import bsrn_limits_flag  # noqa: F401
 from solarpy.quality.comparison import diffuse_fraction_flag  # noqa: F401
+from solarpy.quality.comparison import closure_flag  # noqa: F401

src/solarpy/quality/comparison.py

@@ -1,4 +1,4 @@
-"""Functions for component comparison quality control tests of irradiance measurements."""
+"""Functions for component comparison quality tests of irradiance measurements."""
 
 import numpy as np
 
@@ -88,3 +88,103 @@ def diffuse_fraction_flag(
         flag = flag | np.isnan(dhi) | np.isnan(ghi)
 
     return flag
+
+
+def closure_flag(
+    ghi,
+    dni,
+    dhi,
+    solar_zenith,
+    *,
+    check="both",
+    outside_domain_flag=False,
+    nan_flag=False,
+):
+    """Flag measurements where the three-component closure ratio exceeds plausible limits.
+
+    The closure ratio R = GHI / (DHI + DNI · cos(SZA)) is compared against
+    solar-zenith-dependent limits when the component sum exceeds 50 W/m².
+
+    The limits are:
+
+    - R must be within ±8% of 1.0 for solar zenith < 75°
+    - R must be within ±15% of 1.0 for 75° ≤ solar zenith < 93°
+    - not tested for component sum ≤ 50 W/m² or solar zenith ≥ 93°
+
+    The comparison test is part of the BSRN QC tests [1]_, [2]_.
+
+    Parameters
+    ----------
+    ghi : array-like of float
+        Global horizontal irradiance [W/m²].
+    dni : array-like of float
+        Direct normal irradiance [W/m²].
+    dhi : array-like of float
+        Diffuse horizontal irradiance [W/m²].
+    solar_zenith : array-like of float
+        Solar zenith angle [degrees].
+    check : {'high-zenith', 'low-zenith', 'both'}, optional
+        Which solar zenith angle domain to check. Default is ``'both'``.
+    outside_domain_flag : bool, optional
+        Value to assign to the flag when conditions are outside the
+        valid test boundary. Can be either ``True`` or ``False``.
+        Default is ``False``, which does not flag untested values as
+        suspicious.
+    nan_flag : bool, optional
+        If ``True``, flag values where *ghi*, *dhi*, or *dni* is NaN.
+        Default is ``False``, which does not flag NaN values as suspicious.
+
+    Returns
+    -------
+    flag : same type as *ghi*
+        Boolean array. ``True`` indicates the value failed the test,
+        ``False`` indicates it passed or was outside the test domain.
+
+    See Also
+    --------
+    diffuse_fraction_flag
+
+    References
+    ----------
+    .. [1] C. N. Long and Y. Shi, "An Automated Quality Assessment and Control
+       Algorithm for Surface Radiation Measurements," *The Open Atmospheric
+       Science Journal*, vol. 2, no. 1, pp. 23–37, Apr. 2008.
+       :doi:`10.2174/1874282300802010023`
+    .. [2] `C. N. Long and E. G. Dutton, "BSRN Global Network recommended QC
+       tests, V2.0," BSRN, 2002.
+       <https://bsrn.awi.de/fileadmin/user_upload/bsrn.awi.de/Publications/BSRN_recommended_QC_tests_V2.pdf>`_
+    """
+    mu0 = np.cos(np.radians(solar_zenith))
+    sum_sw = dhi + dni * mu0
+
+    with np.errstate(divide="ignore", invalid="ignore"):
+        R = ghi / sum_sw
+
+    is_sum_50 = sum_sw > 50
+    is_low_zenith = solar_zenith < 75
+    is_high_zenith = (solar_zenith >= 75) & (solar_zenith < 93)
+
+    if check == "high-zenith":
+        flag = is_sum_50 & is_high_zenith & (np.abs(R - 1.0) >= 0.15)
+        outside_domain = np.logical_not(is_sum_50 & is_high_zenith)
+    elif check == "low-zenith":
+        flag = is_sum_50 & is_low_zenith & (np.abs(R - 1.0) >= 0.08)
+        outside_domain = np.logical_not(is_sum_50 & is_low_zenith)
+    elif check == "both":
+        flag = is_sum_50 & (
+            (is_low_zenith & (np.abs(R - 1.0) >= 0.08))
+            | (is_high_zenith & (np.abs(R - 1.0) >= 0.15))
+        )
+        outside_domain = np.logical_not(is_sum_50 & (is_low_zenith | is_high_zenith))
+    else:
+        raise ValueError(
+            f"check must be 'both', 'low-zenith', or 'high-zenith', got '{check}'."
+        )
+
+    if outside_domain_flag:
+        flag = flag | outside_domain
+
+    if nan_flag:
+        flag = flag | np.isnan(ghi) | np.isnan(dhi) | np.isnan(dni)
+
+    return flag

tests/quality/test_comparison.py

@@ -2,8 +2,7 @@
 import pandas as pd
 import pytest
 import warnings
-from solarpy.quality.comparison import diffuse_fraction_flag
-
+from solarpy.quality.comparison import diffuse_fraction_flag, closure_flag
 
 GHI = 100.0
 SZA_LOW = 45.0  # < 75°
@@ -145,3 +144,158 @@ def test_pandas_series():
     sza = pd.Series([SZA_LOW, SZA_LOW])
     flag = diffuse_fraction_flag(ghi, dhi, sza)
     assert isinstance(flag, pd.Series)
+
+
+# ===========================================================
+# closure_flag
+# ===========================================================
+
+DNI = 500.0
+DHI = 100.0
+SUM_SW_LOW = DHI + DNI * np.cos(np.radians(SZA_LOW))
+SUM_SW_HIGH = DHI + DNI * np.cos(np.radians(SZA_HIGH))
+
+
+# --- low-zenith domain, tolerance ±8% ---
+
+
+def test_closure_low_zenith_within_upper_limit_not_flagged():
+    ghi = SUM_SW_LOW * 1.07
+    assert closure_flag(ghi, DNI, DHI, SZA_LOW) == False  # noqa: E712
+
+
+def test_closure_low_zenith_at_upper_limit_flagged():
+    ghi = SUM_SW_LOW * 1.08
+    assert closure_flag(ghi, DNI, DHI, SZA_LOW) == True  # noqa: E712
+
+
+def test_closure_low_zenith_above_upper_limit_flagged():
+    ghi = SUM_SW_LOW * 1.09
+    assert closure_flag(ghi, DNI, DHI, SZA_LOW) == True  # noqa: E712
+
+
+def test_closure_low_zenith_within_lower_limit_not_flagged():
+    ghi = SUM_SW_LOW * 0.93
+    assert closure_flag(ghi, DNI, DHI, SZA_LOW) == False  # noqa: E712
+
+
+def test_closure_low_zenith_below_lower_limit_flagged():
+    ghi = SUM_SW_LOW * 0.91
+    assert closure_flag(ghi, DNI, DHI, SZA_LOW) == True  # noqa: E712
+
+
+# --- high-zenith domain, tolerance ±15% ---
+
+
+def test_closure_high_zenith_within_upper_limit_not_flagged():
+    ghi = SUM_SW_HIGH * 1.14
+    assert closure_flag(ghi, DNI, DHI, SZA_HIGH) == False  # noqa: E712
+
+
+def test_closure_high_zenith_above_upper_limit_flagged():
+    assert closure_flag(116, 0, 100, SZA_HIGH) == True  # noqa: E712
+
+
+def test_closure_high_zenith_below_lower_limit_flagged():
+    ghi = SUM_SW_HIGH * 0.84
+    assert closure_flag(ghi, DNI, DHI, SZA_HIGH) == True  # noqa: E712
+
+
+# --- outside domain ---
+
+
+def test_closure_sum_sw_below_threshold_not_flagged():
+    flag = closure_flag(ghi=40.0, dni=0.0, dhi=20.0, solar_zenith=SZA_LOW)
+    assert flag == False  # noqa: E712
+
+
+def test_closure_nighttime_not_flagged():
+    ghi = SUM_SW_LOW * 1.5
+    assert closure_flag(ghi, DNI, DHI, SZA_NIGHT) == False  # noqa: E712
+
+
+def test_closure_sum_sw_zero_no_warning():
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        flag = closure_flag(np.array(0.0), np.array(0.0), np.array(0.0), SZA_LOW)
+    assert flag == False  # noqa: E712
+
+
+# --- check parameter ---
+
+
+def test_closure_check_low_zenith_ignores_high_zenith_violation():
+    ghi = SUM_SW_HIGH * 1.16  # violation in high-zenith domain
+    flag = closure_flag(ghi, DNI, DHI, SZA_HIGH, check="low-zenith")
+    assert flag == False  # noqa: E712
+
+
+def test_closure_check_high_zenith_ignores_low_zenith_violation():
+    ghi = SUM_SW_LOW * 1.09  # violation in low-zenith domain
+    flag = closure_flag(ghi, DNI, DHI, SZA_LOW, check="high-zenith")
+    assert flag == False  # noqa: E712
+
+
+def test_closure_check_invalid_raises():
+    with pytest.raises(ValueError, match="check must be"):
+        closure_flag(GHI, DNI, DHI, SZA_LOW, check="invalid")
+
+
+# --- outside_domain_flag ---
+
+
+def test_closure_outside_domain_flag_true_flags_nighttime():
+    flag = closure_flag(GHI, DNI, DHI, SZA_NIGHT, outside_domain_flag=True)
+    assert flag == True  # noqa: E712
+
+
+def test_closure_outside_domain_flag_true_flags_low_sum_sw():
+    flag = closure_flag(40.0, 0.0, 20.0, SZA_LOW, outside_domain_flag=True)
+    assert flag == True  # noqa: E712
+
+
+def test_closure_outside_domain_flag_false_does_not_flag_nighttime():
+    flag = closure_flag(GHI, DNI, DHI, SZA_NIGHT, outside_domain_flag=False)
+    assert flag == False  # noqa: E712
+
+
+# --- nan_flag ---
+
+
+def test_closure_nan_not_flagged_by_default():
+    assert closure_flag(float("nan"), DNI, DHI, SZA_LOW) == False  # noqa: E712
+    assert closure_flag(GHI, float("nan"), DHI, SZA_LOW) == False  # noqa: E712
+    assert closure_flag(GHI, DNI, float("nan"), SZA_LOW) == False  # noqa: E712
+
+
+def test_closure_nan_flagged_when_nan_flag_true():
+    ghi_flag = closure_flag(np.nan, DNI, DHI, SZA_LOW, nan_flag=True)
+    assert ghi_flag == True  # noqa: E712
+    dni_flag = closure_flag(GHI, np.nan, DHI, SZA_LOW, nan_flag=True)
+    assert dni_flag == True  # noqa: E712
+    dhi_flag = closure_flag(GHI, DNI, np.nan, SZA_LOW, nan_flag=True)
+    assert dhi_flag == True  # noqa: E712
+
+
+# --- array inputs ---
+
+
+def test_closure_numpy_array():
+    ghi = np.array([SUM_SW_LOW * 1.07, SUM_SW_LOW * 1.09, 40.0, SUM_SW_LOW * 0.91])
+    dni = np.array([DNI, DNI, 0.0, DNI])
+    dhi = np.array([DHI, DHI, 20.0, DHI])
+    sza = np.array([SZA_LOW, SZA_LOW, SZA_LOW, SZA_LOW])
+    flag = closure_flag(ghi, dni, dhi, sza)
+    assert flag[0] == False  # noqa: E712  within ±8%
+    assert flag[1] == True  # noqa: E712  above +8%
+    assert flag[2] == False  # noqa: E712  sum_sw ≤ 50
+    assert flag[3] == True  # noqa: E712  below -8%
+
+
+def test_closure_pandas_series():
+    ghi = pd.Series([SUM_SW_LOW * 1.07, SUM_SW_LOW * 1.09])
+    dni = pd.Series([DNI, DNI])
+    dhi = pd.Series([DHI, DHI])
+    sza = pd.Series([SZA_LOW, SZA_LOW])
+    flag = closure_flag(ghi, dni, dhi, sza)
+    assert isinstance(flag, pd.Series)