Add BSRN limit checks

docs/source/conf.py

@@ -29,11 +29,16 @@
     'myst_nb',  # markdown and jupyter-notebook parsing
     'sphinx.ext.autodoc',  # generate documentation from docstrings
     'sphinx.ext.autosummary',
+    'sphinx.ext.extlinks',  # short external link roles (e.g. :doi:)
     'sphinx.ext.intersphinx',
     'sphinx.ext.napoleon',  # parsing of Numpy docstrings
     'sphinx_copybutton',  # add copy button to code examples
 ]
 
+extlinks = {
+    'doi': ('https://doi.org/%s', 'doi:%s'),
+}
+
 # Configure copybutton to strip >>> when copying code
 copybutton_prompt_text = r">>> |\.\.\. "
 copybutton_prompt_is_regexp = True

docs/source/documentation.rst

@@ -13,4 +13,6 @@ Code documentation
    plotting.plot_intraday_heatmap
    plotting.plot_shading_heatmap
    horizon.get_horizon_mines
+   quality.bsrn_limits
+   quality.bsrn_limits_flag
    iotools.read_t16

src/solarpy/quality/__init__.py

@@ -1 +1,2 @@
-from solarpy.quality.limits import upper_limit_bsrn  # noqa: F401
+from solarpy.quality.limits import bsrn_limits  # noqa: F401
+from solarpy.quality.limits import bsrn_limits_flag  # noqa: F401

src/solarpy/quality/limits.py

@@ -1,4 +1,191 @@
+"""Functions for BSRN quality control irradiance limit tests."""
 
+import numpy as np
 
