detector_temp_flag.py

from __future__ import print_function
from scipy.signal import periodogram
import pandas as pd
import numpy as np
from aeri_tools.io.dmv import housekeeping
from matplotlib import pyplot as plt
import os
WIN_SIZE = 40
DETECTOR_TEMP_KEY = 'detectorTemp'
TIME_WINDOW = (60,600) #seconds
MIN_AMP = .02
WIN_TYPE = 'blackmanharris'
def sig_finder_function(rolling_window,**kwargs):
    """Function applied to detectorTemp dataset with pd.rolling_apply
    Calculates the dominant period and
    """

    #calculate the time difference and average frequency
    nsamples = rolling_window.size
    times = pd.Series(rolling_window.index)
    win_type = kwargs.get('win_type',WIN_TYPE)
    try:
        #get the time span covered by the window
        dt = (times.iloc[-1] - times.iloc[0]).seconds
    except AttributeError:
        #catch NaT indicies
        return np.nan
    else:
        #compute the dominant frequency/period
        fs = nsamples/dt
        f, Pxx_den = periodogram(rolling_window,fs,window=win_type)
        max_i = Pxx_den.argmax()
        if f[max_i] != 0:
            primary_T = 1./f[max_i]
        else:
            primary_T = float("inf")
        return primary_T
def series_rolling_apply(series,window,func,**kwargs):
    """rolling applies function to series by passing subsection of series,
    rather than np array
    """
    #lambda to map the passed function to a rolling window
    map_f=lambda i:np.nan if i<window else func(series.iloc[i-window:i],**kwargs)
    return pd.Series(map(map_f,range(series.size)),index = series.index)

def flag_detector_temp_oscillations(path_or_dmv,window_size = WIN_SIZE,std=10,**kwargs):
    """Performs power spectral density analysis on detector temperature,
    and flags the data if a consistent period is found
    """
    #read in the housekeeping dataframe
    housekeeping_df = housekeeping.get_all_housekeeping(path_or_dmv)
    #get the detector temp variable
    detector_temp = housekeeping_df[DETECTOR_TEMP_KEY]
    #detrend by subtracting the day's mean
    detrend_dtemp = detector_temp - detector_temp.mean()
    #approximate periodicity in data using a rolling periodogram
    periods = series_rolling_apply(detrend_dtemp,window_size,sig_finder_function,**kwargs)
    #check for big spikes in the data (these throw off the periodogram)
    #currently disabled
    dtt_roller = detrend_dtemp.rolling(window_size)
    amps = False & ((dtt_roller.max()-dtt_roller.min()) > .2)
    #look for low standard deviation in the period (= a consistent periodicity)
    stds = series_rolling_apply(periods,window_size,lambda x:x.std() < std)
    #flag as bad if there's a consistent periodicity or large jump in the data
    flags = stds | amps
    return flags, periods

def flatten_out(original_flags, pds):
    '''incrementing windows to help flatten out the 0/1 value as either
    0 or 1; removes small blips
    '''

    flags = original_flags.copy()
    values_dict = flags.value_counts()

    num_half_hour = int(len(flags)/48)
    curr_flags = []
    std_data = {}
    std_data['curr_std'] = []
    std_data['mean_std'] = []
    std_data['mean'] = []
    std_data['std_curr_std'] = []
    std_data['std_mean_std'] = []

    lower = 0
    upper = 1

    while upper <= len(flags):
        values_dict = original_flags[lower:upper].value_counts()
        if True in values_dict.keys():
            mean = values_dict[True] / len(flags[lower:upper])
        else:
            mean = 1 - (values_dict[False] / len(flags[lower:upper]))
        curr_std = np.std(pds[lower:upper])
        mean_std = np.mean(pds[lower:upper])
        std_data['mean'].append(mean)
        std_data['curr_std'].append(curr_std)
        std_data['mean_std'].append(mean_std)
        std_curr_std = np.std(std_data['curr_std'][lower:upper])
        std_data['std_curr_std'].append(std_curr_std)
        std_mean_std = np.std(std_data['mean_std'][lower:upper])
        std_data['std_mean_std'].append(std_mean_std)
        if std_mean_std < 20 or np.isnan(std_mean_std):
            flags[lower:upper] = flags[lower:upper].replace(False, True)
        else:
            flags[lower:upper] = flags[lower:upper].replace(True, False)

        upper += 1
        if upper - lower > num_half_hour:
            lower += 1

    new_flags = flags & (np.asarray(std_data['mean']) > 0.35)
    # new_flags = np.asarray(std_data['mean']) > 0.8
    # fig, ax = plt.subplots(2, figsize=(15,10),sharex=True)
    # plt.plot(flags)
    # plt.plot(new_flags)
    # plt.show()
    # plt.clf()
    return new_flags, std_data

