diff --git a/docs/source/documentation.rst b/docs/source/documentation.rst index 723f92e..1d2520c 100644 --- a/docs/source/documentation.rst +++ b/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 diff --git a/src/solarpy/quality/__init__.py b/src/solarpy/quality/__init__.py index 159b2d9..7e49251 100644 --- a/src/solarpy/quality/__init__.py +++ b/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 diff --git a/src/solarpy/quality/comparison.py b/src/solarpy/quality/comparison.py index 5fbf2cc..a3cd7aa 100644 --- a/src/solarpy/quality/comparison.py +++ b/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. + `_ + """ + 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 diff --git a/tests/quality/test_comparison.py b/tests/quality/test_comparison.py index 72c7670..d14f137 100644 --- a/tests/quality/test_comparison.py +++ b/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)