-def upper_limit_bsrn():
-    return 1
+
+_BSRN_LIMITS = {
+    "ppl-ghi": (1.50, 1.2,  100,   -4.),
+    "erl-ghi": (1.20, 1.2,   50,   -2.),
+    "ppl-dni": (1.00, 1.0,    0,   -4.),
+    "erl-dni": (0.95, 0.2,   10,   -2.),
+    "ppl-dhi": (0.95, 1.2,   50,   -4.),
+    "erl-dhi": (0.75, 1.2,   30,   -2.),
+}
+
+
+def bsrn_limits(solar_zenith, dni_extra, limits):
+    """Calculate the BSRN upper and/or lower irradiance limit values.
+
+    The BSRN upper and lower bound limit checks were developed by Long & Shi
+    (2008) [1]_, [2]_. The upper limit follows the form::
+
+        upper = a * DNI_extra * cos(solar_zenith) ^ b + c
+
+    where *a*, *b*, and *c* are coefficients that depend on the variable
+    and test level. A value is flagged if it lies outside [lower, upper].
+
+    Parameters
+    ----------
+    solar_zenith : array-like of float
+        Solar zenith angle [degrees].
+    dni_extra : array-like of float
+        Extraterrestrial normal irradiance [W/m²].
+    limits : str or tuple of float
+        Either a named limit string or a tuple ``(a, b, c, lower)``.
+
+        Named limit (Long & Shi, 2008):
+
+        - ``"ppl-ghi"`` — Physically Possible Limit for GHI
+        - ``"erl-ghi"`` — Extremely Rare Limit for GHI
+        - ``"ppl-dni"`` — Physically Possible Limit for DNI
+        - ``"erl-dni"`` — Extremely Rare Limit for DNI
+        - ``"ppl-dhi"`` — Physically Possible Limit for DHI
+        - ``"erl-dhi"`` — Extremely Rare Limit for DHI
+
+        When passing a tuple, provide ``(a, b, c, lower)`` where the upper
+        bound is ``a * dni_extra * cos(solar_zenith) ** b + c`` and *lower*
+        is the minimum allowed value.
+
+    Returns
+    -------
+    lower : float
+        Lower limit value [W/m²].
+    upper : same type as input
+        Upper limit values [W/m²].
+
+    See Also
+    --------
+    bsrn_limits_flag : Test irradiance values against these limits.
+
+    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 Y. Shi, "An Automated Quality Assessment and Control
+       Algorithm for Surface Radiation Measurements," BSRN, 2008. [Online].
+       Available: `BSRN recommended QC tests v2
+       <https://bsrn.awi.de/fileadmin/user_upload/bsrn.awi.de/Publications/BSRN_recommended_QC_tests_V2.pdf>`_
+    """
+    if isinstance(limits, str):
+        if limits not in _BSRN_LIMITS:
+            raise ValueError(
+                f"Unknown limit '{limits}'. "
+                f"Valid options are: {list(_BSRN_LIMITS.keys())}."
+            )
+        a, b, c, lower = _BSRN_LIMITS[limits]
+    elif isinstance(limits, tuple):
+        if len(limits) != 4:
+            raise ValueError(
+                f"limit tuple must have 4 elements (a, b, c, lower), "
+                f"got {len(limits)}."
+            )
+        a, b, c, lower = limits
+    else:
+        raise ValueError("limit must be a string or a tuple of 4 floats.")
+
+    cos_sza = np.cos(np.deg2rad(solar_zenith))
+    upper = a * dni_extra * cos_sza ** b + c
+
+    return lower, upper
+
+
+def bsrn_limits_flag(irradiance, solar_zenith, dni_extra, limits, check='both', nan_flag=True):
+    """Flag irradiance values that fall outside the BSRN quality control limits.
+
+    Parameters
+    ----------
+    irradiance : array-like of float
+        Irradiance values to check [W/m²].
+    solar_zenith : array-like of float
+        Solar zenith angle [degrees]. Must be the same length as *irradiance*.
+    dni_extra : array-like of float
+        Extraterrestrial normal irradiance [W/m²]. Must be the same length
+        as *irradiance*.
+    limits : str or tuple of float
+        Either a named limit string or a tuple of coefficients
+        ``(a, b, c, lower)``.
+
+        Named limit (Long & Shi, 2008) [1]_, [2]_:
+
+        - ``"ppl-ghi"`` — Physically Possible Limit for GHI
+        - ``"erl-ghi"`` — Extremely Rare Limit for GHI
+        - ``"ppl-dni"`` — Physically Possible Limit for DNI
+        - ``"erl-dni"`` — Extremely Rare Limit for DNI
+        - ``"ppl-dhi"`` — Physically Possible Limit for DHI
+        - ``"erl-dhi"`` — Extremely Rare Limit for DHI
+
+        When passing a tuple, provide ``(a, b, c, lower)`` where the upper
+        bound is ``a * dni_extra * cos(solar_zenith) ** b + c`` and *lower* is the
+        minimum allowed value.
+    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 ``True``,
+        which flags NaN values as suspicious.
+
+    Returns
+    -------
+    flag : same type as *irradiance*
+        Boolean array of the same length as *irradiance*. ``True`` indicates
+        the value failed the test (outside bounds), ``False`` indicates
+        it passed.
+
+    See Also
+    --------
+    bsrn_limits : Calculate the limit values without testing.
+
+    Examples
+    --------
+    Test GHI measurements against the BSRN limits:
+
+    >>> import pandas as pd
+    >>> import numpy as np
+    >>> import pvlib
+    >>>
+    >>> # One year of hourly timestamps for Copenhagen
+    >>> times = pd.date_range("2023-01-01", periods=8760, freq="h", tz="UTC")
+    >>> latitude, longitude = 55.68, 12.57
+    >>>
+    >>> # Calculate solar position and extraterrestrial irradiance using pvlib
+    >>> solpos = pvlib.solarposition.get_solarposition(times, latitude, longitude)
+    >>> solar_zenith = solpos["apparent_zenith"]
+    >>> dni_extra = pvlib.irradiance.get_extra_radiation(times)
+    >>>
+    >>> # Create synthetic GHI: sine wave clipped to daytime
+    >>> rng = np.random.default_rng(seed=0)
+    >>> cos_sza = np.cos(np.deg2rad(solar_zenith))
+    >>> ghi = np.clip(900 * cos_sza + rng.standard_normal(8760) * 20, 0, None)
+    >>>
+    >>> # Run PPL and ERL tests
+    >>> ppl_flag = bsrn_limits_flag(ghi, solar_zenith, dni_extra, limits="ppl-ghi")
+    >>> erl_flag = bsrn_limits_flag(ghi, solar_zenith, dni_extra, limits="erl-ghi")
+
+    Use custom coefficients:
+
+    >>> flag = bsrn_limits_flag(ghi, solar_zenith, dni_extra, limits=(1.2, 1.2, 50, -4))
+
+    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 Y. Shi, "An Automated Quality Assessment and Control
+       Algorithm for Surface Radiation Measurements," BSRN, 2008. [Online].
+       Available: `BSRN recommended QC tests v2
+       <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)
+    if check == 'upper':
+        flag = irradiance > upper
+    elif check == 'lower':
+        flag = irradiance < lower
+    elif check == 'both':
+        flag = (irradiance < lower) | (irradiance > upper)
+    else:
+        raise ValueError(f"check must be 'both', 'upper', or 'lower', got '{check}'.")
+    if nan_flag:
+        flag = flag | np.isnan(irradiance)
+    return flag

tests/quality/test_limits.py

@@ -0,0 +1,153 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+from solarpy.quality.limits import bsrn_limits, bsrn_limits_flag
+
+
+SZA = 45.0          # solar zenith angle [degrees]
+DNI_EXTRA = 1361.0  # extraterrestrial normal irradiance [W/m²]
+COS_SZA = np.cos(np.deg2rad(SZA))
+
+
+# bsrn_limit — float inputs
+
+def test_bsrn_limit_returns_tuple_of_two():
+    lower, upper = bsrn_limits(SZA, DNI_EXTRA, "ppl-ghi")
+    assert isinstance(lower, float)
+    assert isinstance(upper, float)
+
+
+def test_bsrn_limit_ppl_ghi():
+    lower, upper = bsrn_limits(SZA, DNI_EXTRA, "ppl-ghi")
+    expected = 1.50 * DNI_EXTRA * COS_SZA ** 1.2 + 100
+    assert upper == pytest.approx(expected)
+    assert lower == -4
+
+
+def test_bsrn_limit_erl_ghi():
+    _, upper = bsrn_limits(SZA, DNI_EXTRA, "erl-ghi")
+    expected = 1.20 * DNI_EXTRA * COS_SZA ** 1.2 + 50
+    assert upper == pytest.approx(expected)
+
+
+def test_bsrn_limit_ppl_dni_upper():
+    _, upper = bsrn_limits(SZA, DNI_EXTRA, "ppl-dni")
+    expected = 1.00 * DNI_EXTRA * COS_SZA ** 1.0 + 0
+    assert upper == pytest.approx(expected)
+
+
+def test_bsrn_limit_ppl_dhi_upper():
+    _, upper = bsrn_limits(SZA, DNI_EXTRA, "ppl-dhi")
+    expected = 0.95 * DNI_EXTRA * COS_SZA ** 1.2 + 50
+    assert upper == expected
+
+
+def test_bsrn_limit_zenith_90_upper_is_c():
+    _, upper = bsrn_limits(90.0, DNI_EXTRA, "ppl-ghi")
+    assert upper == 100.0
+
+
+def test_bsrn_limit_custom_tuple():
+    lower, upper = bsrn_limits(SZA, DNI_EXTRA, (1.0, 1.0, 99, -99))
+    assert lower == -99
+    assert upper == pytest.approx(DNI_EXTRA * COS_SZA + 99)
+
+
+def test_bsrn_limit_invalid_string_raises():
+    with pytest.raises(ValueError, match="Unknown limit"):
+        bsrn_limits(SZA, DNI_EXTRA, "invalid")
+
+
+def test_bsrn_limit_invalid_tuple_length_raises():
+    with pytest.raises(ValueError, match="4 elements"):
+        bsrn_limits(SZA, DNI_EXTRA, (1.0, 1.0, 50))
+
+
+def test_bsrn_limit_invalid_type_raises():
+    with pytest.raises(ValueError):
+        bsrn_limits(SZA, DNI_EXTRA, 42)
+
+
+# bsrn_limit — array inputs
+
+def test_bsrn_limit_numpy_array():
+    sza = np.array([0.0, 45.0, 90.0])
+    _, upper = bsrn_limits(sza, DNI_EXTRA, "ppl-ghi")
+    assert upper.shape == (3,)
+
+
+def test_bsrn_limit_pandas_series():
+    sza = pd.Series([0.0, 45.0, 90.0])
+    _, upper = bsrn_limits(sza, DNI_EXTRA, "ppl-ghi")
+    assert isinstance(upper, pd.Series)
+    assert len(upper) == 3
+
+
+# bsrn_limits_flag — float inputs
+
+def test_bsrn_limits_flag_value_within_bounds_not_flagged():
+    lower, upper = bsrn_limits(SZA, DNI_EXTRA, "ppl-ghi")
+    mid = (lower + upper) / 2
+    assert bsrn_limits_flag(mid, SZA, DNI_EXTRA, "ppl-ghi") == False  # noqa: E712
+
+
+def test_bsrn_limits_flag_value_above_upper_flagged():
+    _, upper = bsrn_limits(SZA, DNI_EXTRA, "ppl-ghi")
+    assert bsrn_limits_flag(upper + 1, SZA, DNI_EXTRA, "ppl-ghi") == True  # noqa: E712
+
+
+def test_bsrn_limits_flag_value_below_lower_flagged():
+    lower, _ = bsrn_limits(SZA, DNI_EXTRA, "ppl-ghi")
+    assert bsrn_limits_flag(lower - 1, SZA, DNI_EXTRA, "ppl-ghi") == True  # noqa: E712
+
+
+def test_bsrn_limits_flag_check_upper_only():
+    lower, _ = bsrn_limits(SZA, DNI_EXTRA, "ppl-ghi")
+    # value below lower but only checking upper — should not be flagged
+    result = bsrn_limits_flag(lower - 1, SZA, DNI_EXTRA, "ppl-ghi", check='upper')
+    assert result == False  # noqa: E712
+
+
+def test_bsrn_limits_flag_check_lower_only():
+    _, upper = bsrn_limits(SZA, DNI_EXTRA, "ppl-ghi")
+    # value above upper but only checking lower — should not be flagged
+    result = bsrn_limits_flag(upper + 1, SZA, DNI_EXTRA, "ppl-ghi", check='lower')
+    assert result == False  # noqa: E712
+
+
+def test_bsrn_limits_flag_nan_flagged_by_default():
+    assert bsrn_limits_flag(np.nan, SZA, DNI_EXTRA, "ppl-ghi") == True  # noqa: E712
+
+
+def test_bsrn_limits_flag_nan_not_flagged_when_nan_flag_false():
+    result = bsrn_limits_flag(np.nan, SZA, DNI_EXTRA, "ppl-ghi", nan_flag=False)
+    assert result == False  # noqa: E712
+
+
+def test_bsrn_limits_flag_invalid_check_raises():
+    with pytest.raises(ValueError, match="check must be"):
+        bsrn_limits_flag(500.0, SZA, DNI_EXTRA, "ppl-ghi", check='invalid')
+
+
+# bsrn_limits_flag — array inputs
+
+def test_bsrn_limits_flag_numpy_array():
+    ghi = np.array([-10.0, 500.0, 9999.0])
+    flag = bsrn_limits_flag(ghi, SZA, DNI_EXTRA, "ppl-ghi")
+    assert flag[0] == True   # below lower  # noqa: E712
+    assert flag[1] == False  # within bounds  # noqa: E712
+    assert flag[2] == True   # above upper  # noqa: E712
+
+
+def test_bsrn_limits_flag_pandas_series():
+    ghi = pd.Series([-10.0, 500.0, 9999.0])
+    flag = bsrn_limits_flag(ghi, SZA, DNI_EXTRA, "ppl-ghi")
+    assert isinstance(flag, pd.Series)
+
+
+def test_bsrn_limits_flag_nan_in_array_flagged():
+    ghi = np.array([500.0, np.nan])
+    flag = bsrn_limits_flag(ghi, SZA, DNI_EXTRA, "ppl-ghi")
+    assert flag[0] == False  # noqa: E712
+    assert flag[1] == True   # noqa: E712