Add BSRN diffuse fraction test

docs/source/documentation.rst

@@ -15,4 +15,5 @@ Code documentation
    horizon.get_horizon_mines
    quality.bsrn_limits
    quality.bsrn_limits_flag
+   quality.diffuse_fraction_flag
    iotools.read_t16

src/solarpy/quality/__init__.py

@@ -1,2 +1,3 @@
 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

src/solarpy/quality/comparison.py

@@ -0,0 +1,90 @@
+"""Functions for component comparison quality control tests of irradiance measurements."""
+
+import numpy as np
+
+
+def diffuse_fraction_flag(
+    ghi, dhi, solar_zenith, *, check="both", outside_domain_flag=False, nan_flag=False
+):
+    """Flag measurements where the diffuse fraction exceeds physically plausible limits.
+
+    The diffuse fraction K = DHI / GHI is tested against solar-zenith-dependent
+    upper limits when GHI exceeds 50 W/m². The limits are:
+
+    - K must be < 1.05 for solar zenith < 75°
+    - K must be < 1.10 for 75° ≤ solar zenith < 93°
+    - not tested for GHI ≤ 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²].
+    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* or *dhi* 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
+    --------
+    bsrn_limits_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>`_
+    """
+    # Suppress divide-by-zero warning
+    with np.errstate(divide="ignore", invalid="ignore"):
+        K = dhi / ghi
+    # TODO: Consideer adding option to also test for (dhi > 50) | (ghi > 50)
+    is_ghi_50 = ghi > 50
+
+    is_low_zenith = solar_zenith < 75
+    is_high_zenith = (solar_zenith >= 75) & (solar_zenith < 93)
+
+    if check == "high-zenith":
+        flag = is_ghi_50 & is_high_zenith & (K >= 1.10)
+        outside_domain = np.logical_not(is_ghi_50 & is_high_zenith)
+    elif check == "low-zenith":
+        flag = is_ghi_50 & is_low_zenith & (K >= 1.05)
+        outside_domain = np.logical_not(is_ghi_50 & is_low_zenith)
+    elif check == "both":
+        flag = is_ghi_50 & (is_low_zenith & (K >= 1.05)) | (
+            is_high_zenith & (K >= 1.10)
+        )
+        outside_domain = np.logical_not(is_ghi_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(dhi) | np.isnan(ghi)
+
+    return flag

src/solarpy/quality/limits.py

@@ -93,7 +93,7 @@ def bsrn_limits(solar_zenith, dni_extra, limits):
     return lower, upper
 
 
-def bsrn_limits_flag(irradiance, solar_zenith, dni_extra, limits, check='both', nan_flag=False):
+def bsrn_limits_flag(irradiance, solar_zenith, dni_extra, limits, *, check='both', nan_flag=False):
     """Flag irradiance values that fall outside the BSRN quality control limits.
 
     Parameters
@@ -125,8 +125,9 @@ def bsrn_limits_flag(irradiance, solar_zenith, dni_extra, limits, check='both',
     check : {'both', 'upper', 'lower'}, optional
         Which bounds to check. Default is ``'both'``.
     nan_flag : bool, optional
-        Flag value to assign when *irradiance* is NaN. Default is ``False``,
-        which does not flag NaN values as suspicious.
+        Flag value to assign when *irradiance* is NaN. Value can be either
+        ``True`` or ``False``. Default is ``False``, which does not flag
+        NaN values as suspicious.
 
     Returns
     -------
@@ -138,6 +139,7 @@ def bsrn_limits_flag(irradiance, solar_zenith, dni_extra, limits, check='both',
     See Also
     --------
     bsrn_limits : Calculate the limit values without testing.
+    diffuse_fraction_flag : Flag measurements based on the diffuse fraction.
 
     Examples
     --------
@@ -175,9 +177,8 @@ def bsrn_limits_flag(irradiance, solar_zenith, dni_extra, limits, check='both',
        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 Y. Shi, "An Automated Quality Assessment and Control
-       Algorithm for Surface Radiation Measurements," BSRN, 2002. [Online].
-       Available: `BSRN recommended QC tests v2
+    .. [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>`_
     """
     lower, upper = bsrn_limits(solar_zenith, dni_extra, limits)

tests/quality/test_comparison.py

@@ -0,0 +1,147 @@
+import numpy as np
+import pandas as pd
+import pytest
+import warnings
+from solarpy.quality.comparison import diffuse_fraction_flag
+
+
+GHI = 100.0
+SZA_LOW = 45.0  # < 75°
+SZA_HIGH = 80.0  # 75° ≤ SZA < 93°
+SZA_NIGHT = 95.0  # ≥ 93°
+
+
+# --- low-zenith domain (SZA < 75) ---
+
+
+def test_low_zenith_within_limit_not_flagged():
+    dhi = GHI * 1.04  # K = 1.04 < 1.05
+    assert diffuse_fraction_flag(GHI, dhi, SZA_LOW) == False  # noqa: E712
+
+
+def test_low_zenith_at_limit_flagged():
+    dhi = GHI * 1.05  # K = 1.05, >= threshold
+    assert diffuse_fraction_flag(GHI, dhi, SZA_LOW) == True  # noqa: E712
+
+
+def test_low_zenith_above_limit_flagged():
+    dhi = GHI * 1.06  # K = 1.06 > 1.05
+    assert diffuse_fraction_flag(GHI, dhi, SZA_LOW) == True  # noqa: E712
+
+
+# --- high-zenith domain (75 ≤ SZA < 93) ---
+
+
+def test_high_zenith_within_limit_not_flagged():
+    dhi = GHI * 1.09  # K = 1.09 < 1.10
+    assert diffuse_fraction_flag(GHI, dhi, SZA_HIGH) == False  # noqa: E712
+
+
+def test_high_zenith_at_limit_flagged():
+    dhi = GHI * 1.10  # K = 1.10, >= threshold
+    assert diffuse_fraction_flag(GHI, dhi, SZA_HIGH) == True  # noqa: E712
+
+
+def test_high_zenith_above_limit_flagged():
+    dhi = GHI * 1.11  # K = 1.11 > 1.10
+    assert diffuse_fraction_flag(GHI, dhi, SZA_HIGH) == True  # noqa: E712
+
+
+# --- outside domain ---
+
+
+def test_ghi_below_threshold_not_flagged():
+    assert diffuse_fraction_flag(40.0, 200.0, SZA_LOW) == False  # noqa: E712
+
+
+def test_nighttime_not_flagged():
+    dhi = GHI * 1.5  # K far above limit, but SZA ≥ 93
+    assert diffuse_fraction_flag(GHI, dhi, SZA_NIGHT) == False  # noqa: E712
+
+
+def test_ghi_zero_no_warning():
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        flag = diffuse_fraction_flag(np.array(0.0), np.array(50.0), SZA_LOW)
+    assert flag == False  # noqa: E712
+
+
+# --- check parameter ---
+
+
+def test_check_low_zenith_ignores_high_zenith_violation():
+    dhi = GHI * 1.11  # K > 1.10, violation in high-zenith
+    flag = diffuse_fraction_flag(GHI, dhi, SZA_HIGH, check="low-zenith")
+    assert flag == False  # noqa: E712
+
+
+def test_check_high_zenith_ignores_low_zenith_violation():
+    dhi = GHI * 1.06  # K >= 1.05, violation in low-zenith
+    flag = diffuse_fraction_flag(GHI, dhi, SZA_LOW, check="high-zenith")
+    assert flag == False  # noqa: E712
+
+
+def test_check_invalid_raises():
+    with pytest.raises(ValueError, match="check must be"):
+        diffuse_fraction_flag(GHI, GHI, SZA_LOW, check="invalid")
+
+
+# --- outside_domain_flag ---
+
+
+def test_outside_domain_flag_true_flags_nighttime():
+    flag = diffuse_fraction_flag(GHI, GHI, SZA_NIGHT, outside_domain_flag=True)
+    assert flag == True  # noqa: E712
+
+
+def test_outside_domain_flag_true_flags_low_ghi():
+    flag = diffuse_fraction_flag(40, 40, SZA_LOW, outside_domain_flag=True)
+    assert flag == True  # noqa: E712
+
+
+def test_outside_domain_flag_false_does_not_flag_nighttime():
+    flag = diffuse_fraction_flag(GHI, GHI, SZA_NIGHT, outside_domain_flag=False)
+    assert flag == False  # noqa: E712
+
+
+# --- nan_flag ---
+
+
+def test_nan_ghi_not_flagged_by_default():
+    assert diffuse_fraction_flag(float("nan"), GHI, SZA_LOW) == False  # noqa: E712
+
+
+def test_nan_dhi_not_flagged_by_default():
+    assert diffuse_fraction_flag(GHI, float("nan"), SZA_LOW) == False  # noqa: E712
+
+
+def test_nan_ghi_flagged_when_nan_flag_true():
+    flag = diffuse_fraction_flag(np.nan, GHI, SZA_LOW, nan_flag=True)
+    assert flag == True  # noqa: E712
+
+
+def test_nan_dhi_flagged_when_nan_flag_true():
+    flag = diffuse_fraction_flag(GHI, float("nan"), SZA_LOW, nan_flag=True)
+    assert flag == True  # noqa: E712
+
+
+# --- array inputs ---
+
+
+def test_numpy_array():
+    ghi = np.array([100.0, 100.0, 40.0, 55.0])
+    dhi = np.array([104.0, 106.0, 100.0, 100.0])
+    sza = np.array([SZA_LOW, SZA_LOW, SZA_LOW, SZA_LOW])
+    flag = diffuse_fraction_flag(ghi, dhi, sza)
+    assert flag[0] == False  # noqa: E712
+    assert flag[1] == True  # noqa: E712
+    assert flag[2] == False  # noqa: E712
+    assert flag[3] == True  # noqa: E712
+
+
+def test_pandas_series():
+    ghi = pd.Series([100.0, 100.0])
+    dhi = pd.Series([104.0, 106.0])
+    sza = pd.Series([SZA_LOW, SZA_LOW])
+    flag = diffuse_fraction_flag(ghi, dhi, sza)
+    assert isinstance(flag, pd.Series)