def main(path):
    """Test function that displays detectorTemp alongside the computed flag
    using different windowing functions
    """
    dmv = housekeeping.get_all_housekeeping(path)
    detector_temp = dmv[DETECTOR_TEMP_KEY]
    fig,ax = plt.subplots(5,figsize=(15,10),sharex=True)
    plt.sca(ax[0])
    plt.plot(dmv['Time'],detector_temp,'r')
    # plt.ylim(76.9,77.6)
    plt.ylabel("Detector Temperature (K)")
    plt.grid()

    flags,pds = flag_detector_temp_oscillations(path,win_type=WIN_TYPE)
    flags, std_data = flatten_out(flags.copy(), pds)

    plt.sca(ax[1])
    plt.plot(dmv['Time'],pds,lw=2,label=WIN_TYPE)
    plt.grid(True)
    plt.ylabel("Dominant Period (s)")
    #plt.legend()

    plt.sca(ax[2])
    plt.plot(dmv['Time'],flags,lw=2,label=WIN_TYPE)
    plt.grid(True)
    plt.ylabel("Standard Deviation Flag")
    plt.xlabel("Time (Hours UTC)")
    plt.ylim(-0.1,1.1)
    plt.sca(ax[3])
    plt.plot(dmv['Time'], np.asarray(std_data['mean']))

    title = os.path.split(path)[-1][:6]
    plt.suptitle(title,fontsize=14)

    plt.sca(ax[4])
    plt.plot(dmv['Time'], std_data['curr_std'])     #blue
    plt.plot(dmv['Time'], std_data['mean_std'])     #green
    plt.plot(dmv['Time'], std_data['std_curr_std']) #red
    plt.plot(dmv['Time'], std_data['std_mean_std'])

    print('title: ', title)
    # plt.savefig("sgp_detect_pngs/detector_temp_flag_test_%s.png"%title,bbox_inches='tight')
    # plt.savefig("detector_temp_flag_test_%s_2.png"%title,bbox_inches='tight')
    # plt.show()
    plt.clf()
    plt.close()

if __name__ == '__main__':
    from sys import argv
    import os
    import time

    start_time = time.time()

    pathname = argv[1]

    skip_num = 0
    curr_num = 0
    print('skip_num is ', skip_num)

    if os.path.isdir(pathname):
        for filename in os.listdir(pathname):
            full_path = pathname + '/' + filename
            if os.path.isdir(full_path):
                for fname in os.listdir(full_path):
                    fuller_path = full_path + '/' + fname
                    if os.path.isfile(fuller_path) and 'B1.CXS' in fname:
                        if curr_num >= skip_num:
                            main(fuller_path)
                        curr_num += 1
                        print(curr_num, ': ', fuller_path)
            elif os.path.isfile(full_path) and 'B1.CXS' in filename:
                if curr_num >= skip_num:
                    main(full_path)
                curr_num += 1
                print(curr_num, ': ', filename)
    else:
        main(pathname)

    print('execution time: %d minute(s), %.2f second(s)' % ((time.time()-start_time)//60, (time.time()-start_time)%60))