Skip to content

Commit

Permalink
Adding radiation quality control tests (#874)
Browse files Browse the repository at this point in the history 
* Adding Closure test from BSRN. Also allows test_options to be case insentive.

* Adding other component derivations and standardized variable names to not require shortwave.

* Adding tests for new retrievals and QC tests.

* Adding K tests

* Adding K tests, test, tests.

* Updating limits

* Adding K-tests

* Updated assert for new test limit.

* Renaming k_test to not be confused with other statistical tests available in numpy or other libraries.
  • Loading branch information
@kenkehoe
kenkehoe authored Nov 7, 2024
1 parent b07be1e commit 6adb9be
Show file tree
Hide file tree
Showing 5 changed files with 710 additions and 34 deletions.
310 changes: 300 additions & 10 deletions act/qc/bsrn_tests.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
from act.utils.data_utils import convert_units


def _calculate_solar_parameters(ds, lat_name, lon_name, solar_constant):
def _calculate_solar_parameters(ds, lat_name, lon_name, solar_constant=1360.8):
"""
Function to calculate solar zenith angles and solar constant adjusted
to Earth Sun distance
Expand Down Expand Up @@ -128,7 +128,7 @@ def bsrn_limits_test(
Parameters
----------
test : str
test : str or list
Type of tests to apply. Options include "Physically Possible" or "Extremely Rare"
gbl_SW_dn_name : str
Variable name in the Dataset for global shortwave downwelling radiation
Expand Down Expand Up @@ -415,9 +415,9 @@ def bsrn_comparison_tests(
Parameters
----------
test : str
test : str or list
Type of tests to apply. Options include: 'Global over Sum SW Ratio', 'Diffuse Ratio',
'SW up', 'LW down to air temp', 'LW up to air temp', 'LW down to LW up'
'SW up', 'LW down to air temp', 'LW up to air temp', 'LW down to LW up', 'Closure'
gbl_SW_dn_name : str
Variable name in Dataset for global shortwave downwelling radiation
measured by unshaded pyranometer
Expand Down Expand Up @@ -451,6 +451,8 @@ def bsrn_comparison_tests(
References
----------
Long, Charles N., and Ellsworth G. Dutton. "BSRN Global Network recommended QC tests, V2. x." (2010).
Long, Charles N., and Y. Shi. “An Automated Quality Assessment and Control Algorithm for Surface
Radiation Measurements.” (2008) https://doi.org/10.2174/1874282300802010023
Examples
--------
Expand All @@ -470,20 +472,23 @@ def bsrn_comparison_tests(
if isinstance(test, str):
test = [test]

test = [ii.lower() for ii in test]

test_options = [
'Global over Sum SW Ratio',
'Diffuse Ratio',
'SW up',
'LW down to air temp',
'LW up to air temp',
'LW down to LW up',
'Closure',
]

solar_constant = 1360.8
sza, Sa = _calculate_solar_parameters(self._ds, lat_name, lon_name, solar_constant)

# Ratio of Global over Sum SW
if test_options[0] in test:
if test_options[0].lower() in test:
if (
gbl_SW_dn_name is None
or glb_diffuse_SW_dn_name is None
Expand Down Expand Up @@ -544,7 +549,7 @@ def bsrn_comparison_tests(
)

# Diffuse Ratio
if test_options[1] in test:
if test_options[1].lower() in test:
if gbl_SW_dn_name is None or glb_diffuse_SW_dn_name is None:
raise ValueError(
'Must set keywords gbl_SW_dn_name, glb_diffuse_SW_dn_name '
Expand Down Expand Up @@ -587,7 +592,7 @@ def bsrn_comparison_tests(
)

# Shortwave up comparison
if test_options[2] in test:
if test_options[2].lower() in test:
if (
glb_SW_up_name is None
or glb_diffuse_SW_dn_name is None
Expand Down Expand Up @@ -637,7 +642,7 @@ def bsrn_comparison_tests(
)

# Longwave down to air temperature comparison
if test_options[3] in test:
if test_options[3].lower() in test:
if glb_LW_dn_name is None or air_temp_name is None:
raise ValueError(
'Must set keywords glb_LW_dn_name, air_temp_name '
Expand Down Expand Up @@ -670,7 +675,7 @@ def bsrn_comparison_tests(
)

# Longwave up to air temperature comparison
if test_options[4] in test:
if test_options[4].lower() in test:
if glb_LW_up_name is None or air_temp_name is None:
raise ValueError(
'Must set keywords glb_LW_up_name, air_temp_name '
Expand Down Expand Up @@ -705,7 +710,7 @@ def bsrn_comparison_tests(
)

# Lonwave down to longwave up comparison
if test_options[5] in test:
if test_options[5].lower() in test:
if glb_LW_dn_name is None or glb_LW_up_name is None:
raise ValueError(
'Must set keywords glb_LW_dn_name, glb_LW_up_name '
Expand Down Expand Up @@ -750,3 +755,288 @@ def bsrn_comparison_tests(
test_assessment=test_assessment,
test_meaning=test_meaning,
)

# Closure test
if test_options[6].lower() in test:
if (
direct_normal_SW_dn_name is None
or glb_diffuse_SW_dn_name is None
or gbl_SW_dn_name is None
):
raise ValueError(
'Must set keywords direct_normal_SW_dn_name, glb_diffuse_SW_dn_name, '
f'gbl_SW_dn_name for {test_options[3]} test.'
)

if use_dask and isinstance(self._ds[gbl_SW_dn_name].data, da.Array):
sza = da.array(sza)
component_sum = self._ds[gbl_SW_dn_name].values + (
self._ds[direct_normal_SW_dn_name].values * np.cos(np.radians(sza))
)

index_lsz = da.where((component_sum > 50.0) & (sza <= 75.0), True, False)
index_hsz = da.where(
(component_sum > 50.0) & (sza > 75.0) & (sza < 93.0), True, False
)

value = (self._ds[glb_diffuse_SW_dn_name].values / component_sum) - 1.0

index_1 = da.where((value >= 0.08) & index_lsz, True, False)
index_2 = da.where((value >= 0.15) & index_hsz, True, False)

index = (index_1 | index_2).compute()

else:
component_sum = self._ds[gbl_SW_dn_name].values + (
self._ds[direct_normal_SW_dn_name].values * np.cos(np.radians(sza))
)

index_lsz = (component_sum > 50.0) & (sza <= 75.0)
index_hsz = (component_sum > 50.0) & (sza > 75.0) & (sza < 93.0)

value = (self._ds[glb_diffuse_SW_dn_name].values / component_sum) - 1.0

index_1 = (value >= 0.08) & index_lsz
index_2 = (value >= 0.15) & index_hsz

index = index_1 | index_2

test_meaning = "Closure test indicating value outside of expected range"
self._ds.qcfilter.add_test(
direct_normal_SW_dn_name,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)
self._ds.qcfilter.add_test(
glb_diffuse_SW_dn_name,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)
self._ds.qcfilter.add_test(
gbl_SW_dn_name,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)

def normalized_rradiance_test(
self,
test,
dni=None,
dhi=None,
ghi=None,
test_assessment='Indeterminate',
lat_name='lat',
lon_name='lon',
alt_name='alt',
upper_total_transmittance_limit=1.6,
upper_diffuse_transmittance_limit=1.0,
use_dask=False,
):
"""
Method to apply Normalized Irradiance Tests tests and add results to ancillary quality control variable.
Need to provided variable name for each measurement for the test to be performed. All radiation
data must be in W/m^2 units. All results from tests are set in ancillary quality control varibles
for dni and dhi.
test : list or str
Type of tests to apply. Options include: 'Clearness index', 'Upper total transmittance',
'Upper direct transmittance', 'Upper diffuse transmittance'
dni : str
Variable name in Dataset for direct normal irradiance
dhi : str
Variable name in Dataset for diffuse hemipheric irradiance
ghi : str
Variable name in Dataset for global hemispheric irradiance. Used in
'Upper direct transmittance' test only.
test_assessment : str
Test assessment string value appended to flag_assessments attribute of QC variable.
lat_name : str
Variable name in the Dataset for latitude
lon_name : str
Variable name in the Dataset for longitude
alt_name : str
Variable name in the Dataset for altitude. Used in 'Upper direct transmittance' test.
upper_total_transmittance_limit : float
Limit value used in 'Upper total transmittance' test. Values larger than this will be flagged.
upper_diffuse_transmittance_limit : float
Limit value used in 'Upper diffuse transmittance' test. Values larger than this will be flagged.
use_dask : boolean
Option to use Dask for processing if data is stored in a Dask array
References
----------
Sengupta, Manajit, Habte, Aron, Wilbert, Stefan , Christian Gueymard, Jan Remund, Elke Lorenz,
Wilfried van Sark, and Adam R. Jensen "Best Practices Handbook for the Collection and Use of Solar
Resource Data for Solar Energy Applications: Fourth Edition" (2024).
Examples
--------
.. code-block:: python
ds = act.io.arm.read_arm_netcdf(act.tests.EXAMPLE_BRS, cleanup_qc=True)
ds.qcfilter.normalized_rradiance_test(
['Clearness index', 'Upper total transmittance',
'Upper direct transmittance', 'Upper diffuse transmittance'],
dhi='down_short_diffuse_hemisp',
ghi='down_short_hemisp',
dni='short_direct_normal'
)
"""

if isinstance(test, str):
test = [test]

test = [ii.lower() for ii in test]

solar_constant = 1366
sza, Sa = _calculate_solar_parameters(self._ds, lat_name, lon_name, solar_constant)

test_options = [
'Clearness index',
'Upper total transmittance',
'Upper direct transmittance',
'Upper diffuse transmittance',
]

if test_options[0].lower() in test:
if dni is None:
raise RuntimeError(
"Need to set 'dni' keyword to perform 'Clearness index' test in normalized_rradiance_test()."
)
if dhi is None:
raise RuntimeError(
"Need to set 'dhi' keyword to perform 'Clearness index' test in normalized_rradiance_test()."
)

if use_dask and isinstance(self._ds[dni].data, da.Array):
sza = da.array(sza)
Sa = da.array(Sa)
K_direct = self._ds[dni].values / Sa
K_diffuse = self._ds[dhi].values / (Sa * np.cos(np.radians(sza)))
K_total = K_direct + K_diffuse
index = (K_direct > K_total).compute()

else:
K_direct = self._ds[dni].values / Sa
K_diffuse = self._ds[dhi].values / (Sa * np.cos(np.radians(sza)))
K_total = K_direct + K_diffuse
index = K_direct > K_total

test_meaning = "Normalized direct normal irradiance greater than total transmittance."
self._ds.qcfilter.add_test(
dni,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)
self._ds.qcfilter.add_test(
dhi,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)

if test_options[1].lower() in test:
if dni is None:
raise RuntimeError(
"Need to set 'dni' keyword to perform 'Upper total transmittance' test in normalized_rradiance_test()."
)
if dhi is None:
raise RuntimeError(
"Need to set 'dhi' keyword to perform 'Upper total transmittance' test in normalized_rradiance_test()."
)

if use_dask and isinstance(self._ds[dni].data, da.Array):
sza = da.array(sza)
Sa = da.array(Sa)
K_direct = self._ds[dni].values / Sa
K_diffuse = self._ds[dhi].values / (Sa * np.cos(np.radians(sza)))
K_total = K_direct + K_diffuse
index = (K_total > upper_total_transmittance_limit) & (sza > 0) & (sza < 90)
index = index.compute()

else:
K_direct = self._ds[dni].values / Sa
K_diffuse = self._ds[dhi].values / (Sa * np.cos(np.radians(sza)))
K_total = K_direct + K_diffuse
index = (K_total > upper_total_transmittance_limit) & (sza > 0) & (sza < 90)

test_meaning = f"Total transmittance greater than {upper_total_transmittance_limit}"
self._ds.qcfilter.add_test(
dni,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)
self._ds.qcfilter.add_test(
dhi,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)

if test_options[2].lower() in test:
if dni is None:
raise RuntimeError(
"Need to set 'dni' keyword to perform 'Upper direct transmittance' test in normalized_rradiance_test()."
)
if ghi is None:
raise RuntimeError(
"Need to set 'ghi' keyword to perform 'Upper direct transmittance' test in normalized_rradiance_test()."
)

if use_dask and isinstance(self._ds[dni].data, da.Array):
Sa = da.array(Sa)
K_direct = self._ds[dni].values / Sa
alt = da.array(
convert_units(self._ds[alt_name].values, self._ds[alt_name].attrs['units'], 'm')
)
index1 = (K_direct > 0) & (self._ds[ghi].values > 50.0)
index2 = K_direct > ((1100 + 0.03 * alt) / Sa)
index = (index1 & index2).compute()

else:
K_direct = self._ds[dni].values / Sa
alt = convert_units(
self._ds[alt_name].values, self._ds[alt_name].attrs['units'], 'm'
)
index1 = (K_direct > 0) & (self._ds[ghi].values > 50.0)
index2 = K_direct > ((1100 + 0.03 * alt) / Sa)
index = index1 & index2

test_meaning = "Direct transmittance greater than upper direct transmittance limit"
self._ds.qcfilter.add_test(
dni,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)

if test_options[3].lower() in test:
if dhi is None:
raise RuntimeError(
"Need to set 'dhi' keyword to perform 'Upper diffuse transmittance' test in normalized_rradiance_test()."
)

if use_dask and isinstance(self._ds[dhi].data, da.Array):
sza = da.array(sza)
Sa = da.array(Sa)
K_diffuse = self._ds[dhi].values / (Sa * np.cos(np.radians(sza)))
index = (sza > 0) & (sza < 90) & (K_diffuse > upper_diffuse_transmittance_limit)
index = index.compute()

else:
K_diffuse = self._ds[dhi].values / (Sa * np.cos(np.radians(sza)))
index = (sza > 0) & (sza < 90) & (K_diffuse > upper_diffuse_transmittance_limit)

test_meaning = f"Diffuse transmittance greater than {upper_diffuse_transmittance_limit}"
self._ds.qcfilter.add_test(
dhi,
index=index,
test_assessment=test_assessment,
test_meaning=test_meaning,
)
Loading

0 comments on commit 6adb9be

Please sign in to comment.