Source code for act.qc.radiometer_tests

"""
Tests specific to radiometers.

"""

import warnings

import dask
import numpy as np
import pandas as pd
import xarray as xr
from scipy.fftpack import rfft, rfftfreq

from act.utils.datetime_utils import determine_time_delta
from act.utils.geo_utils import is_sun_visible


[docs]def fft_shading_test( ds, variable='diffuse_hemisp_narrowband_filter4', fft_window=30, shad_freq_lower=[0.008, 0.017], shad_freq_upper=[0.0105, 0.0195], ratio_thresh=[3.15, 1.2], time_interval=None, smooth_window=5, shading_thresh=0.4, ): """ Function to test shadowband radiometer (MFRSR, RSS, etc) instruments for shading related problems. Program was adapted by Adam Theisen from the method defined in Alexandrov et al 2007 to process on a point by point basis using a window of data around that point for the FFT analysis. For ARM data, testing has found that this works the best on narrowband filter4 for MFRSR data. Function has been tested and is in use by the ARM DQ Office for problem detection. It is know to have some false positives at times. Need to run ds.clean.cleanup() ahead of time to ensure proper addition to QC variable Parameters ---------- ds : xarray.Dataset Xarray dataset variable : string Name of variable to process fft_window : int Number of samples to use in the FFT window. Default is +- 30 samples Note: this is +- so the full window will be double shad_freq_lower : list Lower frequency over which to look for peaks in FFT shad_freq_upper : list Upper frequency over which to look for peaks in FFT ratio_thresh : list Threshold for each freq window to flag data. I.e. if the peak is 3.15 times greater than the surrounding area time_interval : float Sampling rate of the instrument smooth_window : int Number of samples to use in smoothing FFTs before analysis shading_thresh : float After smoothing, the value over which is considered a shading signal Returns ------- ds : xarray.Dataset Xarray dataset tested for shading problems References ---------- Alexandrov, Mikhail & Kiedron, Peter & Michalsky, Joseph & Hodges, Gary & Flynn, Connor & Lacis, Andrew. (2007). Optical depth measurements by shadow-band radiometers and their uncertainties. Applied optics. 46. 8027-38. 10.1364/AO.46.008027. """ # Get time and data from variable time = ds['time'].values data = ds[variable].values if 'missing_value' in ds[variable].attrs: missing = ds[variable].attrs['missing_value'] else: missing = -9999.0 # Get time interval between measurements if time_interval is None: dt = determine_time_delta(time) else: dt = time_interval # Compute the FFT for each point +- window samples task = [] sun_up = is_sun_visible(latitude=ds['lat'].values, longitude=ds['lon'].values, date_time=time) for t in range(len(time)): sind = t - fft_window eind = t + fft_window if sind < 0: sind = 0 if eind > len(time): eind = len(time) # Get data and remove all nan/missing values d = data[sind:eind] idx = (d != missing) & (np.isnan(d) is not True) index = np.where(idx) d = d[index] # Add to task for dask processing task.append( dask.delayed(fft_shading_test_process)( time[t], d, shad_freq_lower=shad_freq_lower, shad_freq_upper=shad_freq_upper, ratio_thresh=ratio_thresh, time_interval=dt, is_sunny=sun_up[t], ) ) # Process using dask result = dask.compute(*task) # Run data through a rolling median to filter out singular # false positives shading = [r['shading'] for r in result] shading = pd.Series(shading).rolling(window=smooth_window, min_periods=1).median() # Find indices where shading is indicated idx = np.asarray(shading) > shading_thresh index = np.where(idx) # Add test to QC Variable desc = 'FFT Shading Test' ds.qcfilter.add_test(variable, index=index, test_meaning=desc) # Prepare frequency and fft variables for adding to the dataset fft = np.empty([len(time), fft_window * 2]) fft[:] = np.nan freq = np.empty([len(time), fft_window * 2]) freq[:] = np.nan for i, r in enumerate(result): dummy = r['fft'] fft[i, 0 : len(dummy)] = dummy dummy = r['freq'] freq[i, 0 : len(dummy)] = dummy attrs = { 'units': '', 'long_name': 'FFT Results for Shading Test', 'upper_freq': shad_freq_upper, 'lower_freq': shad_freq_lower, } fft_window = xr.DataArray( range(fft_window * 2), dims=['fft_window'], attrs={'long_name': 'FFT Window', 'units': '1'}, ) da = xr.DataArray( fft, dims=['time', 'fft_window'], attrs=attrs, coords=[ds['time'], fft_window] ) ds['fft'] = da attrs = {'units': '', 'long_name': 'FFT Frequency Values for Shading Test'} da = xr.DataArray( freq, dims=['time', 'fft_window'], attrs=attrs, coords=[ds['time'], fft_window] ) ds['fft_freq'] = da return ds
[docs]def fft_shading_test_process( time, data, shad_freq_lower=None, shad_freq_upper=None, ratio_thresh=None, time_interval=None, is_sunny=None, ): """ Processing function to do the FFT calculations/thresholding Parameters ---------- time : datetime Center time of calculation used for calculating sunrise/sunset data : list Data for run through fft processing shad_freq_lower : list Lower limits of freqencies to look for shading issues shad_freq_upper : list Upper limits of freqencies to look for shading issues ratio_thresh : list Thresholds to apply, corresponding to frequencies chosen time_interval : float Time interval of data Returns ------- shading : int Binary to indicate shading problem (1) or not (0) """ if not is_sunny: return {'shading': 0, 'fft': [np.nan] * len(data), 'freq': [np.nan] * len(data)} # FFT Algorithm fftv = abs(rfft(data)) freq = rfftfreq(fftv.size, d=time_interval) # Get FFT data under threshold with warnings.catch_warnings(): warnings.filterwarnings('ignore', category=RuntimeWarning) idx = fftv > 1.0 index = np.where(idx) fftv[index] = np.nan freq[index] = np.nan # Return if FFT is empty if len(fftv) == 0: return {'shading': 0, 'fft': [np.nan] * len(data), 'freq': [np.nan] * len(data)} # Commented out as it seems to work better without smoothing # fftv=pd.DataFrame(data=fftv).rolling(min_periods=3,window=3,center=True).mean().values.flatten() ratio = [] # Calculates the ratio (size) of the peaks in the FFT to the surrounding # data wind = 3 # Run through each frequency to look for peaks # Calculate threshold of peak value to surrounding values for i in range(len(shad_freq_lower)): with warnings.catch_warnings(): warnings.filterwarnings('ignore', category=RuntimeWarning) idx = np.logical_and(freq > shad_freq_lower[i], freq < shad_freq_upper[i]) index = np.where(idx) if len(index[0]) == 0: continue peak = max(fftv[index]) index = index[0] sind = index[0] - wind if sind < 0: sind = 0 eind = index[-1] + wind if eind > len(fftv): eind = len(fftv) if len(range(sind, index[0])) == 0 or len(range(index[-1], eind)) == 0: ratio.append(0.0) else: # Calculates to the left/right of each peak peak_l = max(fftv[range(sind, index[0])]) peak_r = max(fftv[range(index[-1], eind)]) mean_value = np.mean([peak_l, peak_r]) if mean_value == 0.0: ratio.append(np.nan) else: ratio.append(peak / mean_value) # Checks ratios against thresholds for each freq range shading = 0 if len(ratio) > 0: pass1 = False pass2 = False if ratio[0] > ratio_thresh[0]: pass1 = True if len(ratio) > 1: if ratio[1] > ratio_thresh[1]: pass2 = True else: pass2 = True if pass1 and pass2: shading = 1 return {'shading': shading, 'fft': fftv, 'freq': freq}