From 6adb9be009d3f98830b62211cf411ae121360897 Mon Sep 17 00:00:00 2001 From: Ken Kehoe Date: Thu, 7 Nov 2024 09:14:49 -0700 Subject: [PATCH] Adding radiation quality control tests (#874) * 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. --- act/qc/bsrn_tests.py | 310 ++++++++++++++++++++++++++++- act/retrievals/radiation.py | 245 +++++++++++++++++++++-- examples/qc/plot_qc_bsrn.py | 13 ++ tests/qc/test_bsrn_tests.py | 120 ++++++++++- tests/retrievals/test_radiation.py | 56 ++++++ 5 files changed, 710 insertions(+), 34 deletions(-) diff --git a/act/qc/bsrn_tests.py b/act/qc/bsrn_tests.py index 1e3c5e4ef4..d8ff0de574 100644 --- a/act/qc/bsrn_tests.py +++ b/act/qc/bsrn_tests.py @@ -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 @@ -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 @@ -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 @@ -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 -------- @@ -470,6 +472,8 @@ 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', @@ -477,13 +481,14 @@ def bsrn_comparison_tests( '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 @@ -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 ' @@ -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 @@ -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 ' @@ -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 ' @@ -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 ' @@ -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, + ) diff --git a/act/retrievals/radiation.py b/act/retrievals/radiation.py index 09c5d1ade3..83913b9a95 100644 --- a/act/retrievals/radiation.py +++ b/act/retrievals/radiation.py @@ -8,6 +8,7 @@ from scipy.constants import Stefan_Boltzmann from act.utils.geo_utils import get_solar_azimuth_elevation +from act.utils.data_utils import convert_units def calculate_dsh_from_dsdh_sdn( @@ -19,45 +20,247 @@ def calculate_dsh_from_dsdh_sdn( ): """ - Function to derive the downwelling shortwave hemispheric irradiance from the + Function to derive the downwelling shortwave hemispheric irradiance (dsh) from the downwelling shortwave diffuse hemispheric irradiance (dsdh) and the shortwave - direct normal irradiance (sdn) at a given location (lat,lon) + direct normal irradiance (sdn). The derived values are added the returned + Datasets as a new varible. Parameters ---------- ds : xarray.Dataset - Xarray dataset where variables for these calculations are stored + Xarray dataset where variables for these calculations are stored. dsdh : str Name of the downwelling shortwave diffuse hemispheric irradiance field to use. - Defaults to downwelling_sw_diffuse_hemisp_irradiance. sdn : str Name of shortwave direct normal irradiance field to use. - Defaults to shortwave_direct_normal_irradiance. lat : str - Name of latitude field in dataset to use. Defaults to 'lat'. + Name of latitude variable in dataset to use for deriving solar zenight angle. + lon : str + Name of longitude variable in dataset to use for deriving solar zenight angle. + + Returns + ------- + + ds: xarray.Dataset + ACT Xarray Dataset with calculations included as new variable. + + """ + + ds = calculate_ghi_from_dni_dhi(ds, dni=sdn, dhi=dsdh, lat=lat, lon=lon) + + return ds + + +def calculate_ghi_from_dni_dhi( + ds, + dni='short_direct_normal', + dhi='down_short_diffuse_hemisp', + derived='derived_down_short_hemisp', + long_name='Derived global horizontal irradiance', + solar_zenith=None, + lat='lat', + lon='lon', +): + """ + + Function to derive the Global Horizontal Irradiance (GHI) from + Direct Normal Irradiance (DNI) and Diffuse Horizontal + Irradiance (DHI). The derived values are added the returned Datasets as a new DataArray. + + Parameters + ---------- + ds : xarray.Dataset + Xarray dataset where variables used for these calculations are stored + dni : str + Name of the direct normal irradiance DataArray to use. + dhi : str + Name of diffuse hemispheric irradiance DataArray to use. + derived : str + Name of new diffuse horizontal irradiance DataArray. + long_name : str + Long name used in new DataArray. + solar_zenith : str or None + Name of solar zenith DataArray in Dataset. If set to None will calculate using + location and time variables from Dataset. + lat : str + Name of latitude DataArray in dataset to use for deriving solar zenight angle. + Ignored if solar_zenith provided. lon : str - Name of longitued field in dataset to use. Defaults to 'lon'. + Name of longitued DataArray in dataset to use for deriving solar zenight angle. + Ignored if solar_zenith provided. Returns ------- ds: xarray.Dataset - ACT Xarray Dataset with calculations included as new variables. + ACT Xarray Dataset with global horizontal irradiance included as new DataArray. """ - # Calculating Derived Down Short Hemisp - elevation, _, _ = get_solar_azimuth_elevation(ds[lat].values, ds[lon].values, ds['time'].values) - solar_zenith = np.cos(np.radians(90.0 - elevation)) - dsh = ds[dsdh].values + (solar_zenith * ds[sdn].values) + # Get solar zenith angle + if solar_zenith is not None: + sz = ds[solar_zenith].values + sz = convert_units(sz, ds[solar_zenith].attrs['units'], 'radians') + cos_sz = np.cos(sz) + else: + elevation, _, _ = get_solar_azimuth_elevation( + ds[lat].values, ds[lon].values, ds['time'].values + ) + cos_sz = np.cos(np.radians(90.0 - elevation)) + + ghi = ds[dhi].values + (cos_sz * ds[dni].values) - # Add data back to DataArray - ds['derived_down_short_hemisp'] = xr.DataArray( - dsh, + # Add data into Dataset + ds[derived] = xr.DataArray( + ghi, dims=['time'], attrs={ - 'long_name': 'Derived Downwelling Shortwave Hemispheric Irradiance', - 'units': 'W/m^2', + 'long_name': long_name, + 'units': ds[dhi].attrs['units'], + }, + ) + + return ds + + +def calculate_dni_from_dhi_ghi( + ds, + dhi='down_short_diffuse_hemisp', + ghi='down_short_hemisp', + derived='derived_short_direct_normal', + long_name='Derived direct normal irradiance', + solar_zenith=None, + lat='lat', + lon='lon', +): + """ + + Function to derive the Direct Normal Irradiance (DNI) from + Diffuse Horizontal Irradiance (DHI) and Global Horizontal + Irradiance (GHI). The derived values are added the returned Datasets as a new DataArray. + + Parameters + ---------- + ds : xarray.Dataset + Xarray dataset where variables used for these calculations are stored + dhi : str + Name of the diffuse horizontal irradiance DataArray to use. + ghi : str + Name of global hemispheric irradiance DataArray to use. + derived : str + Name of new direct normal irradiance DataArray. + long_name : str + Long name used in new DataArray. + solar_zenith : str or None + Name of solar zenith DataArray in Dataset. If set to None will calculate using + location and time variables from Dataset. + lat : str + Name of latitude DataArray in dataset to use for deriving solar zenight angle. + Ignored if solar_zenith provided. + lon : str + Name of longitued DataArray in dataset to use for deriving solar zenight angle. + Ignored if solar_zenith provided. + + Returns + ------- + + ds: xarray.Dataset + ACT Xarray Dataset with direct normal irradiance included as new DataArray. + + """ + + # Get solar zenith angle + if solar_zenith is not None: + sz = ds[solar_zenith].values + sz = convert_units(sz, ds[solar_zenith].attrs['units'], 'radians') + cos_sz = np.cos(sz) + else: + elevation, _, _ = get_solar_azimuth_elevation( + ds[lat].values, ds[lon].values, ds['time'].values + ) + cos_sz = np.cos(np.radians(90.0 - elevation)) + + dni = (ds[ghi].values - ds[dhi].values) / cos_sz + + # Add data into Dataset + ds[derived] = xr.DataArray( + dni, + dims=['time'], + attrs={ + 'long_name': long_name, + 'units': ds[dhi].attrs['units'], + }, + ) + + return ds + + +def calculate_dhi_from_dni_ghi( + ds, + dni='short_direct_normal', + ghi='down_short_hemisp', + derived='derived_down_short_diffuse_hemisp', + long_name='Derived diffuse horizontal irradiance', + solar_zenith=None, + lat='lat', + lon='lon', +): + """ + + Function to derive the Diffuse Horizontal Irradiance (DHI) from + Direct Normal Irradiance (DNI) and Global Horizontal + Irradiance (GHI). The derived values are added the returned Datasets as a new DataArray. + + Parameters + ---------- + ds : xarray.Dataset + Xarray dataset where variables used for these calculations are stored + dni : str + Name of the dirret normal irradiance DataArray to use. + ghi : str + Name of global hemispheric irradiance DataArray to use. + derived : str + Name of new diffuse horizontal irradiance DataArray. + long_name : str + Long name used in new DataArray. + solar_zenith : str or None + Name of solar zenith DataArray in Dataset. If set to None will calculate using + location and time variables from Dataset. + lat : str + Name of latitude DataArray in dataset to use for deriving solar zenight angle. + Ignored if solar_zenith provided. + lon : str + Name of longitued DataArray in dataset to use for deriving solar zenight angle. + Ignored if solar_zenith provided. + + Returns + ------- + + ds: xarray.Dataset + ACT Xarray Dataset with diffuse horizontal irradiance included as new DataArray. + + """ + + # Get solar zenith angle + if solar_zenith is not None: + sz = ds[solar_zenith].values + sz = convert_units(sz, ds[solar_zenith].attrs['units'], 'radians') + cos_sz = np.cos(sz) + else: + elevation, _, _ = get_solar_azimuth_elevation( + ds[lat].values, ds[lon].values, ds['time'].values + ) + cos_sz = np.cos(np.radians(90.0 - elevation)) + + dhi = ds[ghi].values - (ds[dni].values * cos_sz) + + # Add data into Dataset + ds[derived] = xr.DataArray( + dhi, + dims=['time'], + attrs={ + 'long_name': long_name, + 'units': ds[ghi].attrs['units'], }, ) @@ -86,10 +289,10 @@ def calculate_irradiance_stats( Name of the second irradiance variable diff_output_variable : str Variable name to store the difference results - Defaults to 'diff[underscore]'+variable + Defaults to 'diff_' + variable ratio_output_variable : str Variable name to store the ratio results - Defaults to 'ratio[underscore]'+variable + Defaults to 'ratio_' + variable Returns ------- @@ -101,8 +304,10 @@ def calculate_irradiance_stats( if variable is None or variable2 is None: return ds + if diff_output_variable is None: diff_output_variable = 'diff_' + variable + if ratio_output_variable is None: ratio_output_variable = 'ratio_' + variable @@ -145,7 +350,7 @@ def calculate_net_radiation( ): """ - Function to calculate the net radiation from upwelling short and long-wave irradiance and + Function to calculate the net radiation from upwelling short and long-wave irradiance and downwelling short and long-wave hemisperic irradiances Parameters diff --git a/examples/qc/plot_qc_bsrn.py b/examples/qc/plot_qc_bsrn.py index 37e3d907e7..6ec1afe15d 100644 --- a/examples/qc/plot_qc_bsrn.py +++ b/examples/qc/plot_qc_bsrn.py @@ -64,6 +64,19 @@ direct_normal_SW_dn_name='short_direct_normal', ) +# Add K-tests QC to ancillary QC variables. +ds.qcfilter.normalized_rradiance_test( + [ + 'Clearness index', + 'Upper total transmittance', + 'Upper direct transmittance', + 'Upper diffuse transmittance', + ], + dni='short_direct_normal', + dhi='down_short_hemisp', + ghi='down_short_diffuse_hemisp', +) + # Creat Plot Display and plot data including embedded QC from data file variable = 'down_short_hemisp' display = act.plotting.TimeSeriesDisplay(ds, figsize=(15, 10), subplot_shape=(2,)) diff --git a/tests/qc/test_bsrn_tests.py b/tests/qc/test_bsrn_tests.py index c83d1c65db..81e0c4b430 100644 --- a/tests/qc/test_bsrn_tests.py +++ b/tests/qc/test_bsrn_tests.py @@ -1,5 +1,5 @@ import copy - +import pytest import dask.array as da import numpy as np import xarray as xr @@ -232,6 +232,11 @@ def test_bsrn_limits_test(): ) ds['up_long_hemisp'].values[0:100] = ds['up_long_hemisp'].values[0:100] - 200 + # Closure test modified data + ds['down_short_diffuse_hemisp'].values[714:814] = ( + ds['down_short_diffuse_hemisp'].values[714:814] + 50.0 + ) + ds.qcfilter.bsrn_comparison_tests( [ 'Global over Sum SW Ratio', @@ -240,6 +245,7 @@ def test_bsrn_limits_test(): 'LW down to air temp', 'LW up to air temp', 'LW down to LW up', + 'Closure', ], gbl_SW_dn_name='down_short_hemisp', glb_diffuse_SW_dn_name='down_short_diffuse_hemisp', @@ -256,11 +262,11 @@ def test_bsrn_limits_test(): # Ratio of Global over Sum SW result = ds.qcfilter.get_qc_test_mask('down_short_hemisp', test_number=5) - assert np.sum(result) == 190 + assert np.sum(result) == 276 # Diffuse Ratio result = ds.qcfilter.get_qc_test_mask('down_short_hemisp', test_number=6) - assert np.sum(result) == 47 + assert np.sum(result) == 103 # Shortwave up comparison result = ds.qcfilter.get_qc_test_mask('up_short_hemisp', test_number=5) @@ -274,6 +280,112 @@ def test_bsrn_limits_test(): result = ds.qcfilter.get_qc_test_mask('down_long_hemisp_shaded', test_number=5) assert np.sum(result) == 976 - # Lonwave down to longwave up comparison + # Longwave down to longwave up comparison result = ds.qcfilter.get_qc_test_mask('down_long_hemisp_shaded', test_number=6) assert np.sum(result) == 100 + + # Closure test + assert ( + ds['qc_down_short_hemisp'].attrs['flag_meanings'][6] + == 'Closure test indicating value outside of expected range' + ) + result = ds.qcfilter.get_qc_test_mask('down_short_hemisp', test_number=7) + assert np.sum(result) == 38 + assert ( + ds['qc_short_direct_normal'].attrs['flag_meanings'][6] + == 'Closure test indicating value outside of expected range' + ) + result = ds.qcfilter.get_qc_test_mask('short_direct_normal', test_number=7) + assert np.sum(result) == 38 + assert ( + ds['qc_down_short_diffuse_hemisp'].attrs['flag_meanings'][7] + == 'Closure test indicating value outside of expected range' + ) + result = ds.qcfilter.get_qc_test_mask('down_short_diffuse_hemisp', test_number=8) + assert np.sum(result) == 38 + + +def test_normalized_rradiance_test(): + keep_vars = [ + 'short_direct_normal', + 'down_short_diffuse_hemisp', + 'down_short_hemisp', + 'lat', + 'lon', + 'alt', + ] + ds = read_arm_netcdf(EXAMPLE_BRS, keep_variables=keep_vars) + tests = [ + 'Clearness index', + 'Upper total transmittance', + 'Upper direct transmittance', + 'Upper diffuse transmittance', + ] + for test in tests: + with pytest.raises(RuntimeError): + ds.qcfilter.normalized_rradiance_test(test) + + for use_dask in [False, True]: + ds = read_arm_netcdf(EXAMPLE_BRS, keep_variables=keep_vars) + data = ds['short_direct_normal'].values + data[1050:1100] = data[1050:1100] + 400 + ds['short_direct_normal'].values = data + ds.qcfilter.normalized_rradiance_test( + tests, + dni='short_direct_normal', + dhi='down_short_diffuse_hemisp', + ghi='down_short_hemisp', + use_dask=use_dask, + upper_total_transmittance_limit=1.4, + upper_diffuse_transmittance_limit=0.6, + ) + + test_number = ( + ds['qc_down_short_diffuse_hemisp'] + .attrs['flag_meanings'] + .index('Normalized direct normal irradiance greater than total transmittance.') + + 1 + ) + result = ds.qcfilter.get_qc_test_mask('down_short_diffuse_hemisp', test_number=test_number) + assert np.sum(np.where(result)) == 15780 + test_number = ( + ds['qc_down_short_diffuse_hemisp'] + .attrs['flag_meanings'] + .index('Total transmittance greater than 1.4') + + 1 + ) + result = ds.qcfilter.get_qc_test_mask('down_short_diffuse_hemisp', test_number=test_number) + assert np.sum(np.where(result)) == 789 + test_number = ( + ds['qc_down_short_diffuse_hemisp'] + .attrs['flag_meanings'] + .index('Diffuse transmittance greater than 0.6') + + 1 + ) + result = ds.qcfilter.get_qc_test_mask('down_short_diffuse_hemisp', test_number=test_number) + assert np.sum(np.where(result)) == 2367 + + test_number = ( + ds['qc_short_direct_normal'] + .attrs['flag_meanings'] + .index('Normalized direct normal irradiance greater than total transmittance.') + + 1 + ) + result = ds.qcfilter.get_qc_test_mask('short_direct_normal', test_number=test_number) + assert np.sum(np.where(result)) == 15780 + test_number = ( + ds['qc_short_direct_normal'] + .attrs['flag_meanings'] + .index('Total transmittance greater than 1.4') + + 1 + ) + result = ds.qcfilter.get_qc_test_mask('short_direct_normal', test_number=test_number) + assert np.sum(np.where(result)) == 789 + test_number = ( + ds['qc_short_direct_normal'] + .attrs['flag_meanings'] + .index('Direct transmittance greater than upper direct transmittance limit') + + 1 + ) + result = ds.qcfilter.get_qc_test_mask('short_direct_normal', test_number=test_number) + assert np.sum(np.where(result)) == 18547 diff --git a/tests/retrievals/test_radiation.py b/tests/retrievals/test_radiation.py index fa8357efd2..40a1b13b35 100644 --- a/tests/retrievals/test_radiation.py +++ b/tests/retrievals/test_radiation.py @@ -7,10 +7,66 @@ def test_calculate_sirs_variable(): sirs_ds = act.io.arm.read_arm_netcdf(act.tests.sample_files.EXAMPLE_SIRS) met_ds = act.io.arm.read_arm_netcdf(act.tests.sample_files.EXAMPLE_MET1) + mfrsr_ds = act.io.arm.read_arm_netcdf(act.tests.sample_files.EXAMPLE_MFRSR) + elevation, _, _ = act.utils.geo_utils.get_solar_azimuth_elevation( + sirs_ds['lat'].values, sirs_ds['lon'].values, sirs_ds['time'].values + ) + sirs_ds['solar_zenith_angle'] = xr.DataArray( + 90.0 - elevation, + dims=['time'], + attrs={ + 'long_name': 'Solar zenith angle', + 'units': 'degrees', + }, + ) + + # old ds = act.retrievals.radiation.calculate_dsh_from_dsdh_sdn(sirs_ds) assert np.isclose(np.nansum(ds['derived_down_short_hemisp'].values), 61157.71, atol=0.1) + # new + ds = act.retrievals.radiation.calculate_ghi_from_dni_dhi(sirs_ds) + assert np.isclose(np.nansum(ds['derived_down_short_hemisp'].values), 61157.71, atol=0.1) + + # Test with providing solar_zenith variable in Dataset + ds = act.retrievals.radiation.calculate_ghi_from_dni_dhi( + sirs_ds, solar_zenith='solar_zenith_angle' + ) + assert np.isclose(np.nansum(ds['derived_down_short_hemisp'].values), 61157.71, atol=0.1) + + ds = act.retrievals.radiation.calculate_dhi_from_dni_ghi(sirs_ds) + assert np.isclose(np.nansum(ds['derived_down_short_diffuse_hemisp'].values), 59825.39, atol=0.1) + + # Test with providing solar_zenith variable in Dataset + ds = act.retrievals.radiation.calculate_dhi_from_dni_ghi( + sirs_ds, solar_zenith='solar_zenith_angle' + ) + assert np.isclose(np.nansum(ds['derived_down_short_diffuse_hemisp'].values), 59825.39, atol=0.1) + + # Test direct normal + mfrsr_ds = act.retrievals.radiation.calculate_dni_from_dhi_ghi( + mfrsr_ds, + ghi='hemisp_narrowband_filter3', + dhi='diffuse_hemisp_narrowband_filter3', + derived='derived_direct_normal_narrowband_filter3', + ) + assert np.isclose( + np.nansum(mfrsr_ds['derived_direct_normal_narrowband_filter3'].values), 2553.92, atol=0.1 + ) + + # Test with providing solar_zenith variable in Dataset + mfrsr_ds = act.retrievals.radiation.calculate_dni_from_dhi_ghi( + mfrsr_ds, + ghi='hemisp_narrowband_filter3', + dhi='diffuse_hemisp_narrowband_filter3', + derived='derived_direct_normal_narrowband_filter3', + solar_zenith='solar_zenith_angle', + ) + assert np.isclose( + np.nansum(mfrsr_ds['derived_direct_normal_narrowband_filter3'].values), 2553.59, atol=0.1 + ) + ds = act.retrievals.radiation.calculate_irradiance_stats( ds, variable='derived_down_short_hemisp',