# Module for Data Analysis of CWFT Experiments # Created by Nicole Leemann (with ChatGPT as Coding Assistant) # Atmospheric Chemistry Lab Course 2026 Edition import numpy as np import pandas as pd import plotly.express as px import plotly.graph_objects as go from plotly.subplots import make_subplots from pathlib import Path from babel.dates import format_date from scipy.optimize import curve_fit import re import ast # some independent helper functions (for background fits) def linfit(x, p1, p0): return p1 * x + p0 def poly2fit(x, p2, p1, p0): return p2 * x**2 +p1 * x + p0 def log_decay(x, p1, p0): return p1*np.log(x) + p0 class experiment: """ A class for the data analysis of CWFT experiments. Input Parameters ---------- exp_name : string Name of the experiment (used for plot titles) instrument_list : string list a list of all instruments used in the experiment (BTOF not included), e.g. ["T360U_CO2", "T200_NOx", "Ecophysics_NOx_899", "ICAD_HONO", "ICAD_NO2", "Humidity_Sensors"] setup_dict: dictionary dictionary containing setup parameters like temperatures, pressure, flows and tube dimensions, format: { "p_room": 975, # ambient pressure [mbar] "T_room": 23, # room temperature [°C] "T_tube": 18, # CWFT temperature [°C] "T_hum": 15, # humidifier temperature [°C] "flow_N2": 1300, # MFC set flow for N2 "flow_O2": 300, # MFC set flow for O2 "dilution": 3000, # MFC set flow for N2 dilution "dil_BTOF": 1000, # MFC set flow for BTOF dilution "tube_length": 32.6, # tube length [cm] "tube_diam": 1.64, # tube inner diameter [cm] "vol_coat": 1000 # volume of coating solution [muL] } exp_start : string start date & time of the experiment, format: YYYY-MM-DD hh:mm:ss exp_end : string start date & time of the experiment, format: YYYY-MM-DD hh:mm:ss time_tube_in : string time when the flow was switched to the tube path, format: YYYY-MM-DD hh:mm:ss time_lights_on : string time when the lights were switched on, format: YYYY-MM-DD hh:mm:ss time_lights_off : string time when the lights were switched off, format: YYYY-MM-DD hh:mm:ss time_tube_out : string time when the flow was switched back to the bypass path, format: YYYY-MM-DD hh:mm:ss time_lights_on_2 : string, default: None time when the lights were switched on for the second time, format: YYYY-MM-DD hh:mm:ss time_lights_off_2 : string, default: None time when the lights were switched off for the second time, format: YYYY-MM-DD hh:mm:ss color_dict: dictionary dictionary containing lists of color hexcodes used for plotting each instrument, format: { "T360U_CO2": ["#474647"], # CO2 "T200_NOx": ["#80B4E5", "#F08080", "#C5A1F5"], # NO, NO2, NOx "Ecophysics_NOx_899": ["#345E85", "#A14242", "#644EA0"], # NO, NO2, NOx "ICAD_HONO": ["#9a66b3", "#91ddcd"], # HONO, NO2 "ICAD_NO2": ["#6cdf72"] # NO2 } savefolder: path folder where the results are saved """ def __init__(self, exp_name, instrument_list, setup_dict, exp_start, exp_end, time_tube_in, time_lights_on, time_lights_off, time_tube_out, color_dict, savefolder, time_lights_on_2=None, time_lights_off_2=None): ## Inputs self.exp_name = exp_name self.savefolder = savefolder self.instruments = instrument_list self.exp_start = pd.Timestamp(exp_start) self.exp_end = pd.Timestamp(exp_end) # approx time (can be adjusted later) self.time_tube_in = pd.Timestamp(time_tube_in) self.time_lights_on = pd.Timestamp(time_lights_on) self.time_lights_off = pd.Timestamp(time_lights_off) self.time_tube_out = pd.Timestamp(time_tube_out) if time_lights_on_2 != None: # if second irradiation time self.time_lights_on_2 = pd.Timestamp(time_lights_on_2) self.time_lights_off_2 = pd.Timestamp(time_lights_off_2) else: self.time_lights_on_2 = None self.time_lights_off_2 = None # color dictionary for plotting self.color_dict = color_dict # setup parameters self.setup_dict = setup_dict self.RH = 80 # RH in flow tube -> will be calculated and replaced later self.dilution_factor = 1 # will be calculated later ## parameters set to defaults / to be filled later # data files (to be filled) self.co2_file = None self.nox_file = None self.nox_eco_file = None self.icad_hono_file = None self.icad_no2_file = None self.rh_file = None # timeshift of instruments (align times to CO2 analyzer) self.timeshift_nox = pd.Timedelta(minutes=0, seconds=0) self.timeshift_nox_eco = pd.Timedelta(minutes=56, seconds=40) self.timeshift_icad_hono = pd.Timedelta(minutes=60, seconds=37) self.timeshift_icad_no2 = pd.Timedelta(minutes=60, seconds=0) self.timeshift_rh = pd.Timedelta(minutes=0, seconds=0) # background correction / parameters of linear fit (to be filled) self.bkgcorr_co2 = [np.array([0,0])] self.bkgcorr_nox = [np.array([0,0]),np.array([0,0]),np.array([0,0])] self.bkgcorr_nox_eco = [np.array([0,0]),np.array([0,0]),np.array([0,0])] self.bkgcorr_icad_hono = [np.array([0,0]),np.array([0,0])] self.bkgcorr_icad_no2 = [np.array([0,0])] # time intervals used for background correction (to be filled) self.bkgint_co2 = [None,None] self.bkgint_nox = [None,None] self.bkgint_nox_eco = [None,None] self.bkgint_icad_hono = [None,None] self.bkgint_icad_no2 = [None,None] # dataframes of raw data self.rawdata_co2 = None self.rawdata_nox = None self.rawdata_nox_eco = None self.rawdata_icad_hono = None self.rawdata_icad_no2 = None self.rawdata_rh = None # dataframes of preprocessed / cut data self.predata_co2 = None self.predata_nox = None self.predata_nox_eco = None self.predata_icad_hono = None self.predata_icad_no2 = None self.predata_rh = None # processed dataframes self.data_co2 = None self.data_nox = None self.data_nox_eco = None self.data_icad_hono = None self.data_icad_no2 = None self.data_rh = None # Constants self.T_std = 273.15 # [K] standard temperature #self.P_std = 1013.25 # [mbar] standard pressure # Antoine equation parameters for vapor pressure (Bridgeman & Aldrich, 1964) from https://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Mask=4&Type=ANTOINE&Plot=on self.A, self.B, self.C = 5.40221, 1838.675, -31.737 # ------------------------------------------------------------------------- # Functions for RH & flow calculation # ------------------------------------------------------------------------- def get_pH2O(self, T): """Calculate water vapor pressure at temperature T [°C].""" log_p = self.A - self.B / (self.C + T + self.T_std) # Antoine equation return 10**log_p def get_RH(self): """Calculate relative humidity in the flow tube.""" self.RH = self.get_pH2O(self.setup_dict["T_hum"]) / self.get_pH2O(self.setup_dict["T_tube"]) * 100 return self.RH def get_dilution_factor(self): """Calculate dilution factor at the instruments.""" flow_tube = self.setup_dict["flow_N2"] + self.setup_dict["flow_O2"] flow_tot = flow_tube + self.setup_dict["dilution"] self.dilution_factor = flow_tube / flow_tot #return self.dilution_factor # ------------------------------------------------------------------------- # Functions for loading & reshaping data # ------------------------------------------------------------------------- # def timestamp_from_filename(fname: Path) -> pd.Timestamp: # """extract timestamp from filename (by ChatGPT)""" # ts_str = fname.name.split("_instrument")[0] # need to implement this for all instruments... # return pd.to_datetime(ts_str, format="%Y-%m-%d_%H.%M.%S") # def find_files(self, exp_start, exp_end, instrument="T360U_CO2"): # """find all files belonging to the experiment for one instrument (based on code by ChatGPT)""" # path = f"H:/Data/{instrument}" # base_path = ( # Path(path) # / f"{exp_start.year }" # / f"{exp_start.month:02d}_{exp_start.strftime("%B")}" # note: not all instruments follow this system... # ) # files = list(base_path.glob(f"{exp_start.strftime("%Y-%m-%d")}_*_instrument.dat")) # find all files from the same day (what to do when exp overnight...) # matching_files = [ # f for f in files # if exp_start <= self.timestamp_from_filename(f) <= exp_end # ] # return matching_files ### assuming for now that there is one file per experiment day for each instrument (will need to correct for this later...) def get_rawdata_T360U_CO2(self, ind=0): """Get the rawdata of T360U_CO2. If multiple files found: load file with index ind. Default: first file""" #file_co2 = f"{path_co2}/2025/11_November/2025-11-07_13.27.38_PC14750_T360U_CO2.dat" #year, month = self.exp_start.year, self.exp_start.month #monthname = format_date(self.exp_start, format="MMMM", locale="de") base_path = Path("ACLC26_Data/T360U_CO2") date = self.exp_start.strftime("%Y-%m-%d") files = list(base_path.glob(f"{date}_*.dat")) # find all relevant files files = sorted(files) self.co2_file = files if not files: raise FileNotFoundError("No instrument files found for this date") if len(files) > 1: print("Multiple files found") #, loading first one") #all of them") print(files) file = files[ind] # take only first file for now, will deal with multiple files later self.rawdata_co2 = pd.read_csv(file, sep = '\t', header=1) # def get_rawdata_T200_NOx(self, ind=0): # """Get the rawdata of T200_NOx.""" # #file_nox = f"{path_nox}/2025/11_November/2025-11-07_13.27.28_PC14750_T200_NOx.dat" # year, month = self.exp_start.year, self.exp_start.month # monthname = format_date(self.exp_start, format="MMMM", locale="de") # path = "H:/Data/T200_NOx" # base_path = ( # Path(path) # / f"{year }" # / f"{month:02d}_{monthname}" # ) # date = self.exp_start.strftime("%Y-%m-%d") # files = list(base_path.glob(f"{date}_*.dat")) # files = sorted(files) # self.nox_file = files # if not files: # raise FileNotFoundError("No instrument files found for this date") # if len(files) > 1: # print("Multiple files found") #, loading first one") #all of them") # print(files) # file = files[ind] # take only first file for now, will deal with multiple files later # self.rawdata_nox = pd.read_csv(file, sep = '\t', header=1) def get_rawdata_RH(self, ind=0): """Get the rawdata of RH sensor.""" #file_nox = f"{path_nox}/2025/11_November/2025-11-07_13.26.51_PC14750_Humidity_Sensors.dat" #year, month = self.exp_start.year, self.exp_start.month #monthname = format_date(self.exp_start, format="MMMM", locale="de") base_path = Path("ACLC26_Data/Humidity_Sensors") date = self.exp_start.strftime("%Y-%m-%d") files = list(base_path.glob(f"{date}_*.dat")) files = sorted(files) self.rh_file = files if not files: raise FileNotFoundError("No instrument files found for this date") if len(files) > 1: print("Multiple files found") #, loading first one") #all of them") print(files) file = files[ind] # take only first file for now, will deal with multiple files later self.rawdata_rh = pd.read_csv(file, sep = '\t', header=1, encoding= 'unicode_escape') def get_rawdata_ICAD_HONO(self): """Get the rawdata of ICAD HONO.""" #file_hono = f"{path_hono}/2025_11_07_Channel1_Data_HONO.dat" #year, month, monthname = self.exp_start.year, self.exp_start.month, self.exp_start.strftime("%B") base_path = Path("ACLC26_Data/ICAD_HONO") date = self.exp_start.strftime("%Y_%m_%d") files = list(base_path.glob(f"{date}_*.dat")) self.icad_hono_file = files file = files[0] self.rawdata_icad_hono = pd.read_csv(file, sep = '\t') def get_rawdata_ICAD_NO2(self): """Get the rawdata of ICAD NO2.""" #file_no2 = f"{path_no2}/2025_11_07_Channel1_Data_NO2.dat" #year, month, monthname = self.exp_start.year, self.exp_start.month, self.exp_start.strftime("%B") base_path = Path("ACLC26_Data/ICAD_NO2") date = self.exp_start.strftime("%Y_%m_%d") files = list(base_path.glob(f"{date}_*.dat")) self.icad_no2_file = files file = files[0] self.rawdata_icad_no2 = pd.read_csv(file, sep = '\t') def get_rawdata_NOx_eco(self): """Get the rawdata of Ecophysics_NOx_899.""" #file_nox_eco = f"{path_nox_eco}/251107_FeCit_Nitrate_rep_NO_NOx_NO2_65_1.txt" #year, month, monthname = self.exp_start.year, self.exp_start.month, self.exp_start.strftime("%B") base_path = Path("ACLC26_Data/Ecophysics_NOx_899") date = self.exp_start.strftime("%Y%m%d")[2:] files = list(base_path.glob(f"{date}_*_1.txt")) self.nox_eco_file = files if len(files) > 1: for f in files: fs = str(f) # convert path to string try: recstart = pd.to_datetime(fs[fs.index(date):fs.index(date)+11], format='%y%m%d_%H%M') # find & convert datetime except: print("No time in filename. Loading first file.") file = files[0] break if self.exp_start < recstart < self.exp_end: file = f break else: file = files[0] self.rawdata_nox_eco = pd.read_csv(file, sep = '\t', skiprows=26) # ------------------------------------------------------------------------- # Data Pre-Processing & Preview Plots # ------------------------------------------------------------------------- def preprocess_data(self, instr): """Create Date/Time column for data and cut unnecessary data. Saves in new dataframe predata_[instr].""" if instr == "T360U CO2": self.rawdata_co2["Date/Time"] = pd.to_datetime(self.rawdata_co2["Date"] + " " + self.rawdata_co2["Time"], # Add date/time column format='%d.%m.%Y %H:%M:%S') # keep only data points from the experiment mask = ( (self.rawdata_co2["Date/Time"] >= self.exp_start) & (self.rawdata_co2["Date/Time"] <= self.exp_end) ) self.predata_co2 = self.rawdata_co2.copy().loc[mask] self.predata_co2 = self.predata_co2[["Date/Time","CO2"]] # keep only relevant data self.predata_co2.columns = ["Date/Time", "raw CO2 [ppm]"] # rename the columns elif instr == "ICAD HONO": self.rawdata_icad_hono["Date/Time"] = pd.to_datetime(self.rawdata_icad_hono["Start Date/Time (UTC)"], # Add date/time column format='%Y-%m-%d %H:%M:%S.%f') mask = ( (self.rawdata_icad_hono["Date/Time"] >= self.exp_start) & (self.rawdata_icad_hono["Date/Time"] <= self.exp_end) ) self.predata_icad_hono = self.rawdata_icad_hono.copy().loc[mask] self.predata_icad_hono = self.predata_icad_hono[["Date/Time","HONO (ppb)", "HONO Uncertainty (ppb)","NO2 (ppb)","NO2 Uncertainty (ppb)"]] # keep only relevant data self.predata_icad_hono.columns = ["Date/Time", "raw HONO [ppb]", "raw HONO Uncertainty [ppb]", "raw NO2 [ppb]", "raw NO2 Uncertainty [ppb]"] # rename the columns elif instr == "ICAD NO2": self.rawdata_icad_no2["Date/Time"] = pd.to_datetime(self.rawdata_icad_no2["Start Date/Time (UTC)"], # Add date/time column format='%Y-%m-%d %H:%M:%S.%f') mask = ( (self.rawdata_icad_no2["Date/Time"] >= self.exp_start) & (self.rawdata_icad_no2["Date/Time"] <= self.exp_end) ) self.predata_icad_no2 = self.rawdata_icad_no2.copy().loc[mask] self.predata_icad_no2 = self.predata_icad_no2[["Date/Time","NO2 (ppb)","NO2 Uncertainty (ppb)"]] # keep only relevant data self.predata_icad_no2.columns = ["Date/Time", "raw NO2 [ppb]", "raw NO2 Uncertainty [ppb]"] # rename the columns elif instr == "T200 NOx": self.rawdata_nox["Date/Time"] = pd.to_datetime(self.rawdata_nox["Date"] + " " + self.rawdata_nox["Time"], # Add date/time column format='%d.%m.%Y %H:%M:%S') mask = ( (self.rawdata_nox["Date/Time"] >= self.exp_start) & (self.rawdata_nox["Date/Time"] <= self.exp_end) ) self.predata_nox = self.rawdata_nox.copy().loc[mask] self.predata_nox = self.predata_nox[["Date/Time","NO","NO2","NOx"]] # keep only relevant data self.predata_nox.columns = ["Date/Time", "raw NO [ppb]", "raw NO2 [ppb]", "raw NOx [ppb]"] # rename the columns elif instr == "Ecophysics NOx 899": self.rawdata_nox_eco["Date/Time"] = pd.to_datetime(self.rawdata_nox_eco["TimeStamp"], # Add date/time column format='%d/%m/%Y %H:%M:%S.%f') mask = ( (self.rawdata_nox_eco["Date/Time"] >= self.exp_start) & (self.rawdata_nox_eco["Date/Time"] <= self.exp_end) ) self.predata_nox_eco = self.rawdata_nox_eco.copy().loc[mask] self.predata_nox_eco = self.predata_nox_eco[["Date/Time","NO, ppb","NO2, ppb","NOx, ppb"]] # keep only relevant data self.predata_nox_eco.columns = ["Date/Time", "raw NO [ppb]", "raw NO2 [ppb]", "raw NOx [ppb]"] # rename the columns elif instr == "Humidity Sensors": self.rawdata_rh["Date/Time"] = pd.to_datetime(self.rawdata_rh["Date"] + " " + self.rawdata_rh["Time"], format='%d.%m.%Y %H:%M:%S') mask = ( (self.rawdata_rh["Date/Time"] >= self.exp_start) & (self.rawdata_rh["Date/Time"] <= self.exp_end) ) self.predata_rh = self.rawdata_rh.copy().loc[mask] self.predata_rh = self.predata_rh[["Date/Time","RH5[%]","T5[°C]","DP5[°C]"]] # get only relevant data self.predata_rh.columns = ["Date/Time","RH [%]","T [°C]","DP [°C]"] # rename the columns def preview_plot(self, instr): """Preview Plot to align time series & irradiation time for instruments in instr list.""" #timeshift_co2 = pd.Timedelta(hours=1, seconds=30) #fig = go.Figure() fig = make_subplots(specs=[[{"secondary_y": True}]]) fig.update_layout(template="plotly_white", title_text=self.exp_name) # Irradiation -> add rectangle to mark it fig.add_shape( type="rect", x0=self.time_lights_on, x1=self.time_lights_off, y0=0, y1=1, xref="x1", yref="paper", fillcolor="gold", opacity=0.3, line_width=0, layer="below", ) if self.time_lights_on_2 != None: # add rectangle for second irradiation period fig.add_shape( type="rect", x0=self.time_lights_on_2, x1=self.time_lights_off_2, y0=0, y1=1, xref="x1", yref="paper", fillcolor="gold", opacity=0.3, line_width=0, layer="below", ) # Vertical lines at tube in and tube out times fig.add_shape( type="line", x0=self.time_tube_in, x1=self.time_tube_in, y0=0, y1=1, xref="x1", yref="paper", line=dict( color="darkgray", width=2, dash="dash" ), layer="below" ) fig.add_shape( type="line", x0=self.time_tube_out, x1=self.time_tube_out, y0=0, y1=1, xref="x1", yref="paper", line=dict( color="darkgray", width=2, dash="dash" ), layer="below" ) ## Add traces # CO2 self.preprocess_data(instr="T360U CO2") fig.add_trace(go.Scatter(x=self.predata_co2["Date/Time"], y=self.predata_co2["raw CO2 [ppm]"], # + self.bkgcorr_co2, mode='lines', name='CO2', line=dict(color=self.color_dict["T360U_CO2"][0])), secondary_y=False ) if "ICAD_HONO" in instr: self.preprocess_data(instr="ICAD HONO") fig.add_trace(go.Scatter(x=self.predata_icad_hono["Date/Time"] + self.timeshift_icad_hono, y=self.predata_icad_hono["raw HONO [ppb]"], # + self.bkgcorr_icad_hono, mode='lines', name='HONO', line=dict(color=self.color_dict["ICAD_HONO"][0])), secondary_y=True ) fig.add_trace(go.Scatter(x=self.predata_icad_hono["Date/Time"] + self.timeshift_icad_hono, y=self.predata_icad_hono["raw NO2 [ppb]"], # + self.bkgcorr_icad_hono, mode='lines', name='NO2 (HONO)', line=dict(color=self.color_dict["ICAD_HONO"][1])), secondary_y=True ) if "ICAD_NO2" in instr: self.preprocess_data(instr="ICAD NO2") fig.add_trace(go.Scatter(x=self.predata_icad_no2["Date/Time"] + self.timeshift_icad_no2, y=self.predata_icad_no2["raw NO2 [ppb]"], # + self.bkgcorr_icad_no2, mode='lines', name='NO2', line=dict(color=self.color_dict["ICAD_NO2"][0])), secondary_y=True ) if "T200_NOx" in instr: self.preprocess_data(instr="T200 NOx") fig.add_trace(go.Scatter(x=self.predata_nox["Date/Time"] - self.timeshift_nox, y=self.predata_nox["raw NO [ppb]"], # + self.bkgcorr_nox, mode='lines', name='NO (T200)', line=dict(color=self.color_dict["T200_NOx"][0])), secondary_y=True ) fig.add_trace(go.Scatter(x=self.predata_nox["Date/Time"] - self.timeshift_nox, y=self.predata_nox["raw NO2 [ppb]"], # + self.bkgcorr_nox, mode='lines', name='NO2 (T200)', line=dict(color=self.color_dict["T200_NOx"][1])), secondary_y=True ) fig.add_trace(go.Scatter(x=self.predata_nox["Date/Time"] - self.timeshift_nox, y=self.predata_nox["raw NOx [ppb]"], # + self.bkgcorr_nox, mode='lines', name='NOx (T200)', line=dict(color=self.color_dict["T200_NOx"][2])), secondary_y=True ) if "Ecophysics_NOx_899" in instr: self.preprocess_data(instr="Ecophysics NOx 899") fig.add_trace(go.Scatter(x=self.predata_nox_eco["Date/Time"] - self.timeshift_nox_eco, y=self.predata_nox_eco["raw NO [ppb]"], # + self.bkgcorr_nox_eco, mode='lines', name='NO (eco)', line=dict(color=self.color_dict["Ecophysics_NOx_899"][0])), secondary_y=True ) fig.add_trace(go.Scatter(x=self.predata_nox_eco["Date/Time"] - self.timeshift_nox_eco, y=self.predata_nox_eco["raw NO2 [ppb]"], # + self.bkgcorr_nox, mode='lines', name='NO2 (eco)', line=dict(color=self.color_dict["Ecophysics_NOx_899"][1])), secondary_y=True ) fig.add_trace(go.Scatter(x=self.predata_nox_eco["Date/Time"] - self.timeshift_nox_eco, y=self.predata_nox_eco["raw NOx [ppb]"], # + self.bkgcorr_nox, mode='lines', name='NOx (eco)', line=dict(color=self.color_dict["Ecophysics_NOx_899"][2])), secondary_y=True ) if "Humidity_Sensors" in instr: self.preprocess_data(instr="Humidity Sensors") fig.add_trace(go.Scatter(x=self.predata_rh["Date/Time"] - self.timeshift_rh, y=self.predata_rh["RH [%]"], mode='lines', name='RH [%]') , secondary_y=True ) # axes and layout fig.update_xaxes(type="date") # make sure it is a time axis fig.update_yaxes(title_text="CO2 [ppm]", secondary_y=False) fig.update_yaxes(title_text="concentration [ppb]", secondary_y=True) #fig.show() return fig def RHTDP_plot(self, num_timeaxis=False): """Plot RH, Temperature and DP. If num_timeaxis=True the time axis will be Time since Irradiation [min].""" self.preprocess_data(instr="Humidity Sensors") if num_timeaxis: self.get_num_timeaxis(self.predata_rh, ts=self.timeshift_rh) xdata = self.predata_rh["Time since Irradiation [min]"] else: xdata = self.predata_rh["Date/Time"] fig = make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.03) #subplot_titles=("relative humidity [%]", "temperature [°C]", "dew point temperature [°C]")) fig.add_trace( go.Scatter(x=xdata, y=self.predata_rh["RH [%]"], mode='lines', name='RH', #line=dict(color="#345E85") ), row=1, col=1 ) fig.add_trace( go.Scatter(x=xdata, y=self.predata_rh["T [°C]"], mode='lines', name='T', #line=dict(color="#345E85") ), row=2, col=1 ) fig.add_trace( go.Scatter(x=xdata, y=self.predata_rh["DP [°C]"], mode='lines', name='DP', #line=dict(color="#345E85") ), row=3, col=1 ) fig.update_yaxes(title_text="RH [%]", range=[10, 25], row=1, col=1) fig.update_yaxes(title_text="T [°C]", range=[20, 30], row=2, col=1) fig.update_yaxes(title_text="DP [°C]", range=[-5, 5], row=3, col=1) if num_timeaxis: fig.update_xaxes(title_text="Time since Irradiation [min]", row=3, col=1) else: fig.update_xaxes(title_text="Time", row=3, col=1) fig.update_layout(height=600, showlegend=False) # , width=1800 #, title_text="Side By Side Subplots") #fig.show() return fig # ------------------------------------------------------------------------- # Background Fitting # ------------------------------------------------------------------------- # fit dunctions defined before class def get_num_timeaxis(self, df_data, ts, ts_sign ="-"): """Corrects for the instrument timeshift ts with sign + or - and converts the time axis (df_data["Date/Time]) from datatime to numeric time and save in new column df_data["Time since Irradiation [min]"].""" if ts_sign == "-": df_data["Date/Time"] = df_data["Date/Time"] - ts if ts_sign == "+": # ICADs df_data["Date/Time"] = df_data["Date/Time"] + ts df_data["Time since Irradiation [min]"] = (df_data["Date/Time"] - self.time_lights_on).dt.total_seconds() / 60 def get_bkg_fit(self, df_data, pred_start, pred_end, param, mode = "lin"): """Linear Fit to determine background of parameter param (str of colname) in df_data between Timestamps pred_start and pred_end. mode="lin" for linear fit, mode="poly" for second order polynomial fit, mode="log_decay" for logarithmic de. Make sure that the numeric timeaxis column already exists. Returns a plot to check the results and a Fit Summary (?)""" # subset data for fit mask = ( (df_data["Date/Time"] >= pred_start) & (df_data["Date/Time"] <= pred_end) ) df_subset = df_data.copy().loc[mask] # convert to numpy x_data_np = df_subset["Time since Irradiation [min]"].to_numpy() y_data_np = df_subset[param].to_numpy() if mode == "lin": popt, pcov = curve_fit(linfit, x_data_np, y_data_np) bkg_pred = linfit(df_data["Time since Irradiation [min]"], *popt) label = f'{param} ' + 'fit: %5.3f * x + %5.3f' % tuple(popt) residuals = df_data[param].to_numpy() - bkg_pred if mode == "poly": popt, pcov = curve_fit(poly2fit, x_data_np, y_data_np) bkg_pred = poly2fit(df_data["Time since Irradiation [min]"], *popt) label = f'{param} ' + 'fit: %5.3f * x^2 + %5.3f *x + %5.3f' % tuple(popt) residuals = df_data[param].to_numpy() - bkg_pred if mode == "log_decay": x_fit = x_data_np + np.abs(df_data["Time since Irradiation [min]"][0]) + 1 # log decay only valid for x>0 x_pred = np.asarray(df_data["Time since Irradiation [min]"]) + np.abs(df_data["Time since Irradiation [min]"][0]) + 1 popt, pcov = curve_fit(log_decay, x_fit, y_data_np) bkg_pred = log_decay(x_pred, *popt) label = f'{param} ' + 'fit: %5.3f * x + %5.3f' % tuple(popt) residuals = df_data[param].to_numpy() - bkg_pred fig = go.Figure() fig.update_layout(template="plotly_white") # Vertical lines at tube in and tube out times fig.add_shape( type="line", x0=(self.time_tube_in - self.time_lights_on).total_seconds() / 60, x1=(self.time_tube_in - self.time_lights_on).total_seconds() / 60, y0=0, y1=1, xref="x1", yref="paper", line=dict( color="red", width=2, dash="dash" ), layer="below" ) fig.add_shape( type="line", x0=(self.time_tube_out - self.time_lights_on).total_seconds() / 60, x1=(self.time_tube_out - self.time_lights_on).total_seconds() / 60, y0=0, y1=1, xref="x1", yref="paper", line=dict( color="red", width=2, dash="dash" ), layer="below" ) # Add traces fig.add_trace(go.Scatter(x=df_data["Time since Irradiation [min]"], y=df_data[param], mode='lines', name='data', line=dict(color="darkgrey", dash="dot")) ) fig.add_trace(go.Scatter(x=x_data_np, y=y_data_np, mode='lines', name='data for fit', line=dict(color="black")) ) fig.add_trace(go.Scatter(x=df_data["Time since Irradiation [min]"], y=bkg_pred, mode='lines', name='bkg prediction', line=dict(color="#DD5800")) ) fig.update_layout(title_text = label) fig.update_xaxes(title_text="Time since Irradiation [min]") fig.update_yaxes(title_text="concentration") return fig, popt, pcov, residuals # ------------------------------------------------------------------------- # Data Processing # ------------------------------------------------------------------------- def process_data(self, savefolder, savename): """Finish processing data (remove background) & save csv file for each instrument's data in folder savefolder.""" self.get_dilution_factor() print("dilution factor:", self.dilution_factor) if "T360U_CO2" in self.instruments: self.data_co2 = self.predata_co2.copy() bkg = self.bkgcorr_co2[0][0] * self.data_co2["Time since Irradiation [min]"] + self.bkgcorr_co2[0][1] # get background self.data_co2["CO2 [ppm]"] = self.data_co2["raw CO2 [ppm]"] - bkg # correct for background self.data_co2["CO2 [ppm]"] = self.data_co2["CO2 [ppm]"] / self.dilution_factor # correct for dilution self.data_co2.to_csv(f"{savefolder}/Processed_Data/{savename}_Data_T360U_CO2.csv", sep=';') # save data in csv if "T200_NOx" in self.instruments: self.data_nox = self.predata_nox.copy() bkg_no = self.bkgcorr_nox[0][0] * self.data_nox["Time since Irradiation [min]"] + self.bkgcorr_nox[0][1] # get background self.data_nox["NO [ppb]"] = self.data_nox["raw NO [ppb]"] - bkg_no # correct for background bkg_no2 = self.bkgcorr_nox[1][0] * self.data_nox["Time since Irradiation [min]"] + self.bkgcorr_nox[1][1] self.data_nox["NO2 [ppb]"] = self.data_nox["raw NO2 [ppb]"] - bkg_no2 bkg_nox = self.bkgcorr_nox[2][0] * self.data_nox["Time since Irradiation [min]"] + self.bkgcorr_nox[2][1] self.data_nox["NOx [ppb]"] = self.data_nox["raw NOx [ppb]"] - bkg_nox self.data_nox[["NO [ppb]", "NO2 [ppb]", "NOx [ppb]"]] = self.data_nox[["NO [ppb]", "NO2 [ppb]", "NOx [ppb]"]] / self.dilution_factor # correct for dilution self.data_nox.to_csv(f"{savefolder}/Processed_Data/{savename}_Data_T200_NOx.csv", sep=';') if "Ecophysics_NOx_899" in self.instruments: self.data_nox_eco = self.predata_nox_eco.copy() bkg_no = self.bkgcorr_nox_eco[0][0] * self.data_nox_eco["Time since Irradiation [min]"] + self.bkgcorr_nox_eco[0][1] self.data_nox_eco["NO [ppb]"] = self.data_nox_eco["raw NO [ppb]"] - bkg_no bkg_no2 = self.bkgcorr_nox_eco[1][0] * self.data_nox_eco["Time since Irradiation [min]"] + self.bkgcorr_nox_eco[1][1] self.data_nox_eco["NO2 [ppb]"] = self.data_nox_eco["raw NO2 [ppb]"] - bkg_no2 bkg_nox = self.bkgcorr_nox_eco[2][0] * self.data_nox_eco["Time since Irradiation [min]"] + self.bkgcorr_nox_eco[2][1] self.data_nox_eco["NOx [ppb]"] = self.data_nox_eco["raw NOx [ppb]"] - bkg_nox self.data_nox_eco[["NO [ppb]", "NO2 [ppb]", "NOx [ppb]"]] = self.data_nox_eco[["NO [ppb]", "NO2 [ppb]", "NOx [ppb]"]] / self.dilution_factor # correct for dilution self.data_nox_eco.to_csv(f"{savefolder}/Processed_Data/{savename}_Data_Ecophysics_NOx_899.csv", sep=';') if "ICAD_HONO" in self.instruments: self.data_icad_hono = self.predata_icad_hono.copy() bkg_hono = self.bkgcorr_icad_hono[0][0] * self.data_icad_hono["Time since Irradiation [min]"] + self.bkgcorr_icad_hono[0][1] self.data_icad_hono["HONO [ppb]"] = self.data_icad_hono["raw HONO [ppb]"] - bkg_hono self.data_icad_hono["HONO Uncertainty [ppb]"] = self.data_icad_hono["raw HONO Uncertainty [ppb]"] - bkg_hono bkg_no2 = self.bkgcorr_icad_hono[1][0] * self.data_icad_hono["Time since Irradiation [min]"] + self.bkgcorr_icad_hono[1][1] self.data_icad_hono["NO2 [ppb]"] = self.data_icad_hono["raw NO2 [ppb]"] - bkg_no2 self.data_icad_hono["NO2 Uncertainty [ppb]"] = self.data_icad_hono["raw NO2 Uncertainty [ppb]"] - bkg_no2 self.data_icad_hono[["HONO [ppb]", "HONO Uncertainty [ppb]", "NO2 [ppb]", "NO2 Uncertainty [ppb]"]] = self.data_icad_hono[["HONO [ppb]", "HONO Uncertainty [ppb]", "NO2 [ppb]", "NO2 Uncertainty [ppb]"]] / self.dilution_factor # correct for dilution self.data_icad_hono.to_csv(f"{savefolder}/Processed_Data/{savename}_Data_ICAD_HONO.csv", sep=';') if "ICAD_NO2" in self.instruments: self.data_icad_no2 = self.predata_icad_no2.copy() bkg_no2 = self.bkgcorr_icad_no2[0][0] * self.data_icad_no2["Time since Irradiation [min]"] + self.bkgcorr_icad_no2[0][1] self.data_icad_no2["NO2 [ppb]"] = self.data_icad_no2["raw NO2 [ppb]"] - bkg_no2 self.data_icad_no2["NO2 Uncertainty [ppb]"] = self.data_icad_no2["raw NO2 Uncertainty [ppb]"] - bkg_no2 self.data_icad_no2[["NO2 [ppb]", "NO2 Uncertainty [ppb]"]] = self.data_icad_no2[["NO2 [ppb]", "NO2 Uncertainty [ppb]"]] / self.dilution_factor # correct for dilution self.data_icad_no2.to_csv(f"{savefolder}/Processed_Data/{savename}_Data_ICAD_NO2.csv", sep=';') if "Humidity_Sensors" in self.instruments: self.data_rh = self.predata_rh.copy() # already in final form print("Note that RH and DP are not dilution corrected.") self.data_rh.to_csv(f"{savefolder}/Processed_Data/{savename}_Data_Humidity_Sensors.csv", sep=';') # ------------------------------------------------------------------------- # Plotting Functions # ------------------------------------------------------------------------- def get_final_data(self, instr): """Returns final processed data of instrument instr.""" if instr == "T360U_CO2": return self.data_co2.copy() if instr == "T200_NOx": return self.data_nox.copy() if instr == "Ecophysics_NOx_899": return self.data_nox_eco.copy() if instr == "ICAD_HONO": return self.data_icad_hono.copy() if instr == "ICAD_NO2": return self.data_icad_no2.copy() def plot_one_instrument(self, instr, data = None, plot_raw = False): """Plot all variables measured by the chosen instrument. If data = None, the data will be loaded from the class object. If plot_raw = True, also plot the rawdata.""" if data is None: data = self.get_final_data(instr) fig = go.Figure() fig.update_layout(template="plotly_white") self.add_timemarkers(fig) unit = "ppb" # unit for all parameters except CO2 # colors for the different instruments (maybe have this as input dictionary...) if instr == "T360U_CO2": params = ["CO2"] colors = self.color_dict["T360U_CO2"] unit = "ppm" # change the unit group = "T360 CO2" if instr == "T200_NOx": params = ["NO", "NO2", "NOx"] colors = self.color_dict["T200_NOx"] group = "T200 NOx" if instr == "Ecophysics_NOx_899": params = ["NO", "NO2", "NOx"] colors = self.color_dict["Ecophysics_NOx_899"] group = "Ecophysics NOx 899" if instr == "ICAD_HONO": params = ["HONO", "NO2"] colors = self.color_dict["ICAD_HONO"] group = "ICAD HONO" if instr == "ICAD_NO2": params = ["NO2"] colors = self.color_dict["ICAD_NO2"] group = "ICAD NO2" # Add traces for n, param in enumerate(params): if plot_raw: fig.add_trace(go.Scatter(x=data["Time since Irradiation [min]"], y=data[f"raw {param} [{unit}]"], mode='lines', name=f'{param} Raw Data', legendgroup=group, legendgrouptitle_text=group, line=dict(color=colors[n], dash="dot"), opacity=0.7) ) fig.add_trace(go.Scatter(x=data["Time since Irradiation [min]"], y=data[f"{param} [{unit}]"], mode='lines', name=f'{param} Data', legendgroup=group, legendgrouptitle_text=group, line=dict(color=colors[n])) ) if instr == "ICAD_HONO" or instr == "ICAD_NO2": upper = data[f"{param} [{unit}]"] + data[f"{param} Uncertainty [{unit}]"] lower = data[f"{param} [{unit}]"] - data[f"{param} Uncertainty [{unit}]"] fig.add_trace(go.Scatter( x = pd.concat([data["Time since Irradiation [min]"], data["Time since Irradiation [min]"][::-1]]), # x, then x reversed y = pd.concat([upper, lower[::-1]]), # upper, then lower reversed name = f"{param} Uncertainty", legendgroup=group, legendgrouptitle_text=group, fill='toself', opacity = 0.2, fillcolor = colors[n], #'rgba(0,100,80,0.2)', line=dict(color='rgba(255,255,255,0)'), hoverinfo="skip", )) fig.update_layout(title_text = self.exp_name) fig.update_xaxes(title_text="Time since Irradiation [min]") fig.update_yaxes(title_text=f"Concentration [{unit}]") # # dummy legend for irradiation # fig.add_trace( # go.Scatter( # x=[None], # y=[None], # mode="markers", # marker=dict( # size=12, # symbol="square", # color="gold", # opacity=0.3 # ), # name="Irradiation", # legendgroup="Time Markers", # showlegend=True # ) # ) return fig def add_timemarkers(self, fig): """Add vertical lines at tube in/out times and rectangles for irradiation to figure fig (for numeric time axis).""" # Irradiation -> add rectangle to mark it fig.add_shape( type="rect", x0=0, x1=(self.time_lights_off - self.time_lights_on).total_seconds() / 60, y0=0, y1=1, xref="x1", yref="paper", fillcolor="gold", opacity=0.3, line_width=0, layer="below" ) if self.time_lights_on_2 != None: # add rectangle for second irradiation period fig.add_shape( type="rect", x0=(self.time_lights_on_2 - self.time_lights_on).total_seconds() / 60, x1=(self.time_lights_off_2 - self.time_lights_on).total_seconds() / 60, y0=0, y1=1, xref="x1", yref="paper", fillcolor="gold", opacity=0.3, line_width=0, layer="below", ) # Vertical lines at tube in and tube out times fig.add_shape( type="line", x0=(self.time_tube_in - self.time_lights_on).total_seconds() / 60, x1=(self.time_tube_in - self.time_lights_on).total_seconds() / 60, y0=0, y1=1, xref="x1", yref="paper", line=dict( color="darkgrey", width=2, dash="dash" ), layer="below" ) fig.add_shape( type="line", x0=(self.time_tube_out - self.time_lights_on).total_seconds() / 60, x1=(self.time_tube_out - self.time_lights_on).total_seconds() / 60, y0=0, y1=1, xref="x1", yref="paper", line=dict( color="darkgrey", width=2, dash="dash" ), layer="below" ) return fig def plot_multiple(self, title, instrs, figures=None): """Plot traces from several instruments (in list instrs). If figures already exist input a list of figures, otherwise figures=None. Title will show up in the figure title after the experiment name.""" # get figures if they don't exist already if figures == None: figures = [] for i in instrs: figures.append(self.plot_one_instrument(i)) # make the figure fig = go.Figure() fig.update_layout(template="plotly_white") self.add_timemarkers(fig) # add irradiation rectangles and vertical lines for tube in/out # add traces for f in figures: for trace in f.data: fig.add_trace(trace) fig.update_layout(title_text = self.exp_name + f" - {title}") fig.update_xaxes(title_text="Time since Irradiation [min]") fig.update_yaxes(title_text=f"Concentration [ppb]") return fig def plot_multiple_2axes(self, title, instrs1, instrs2, figures1=None, figures2=None): """Plot traces from several instruments in one figure with two y-axes. List instrs1 will use the left yaxis, instrs2 the right yaxis. If figures1 & figures2 already exist input a lists of figures, otherwise figures=None. Title will show up in the figure title after the experiment name.""" # get figures if they don't exist already if figures1 == None: figures1 = [] for i in instrs1: figures1.append(self.plot_one_instrument(i)) if figures2 == None: figures2 = [] for j in instrs2: figures2.append(self.plot_one_instrument(j)) # make the figure fig = make_subplots(specs=[[{"secondary_y": True}]]) fig.update_layout(template="plotly_white") self.add_timemarkers(fig) # add irradiation rectangles and vertical lines for tube in/out # add traces to axis 1 for f1 in figures1: for trace in f1.data: fig.add_trace(trace, secondary_y=False) # add traces to axis 2 for f2 in figures2: for trace in f2.data: fig.add_trace(trace, secondary_y=True) fig.update_layout(title_text = self.exp_name + f" - {title}") fig.update_xaxes(title_text="Time since Irradiation [min]") if instrs1==["T360U_CO2"]: fig.update_yaxes(title_text=f"CO2 [ppm]", secondary_y=False) else: if "T360U_CO2" in instrs1: print("WARNING: Plotting CO2 and other data on the same y-axis. -> Always put CO2 data separately or the y-axis label is wrong!") fig.update_yaxes(title_text=f"Concentration [ppb]", secondary_y=False) if instrs2==["T360U_CO2"]: fig.update_yaxes(title_text=f"CO2 [ppm]", secondary_y=True) else: if "T360U_CO2" in instrs2: print("WARNING: Plotting CO2 and other data on the same y-axis. -> Always put CO2 data separately or the y-axis label is wrong!") fig.update_yaxes(title_text=f"Concentration [ppb]", secondary_y=True) return fig # ------------------------------------------------------------------------- # Summary Function # ------------------------------------------------------------------------- def get_summary(self, savefile): """Saves a summary of the data analysis in a file.""" # open file f = open(savefile, 'w', encoding="utf-8") # write summary f.write('_________________________________________________________________________\n' + 'SUMMARY DATA ANALYSIS\n' + '_________________________________________________________________________\n\n') f.write(f"Experiment Name: {self.exp_name}\n") f.write(f"Results saved in: {self.savefolder}\n\n") f.write("SETUP PARAMETERS\n-------------------------------------\n") f.write(f"relative humidity {self.get_RH():.2f} %\n" + f"humidifier temperature {self.setup_dict['T_hum']} °C\n" + f"CWFT temperature {self.setup_dict['T_tube']} °C\n" + f"ambient temperature {self.setup_dict['T_room']} °C\n" + f"ambient pressure {self.setup_dict['p_room']} mbar\n" + f"-------------------------------------\n" + f"N2 flow {self.setup_dict['flow_N2']} mL/min\n" + f"O2 flow {self.setup_dict['flow_O2']} mL/min\n" + f"dilution flow (N2) {self.setup_dict['dilution']} mL/min\n" + f"dilution flow BTOF (N2) {self.setup_dict['dil_BTOF']} mL/min\n" + f"-------------------------------------\n" + f"tube length {self.setup_dict['tube_length']} cm\n" + f"tube inner diameter {self.setup_dict['tube_diam']} cm\n" + f"volume coating {self.setup_dict['vol_coat']} muL\n" ) f.write("-------------------------------------\n") f.write(f"Instrument List: {self.instruments}\n\n") f.write("TIMELINE\n--------------------------------------\n") f.write( f"experiment start {self.exp_start}\n" + f"tube in {self.time_tube_in}\n" + f"lights on {self.time_lights_on}\n" + f"lights off {self.time_lights_off}\n" + f"lights on 2 {self.time_lights_on_2}\n" + f"lights off 2 {self.time_lights_off_2}\n" + f"tube out {self.time_tube_out}\n" + f"experiment end {self.exp_end}\n" ) f.write("--------------------------------------\n\n") f.write("TIME ALIGNMENT TIMESHIFTS (ABSOLUTE)\n" + "---------------------------------------------\n" + "INSTRUMENTS TIMESHIFTS (ABSOLUTE)\n" + "---------------------------------------------\n" ) ts_nox = self.timeshift_nox.components ts_nox_eco = self.timeshift_nox_eco.components ts_hono = self.timeshift_icad_hono.components ts_no2 = self.timeshift_icad_no2.components ts_rh =self.timeshift_rh.components f.write( f"T360U CO2 REFERENCE\n" + f"T200 NOx {ts_nox.hours} h {ts_nox.minutes:02d} min {ts_nox.seconds:02d} s\n" + f"Ecophysics NOx 899 {ts_nox_eco.hours} h {ts_nox_eco.minutes:02d} min {ts_nox_eco.seconds:02d} s\n" + f"ICAD HONO {ts_hono.hours} h {ts_hono.minutes:02d} min {ts_hono.seconds:02d} s\n" + f"ICAD NO2 {ts_no2.hours} h {ts_no2.minutes:02d} min {ts_no2.seconds:02d} s\n" + f"Humidity Sensors {ts_rh.hours} h {ts_rh.minutes:02d} min {ts_rh.seconds:02d} s\n" ) f.write("---------------------------------------------\n\n") f.write("INSTRUMENT INTERVAL BACKGROUND CORRECTIONS (LIN. FIT)\n" + "---------------------------------------------------------------------------\n" + " PARAMETER SLOPE INTERCEPT \n" + "---------------------------------------------------------------------------\n" ) f.write( f"T360U CO2 {self.bkgint_co2}\n" + f" CO2 {self.bkgcorr_co2[0]} ppm\n" + f"T200 NOx {self.bkgint_nox}\n" + f" NO {self.bkgcorr_nox[0]} ppb\n" + f" NO2 {self.bkgcorr_nox[1]} ppb\n" + f" NOx {self.bkgcorr_nox[2]} ppb\n" + f"Ecophysics NOx 899 {self.bkgint_nox_eco}\n" + f" NO {self.bkgcorr_nox_eco[0]} ppb\n" + f" NO2 {self.bkgcorr_nox_eco[1]} ppb\n" + f" NOx {self.bkgcorr_nox_eco[2]} ppb\n" + f"ICAD HONO {self.bkgint_icad_hono}\n" + f" HONO {self.bkgcorr_icad_hono[0]} ppb\n" + f" NO2 {self.bkgcorr_icad_hono[1]} ppb\n" + f"ICAD NO2 {self.bkgint_icad_no2}\n" + f" NO2 {self.bkgcorr_icad_no2[0]} ppb\n" ) f.write("---------------------------------------------------------------------------\n\n") f.write("ANALYSED DATA FILES\n-----------------------------------------------------------------------------------------------------------\n") f.write(f"{self.co2_file}\n{self.nox_file}\n{self.nox_eco_file}\n{self.icad_hono_file}\n{self.icad_no2_file}\n{self.rh_file}\n") f.write("-----------------------------------------------------------------------------------------------------------\n\n") # close file f.close() ### initialize the class from a summary file ### def find_summary_blocks(savefolder, savename): """Find the lines where the blocks of the summary file start and return them.""" summary_file = f"{savefolder}/{savename}_DA_Summary.txt" # summary of data analysis file f = open(summary_file, 'r', encoding="utf-8") for n, line in enumerate(f, start=1): # go through all lines if line.lstrip().startswith("SETUP PARAMETERS"): n_setup = n elif line.lstrip().startswith("TIMELINE"): n_timeline = n elif line.lstrip().startswith("TIME ALIGNMENT"): n_alignment = n elif line.lstrip().startswith("INSTRUMENT"): # background correction n_bkg = n elif line.lstrip().startswith("ANALYSED DATA FILES"): n_files = n f.close() return [n_setup, n_timeline, n_alignment, n_bkg, n_files] def init_from_summary(savefolder, savename, color_dict): """Reads in parameters from summary file (savefolder/savename_DA_Summary.txt) and initialises the experiment class objects with these parameters.""" summary_file = f"{savefolder}/{savename}_DA_Summary.txt" # summary of data analysis file f = open(summary_file, 'r', encoding="utf-8") lines = list(f) n_setup, n_timeline, n_alignment, n_bkg, n_files = find_summary_blocks(savefolder, savename) # find starting line of the summary blocks # create dictionary / parameters to be assigned setup_dict = {"p_room": 1, "T_room": 1, "T_tube": 1, "T_hum": 1, "flow_N2": 1, "flow_O2": 1, "dilution": 1, "dil_BTOF": 1, "tube_length": 1, "tube_diam": 1, "vol_coat": 1} rh = 1 ts_nox, ts_nox_eco, ts_hono, ts_no2, ts_rh = [],[],[],[],[] bkg_interval_co2, bkg_interval_nox, bkg_interval_nox_eco, bkg_interval_hono, bkg_interval_no2 = [],[],[],[],[] bkg_corr_co2, bkg_corr_nox, bkg_corr_nox_eco, bkg_corr_hono, bkg_corr_no2 = [],[],[],[],[] co2_file, nox_file, nox_eco_file, hono_file, no2_file = '', '', '', '', '' # assign so the code will not fail if an instrument is missing number = r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?" # pattern for numbers (optional sign, decimal point and scientific notation) timestamp = r'[\d]{4}-[\d]{2}-[\d]{2} [\d]{2}:[\d]{2}:[\d]{2}' # pattern for timestamps YYYY-MM-DD hh:mm:ss for n, line in enumerate(lines, start=1): # go through all lines if line.lstrip().startswith("Experiment Name"): exp_name = line[line.index(":")+2:].rstrip() if line.lstrip().startswith("Results saved in"): # no need to import path because it is an input of the function savefolder2 = Path(line[line.index(":")+2:].rstrip()) if savefolder != savefolder2: print("Note: The savefolder given in the input (summary file location) doesn't match the savefolder path in the summary file!") # SETUP PARAMETERS if n > n_setup and n < n_timeline-2: if line.lstrip().startswith("relative humidity"): rh = float(re.search(number, line).group()) # extract the number from the line if line.lstrip().startswith("humidifier temperature"): setup_dict["T_hum"] = float(re.search(number, line).group()) if line.lstrip().startswith("CWFT temperature"): setup_dict["T_tube"] = float(re.search(number, line).group()) if line.lstrip().startswith("ambient temperature"): setup_dict["T_room"] = float(re.search(number, line).group()) if line.lstrip().startswith("ambient pressure"): setup_dict["p_room"] = float(re.search(number, line).group()) if line.lstrip().startswith("N2 flow"): setup_dict["flow_N2"] = float(re.findall(number, line)[1]) # get second number, ignore 2 in N2 if line.lstrip().startswith("O2 flow"): setup_dict["flow_O2"] = float(re.findall(number, line)[1]) if line.lstrip().startswith("dilution flow (N2)"): setup_dict["dilution"] = float(re.findall(number, line)[1]) if line.lstrip().startswith("dilution flow BTOF"): setup_dict["dil_BTOF"] = float(re.findall(number, line)[1]) if line.lstrip().startswith("tube length"): setup_dict["tube_length"] = float(re.search(number, line).group()) if line.lstrip().startswith("tube inner diameter"): setup_dict["tube_diam"] = float(re.search(number, line).group()) if line.lstrip().startswith("volume coating"): setup_dict["vol_coat"] = float(re.search(number, line).group()) # Instrument List if line.lstrip().startswith("Instrument List"): instrument_list = ast.literal_eval(line[line.index(":")+2:].rstrip()) # TIMELINE if n > n_timeline and n < n_alignment: if line.lstrip().startswith("experiment start"): exp_start = re.search(timestamp, line).group() if line.lstrip().startswith("tube in"): time_tube_in = re.search(timestamp, line).group() if line.lstrip().startswith("lights on "): time_lights_on = re.search(timestamp, line).group() if line.lstrip().startswith("lights off "): time_lights_off = re.search(timestamp, line).group() if line.lstrip().startswith("tube out"): time_tube_out = re.search(timestamp, line).group() if line.lstrip().startswith("experiment end"): exp_end = re.search(timestamp, line).group() if line.lstrip().startswith("lights on 2"): try: time_lights_on_2 = re.search(timestamp, line).group() except: time_lights_on_2 = None if line.lstrip().startswith("lights off 2"): try: time_lights_off_2 = re.search(timestamp, line).group() except: time_lights_off_2 = None # ALIGNMENT if n > n_alignment and n < n_bkg: if line.lstrip().startswith("T200"): ts_nox = np.array(re.findall(number, line)[1:], dtype=float) # ignore numbers in instrument name if line.lstrip().startswith("Ecophysics NOx"): ts_nox_eco = np.array(re.findall(number, line)[1:], dtype=float) if line.lstrip().startswith("ICAD HONO"): ts_hono = np.array(re.findall(number, line)[:], dtype=float) if line.lstrip().startswith("ICAD NO2"): ts_no2 = np.array(re.findall(number, line)[1:], dtype=float) if line.lstrip().startswith("Humidity Sensors"): ts_rh = np.array(re.findall(number, line)[:], dtype=float) # BACKGROUND if n > n_bkg and n < n_files: if line.lstrip().startswith("T360U"): bkg_interval_co2 = ast.literal_eval(line[line.index("["):].rstrip()) line1 = lines[n].rstrip() bkg_corr_co2.append(np.asarray(re.findall(number, line1)[1:], dtype=float)) if line.lstrip().startswith("T200"): bkg_interval_nox = ast.literal_eval(line[line.index("["):].rstrip()) line1, line2, line3 = lines[n].rstrip(), lines[n+1].rstrip(), lines[n+2].rstrip() bkg_corr_nox.append(np.asarray(re.findall(number, line1)[:], dtype=float)) bkg_corr_nox.append(np.asarray(re.findall(number, line2)[1:], dtype=float)) bkg_corr_nox.append(np.asarray(re.findall(number, line3)[:], dtype=float)) if line.lstrip().startswith("Ecophysics NOx"): bkg_interval_nox_eco = ast.literal_eval(line[line.index("["):].rstrip()) line1, line2, line3 = lines[n].rstrip(), lines[n+1].rstrip(), lines[n+2].rstrip() bkg_corr_nox_eco.append(np.asarray(re.findall(number, line1)[:], dtype=float)) bkg_corr_nox_eco.append(np.asarray(re.findall(number, line2)[1:], dtype=float)) bkg_corr_nox_eco.append(np.asarray(re.findall(number, line3)[:], dtype=float)) if line.lstrip().startswith("ICAD HONO"): bkg_interval_hono = ast.literal_eval(line[line.index("["):].rstrip()) line1, line2 = lines[n].rstrip(), lines[n+1].rstrip() bkg_corr_hono.append(np.asarray(re.findall(number, line1)[:], dtype=float)) bkg_corr_hono.append(np.asarray(re.findall(number, line2)[1:], dtype=float)) if line.lstrip().startswith("ICAD NO2"): bkg_interval_no2 = ast.literal_eval(line[line.index("["):].rstrip()) line1 = lines[n].rstrip() bkg_corr_no2.append(np.asarray(re.findall(number, line1)[1:], dtype=float)) # FILES if n > n_files: if "T360U_CO2" in line: co2_file = line[line.index("'")+1:-4].rstrip() if "T200_NOx" in line: nox_file = line[line.index("'")+1:-4].rstrip() if "Ecophysics_NOx_899" in line: nox_eco_file = line[line.index("'")+1:-4].rstrip() if "ICAD/HONO" in line: hono_file = line[line.index("'")+1:-4].rstrip() if "ICAD/NO2" in line: no2_file = line[line.index("'")+1:-4].rstrip() f.close() # initialise the experiment class FeCitN = experiment(exp_name, instrument_list, setup_dict, exp_start, exp_end, time_tube_in, time_lights_on, time_lights_off, time_tube_out, color_dict, savefolder, time_lights_on_2, time_lights_off_2) # assign extracted values FeCitN.RH = rh FeCitN.timeshift_nox = pd.Timedelta(hours=ts_nox[0] ,minutes=ts_nox[1], seconds=ts_nox[2]) FeCitN.timeshift_nox_eco = pd.Timedelta(hours=ts_nox_eco[0] ,minutes=ts_nox_eco[1], seconds=ts_nox_eco[2]) FeCitN.timeshift_icad_hono = pd.Timedelta(hours=ts_hono[0] ,minutes=ts_hono[1], seconds=ts_hono[2]) FeCitN.timeshift_icad_no2 = pd.Timedelta(hours=ts_no2[0] ,minutes=ts_no2[1], seconds=ts_no2[2]) FeCitN.timeshift_rh = pd.Timedelta(hours=ts_rh[0] ,minutes=ts_rh[1], seconds=ts_rh[2]) FeCitN.bkgint_co2 = bkg_interval_co2 FeCitN.bkgint_nox = bkg_interval_nox FeCitN.bkgint_nox_eco = bkg_interval_nox_eco FeCitN.bkgint_icad_hono = bkg_interval_hono FeCitN.bkgint_icad_no2 = bkg_interval_no2 FeCitN.bkgcorr_co2 = bkg_corr_co2 FeCitN.bkgcorr_nox = bkg_corr_nox FeCitN.bkgcorr_nox_eco = bkg_corr_nox_eco FeCitN.bkgcorr_icad_hono = bkg_corr_hono FeCitN.bkgcorr_icad_no2 = bkg_corr_no2 FeCitN.co2_file = Path(co2_file) FeCitN.nox_file = Path(nox_file) FeCitN.nox_eco_file = Path(nox_eco_file) FeCitN.icad_hono_file = Path(hono_file) FeCitN.icad_no2_file = Path(no2_file) return FeCitN