opc_processing_code.py

#####################################################################################################
# Copyright 2022 by Bo Chen, Texas A&M University in collaboration with Sandia National Laboratories.
# All rights reserved.
#####################################################################################################

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import datetime

class OPC_Measurement:

    @classmethod
    def load_file(cls, file_dir, start_time = None, end_time = None, diagnose = False):
        # initialize object through file input

        lower_bin_boundary = np.loadtxt(file_dir, dtype='str', delimiter=',', skiprows=8, max_rows=1)[1:].astype(float) # last number is the size of the upper boundary of the last bin
        bin_mean = np.loadtxt(file_dir, dtype='str', delimiter=',', skiprows=9, max_rows=1)[1:].astype(float)
        opc_data = np.loadtxt(file_dir, delimiter=',', skiprows=14)

        return cls(opc_data, lower_bin_boundary, bin_mean, start_time = start_time, end_time = end_time, diagnose = diagnose)

    def __init__(self, opc_data, lower_bin_boundary, bin_mean, start_time = None, end_time = None, diagnose = False):

        self.lower_bin_boundary = lower_bin_boundary
        self.bin_mean = bin_mean

        self.OLE_TIME_ZERO = datetime.datetime(1899, 12, 30, 0, 0, 0)
        self.OA_datetime = opc_data[:, 0] # temporary, unclipped

        if start_time is not None:
            starttime_oa = start_time - self.OLE_TIME_ZERO
            starttime_oa = starttime_oa.total_seconds()/86400

            OA_minus_start = self.OA_datetime-starttime_oa
            self.starttime_num = np.abs(OA_minus_start).argmin()
        else:
            self.starttime_num = 0

        if end_time is not None:
            endtime_oa = end_time - self.OLE_TIME_ZERO
            endtime_oa = endtime_oa.total_seconds()/86400

            OA_minus_end = self.OA_datetime-endtime_oa
            self.endtime_num = np.abs(OA_minus_end).argmin()
        else:
            self.endtime_num = None
        
        # accquire Datetime
        self.opc_data_clipped = opc_data[self.starttime_num : self.endtime_num]
        self.OA_datetime =  self.OA_datetime[self.starttime_num : self.endtime_num] # clipped
        self.oa_difference = self.OA_datetime[1:] - self.OA_datetime[:-1]
    
        if diagnose:
            print("start time num is "+str(self.starttime_num))
            print("end time num is "+str(self.endtime_num))
            print("OA_datetime shape is: " + str(self.OA_datetime.shape))
            fig = plt.figure()
            plt.plot(self.oa_difference)
        
        def ole2datetime(oledt):
            return  self.OLE_TIME_ZERO + datetime.timedelta(days=float(oledt))
        vfunc_time = np.vectorize(ole2datetime)
        self.current_datetime = vfunc_time( self.OA_datetime)

        # accquire Data
        self.bins_counts_per_s = self.opc_data_clipped[:, 1:1+np.shape(self.bin_mean)[0]] #/s
        self.mean_of_tof1 = self.opc_data_clipped[:, 25] #us
        self.mean_of_tof3 = self.opc_data_clipped[:, 26] #us
        self.mean_of_tof5 = self.opc_data_clipped[:, 27] #us
        self.mean_of_tof7 = self.opc_data_clipped[:, 28] #us
        self.counts_per_s = self.opc_data_clipped[:, 29] #/s
        self.sample_flowrate = self.opc_data_clipped[:, 31] #ml/s
        self.temperature = self.opc_data_clipped[:, 32] #degree C
        self.relative_humidity = self.opc_data_clipped[:, 33] # %
        self.laser_status = self.opc_data_clipped[:, 39] # %
        self.pm_1 = self.opc_data_clipped[:, 40] # %
        self.pm_2_5 = self.opc_data_clipped[:, 41] # %
        self.pm_10 = self.opc_data_clipped[:, 42] # %
        self.rollmean_pm_1 = self.opc_data_clipped[:, 43] # %
        self.rollmean_pm_2_5 = self.opc_data_clipped[:, 44] # %
        self.rollmean_pm_10 = self.opc_data_clipped[:, 45] # %

        self.concentration = self.counts_per_s / self.sample_flowrate * 1000 #/L
        flowrate_stacked = np.vstack([self.sample_flowrate] * self.bins_counts_per_s.shape[1])
        flowrate_stacked = np.transpose(flowrate_stacked)
        self.bins_concentration = np.divide(self.bins_counts_per_s, flowrate_stacked) * 1000 #/L

        dlogDp = np.vstack([np.log10(self.lower_bin_boundary[1:]/self.lower_bin_boundary[:-1])] * self.bins_concentration.shape[0])

        self.bins_dNdlogDP = np.divide(self.bins_concentration, dlogDp)

    def timeseries_plot(self, fig, axes, data, unit, x_gap, moving_average_window = None, c = "k"):
        time = self.current_datetime

        x_nums = np.arange(0, np.shape(time)[0])
        if moving_average_window is not None:
            axes.plot(x_nums, data, c = "lightgray")
            moving_average = OPC_Measurement.movingaverage(data, moving_average_window)
            axes.plot(x_nums, moving_average, c = c)
            axes.legend(['measurement', 'moving average window: '+ str(moving_average_window)])
        else:
            axes.plot(x_nums, data, c = c)
            axes.legend(['measurement'])
        vfunc = np.vectorize(datetime.datetime.strftime)
        axes.set_xticks(np.arange(0,np.shape(time)[0],x_gap))
        axes.set_xticklabels(vfunc(time[0:np.shape(time)[0]:x_gap], "%H:%M:%S"))
        axes.margins(x=0)
        axes.set_ylabel(unit)
        
    def timeseries_plot_counts_per_s(self, fig, axes, x_gap, moving_average_window = None):
        data = self.counts_per_s
        unit = "#/s"
        self.timeseries_plot(fig, axes, data, unit, x_gap, moving_average_window = moving_average_window)

    def timeseries_plot_bin_counts_per_s(self, fig, axes, bin_num, x_gap, moving_average_window = None):
        data = self.bins_counts_per_s[:, bin_num]
        unit = "#/s"
        self.timeseries_plot(fig, axes, data, unit, x_gap, moving_average_window = moving_average_window)

    def timeseries_plot_sample_flowrate(self, fig, axes, x_gap, moving_average_window = None, diagnose = False):
        data = self.sample_flowrate
        unit = "ml/s"
        if diagnose:
            ax2 = axes.twinx()
            ax2.plot(self.oa_difference)

        self.timeseries_plot(fig, axes, data, unit, x_gap, moving_average_window = moving_average_window)


    def timeseries_plot_concentration(self, fig, axes, x_gap, moving_average_window = None):
        data = self.concentration
        unit = "#/L"
        self.timeseries_plot(fig, axes, data, unit, x_gap, moving_average_window = moving_average_window)

    def timeseries_plot_bin_concentration(self, fig, axes, bin_num, x_gap, moving_average_window = None):
        data = self.bins_concentration[:, bin_num]
        unit = "#/L"
        self.timeseries_plot(fig, axes, data, unit, x_gap, moving_average_window = moving_average_window)

    def size_distribution_timeseries_plot(self, fig, axes, data, unit, x_gap, vmin = 0.1, vmax = None, colorbar = False, norm = None):
        bin_mean = self.bin_mean
        lower_bin_boundary = self.lower_bin_boundary

        time = self.current_datetime

        X, Y = np.mgrid[0:(time.shape[0]+1), 0:(bin_mean.shape[0]+1)]
        if vmax is None:
            vmax = data.max()
        if norm == "linear":
            handle = axes.pcolor(X, Y, data, norm=colors.Normalize(vmin=vmin, vmax=vmax))
            if colorbar:
                fig.colorbar(handle, ax=axes)
        elif norm == "log":
            handle = axes.pcolor(X, Y, data, norm=colors.LogNorm(vmin=vmin, vmax=vmax))
            if colorbar:
                fig.colorbar(handle, ax=axes, extend='min')
        
        axes.set_facecolor((80/255, 25/255, 101/255))
        vfunc = np.vectorize(datetime.datetime.strftime)
        axes.set_xticks(np.arange(0,np.shape(time)[0],x_gap)+0.5)
        axes.set_xticklabels(vfunc(time[0:np.shape(time)[0]:x_gap], "%H:%M:%S"))
        axes.set_yticks(np.arange(lower_bin_boundary.shape[0]))
        axes.set_yticklabels(lower_bin_boundary)
        axes.set_ylabel("$\mathrm{\mu}$m")
    
    def size_distribution_timeseries_plot_counts_per_s(self, fig, axes, x_gap, vmin = 0.1, vmax = None, colorbar = False, norm = "linear"):
        data = self.bins_counts_per_s
        unit = "#/s"
        self.size_distribution_timeseries_plot(fig, axes, data, unit, x_gap, vmin = vmin, vmax = vmax, norm = norm)

    def size_distribution_timeseries_plot_concentration(self, fig, axes, x_gap, vmin = 0.1, vmax = None, colorbar = False, norm = "linear"):
        data = self.bins_concentration
        unit = "#/L"
        self.size_distribution_timeseries_plot(fig, axes, data, unit, x_gap, vmin = vmin, vmax = vmax, colorbar = colorbar, norm = norm)

    def size_distribution_timeseries_plot_dNdlogDP(self, fig, axes, x_gap, vmin = 0.1, vmax = None, colorbar = False, norm = "linear"):
        data = self.bins_dNdlogDP
        unit = "dNdlogDP"
        self.size_distribution_timeseries_plot(fig, axes, data, unit, x_gap, vmin = vmin, vmax = vmax, colorbar = colorbar, norm = norm)

    def vertical_line(self, fig, axes, line_time, c='r', alpha = 1):
        line_time_oa = line_time - self.OLE_TIME_ZERO
        line_time_oa = line_time_oa.total_seconds()/86400
    
        axes.axvline(np.abs(self.OA_datetime-line_time_oa).argmin(), color = c, alpha = alpha)

    def fill_rectangle(self, fig, axes, start_time, end_time, c = None, alpha = None, hatch = None, facecolor = None, edgecolor = None):
        start_time_oa = start_time - self.OLE_TIME_ZERO
        start_time_oa = start_time_oa.total_seconds()/86400
        start_idx = np.abs(self.OA_datetime-start_time_oa).argmin()

        end_time_oa = end_time - self.OLE_TIME_ZERO
        end_time_oa = end_time_oa.total_seconds()/86400
        end_idx = np.abs(self.OA_datetime-end_time_oa).argmin()

        axes.axvspan(start_idx, end_idx, color = c, alpha = alpha, hatch = hatch, facecolor = facecolor, edgecolor = edgecolor)

    def movingaverage(values, window):
        if window%2 == 0:
            raise ValueError('Use odd moving average window')
        weights = np.repeat(1.0, window)/window
        sma = np.convolve(values, weights, mode='valid')
        pre_sma = np.empty(window//2)
        pre_sma[:] = np.nan
        sma = np.hstack([pre_sma, sma, pre_sma])

        return sma

    def down_sample(self, down_sample_ratio):
        if not OPC_Measurement.is_positive_integer(down_sample_ratio):
            raise ValueError("Please use a positive integer for down_sample_ratio")
        if self.opc_data_clipped.shape[0] % down_sample_ratio !=0:
            print(self.opc_data_clipped.shape[0])
            raise ValueError("The number of rows of opc_data is not divisible by down_sample_ratio")
        
        split_num = self.opc_data_clipped.shape[0] / down_sample_ratio
        new_opc_data = np.array(np.array_split(self.opc_data_clipped, split_num, axis = 0))
        new_opc_data = np.mean(new_opc_data, axis = 1)

        new_opc_measurement = OPC_Measurement(new_opc_data, self.lower_bin_boundary, self.bin_mean)
        return new_opc_measurement
        
    def is_positive_integer(N):
        if N < 0 or N == 0:
            return False
        X = int(N)
        temp2 = N - X
        if (temp2 > 0):
            return False
        return True