% Ozone monitor (Thermo 49i) data processing
% Read raw data files, remove repetitions, calculate averages, make plots
% Functions used: frepeated; export_LFA_GoAmazon_jan2015
% Outputs: O3_Level0 (mat_time(utc); o3; hio3; cellai; cellbi; bncht; lmpt; o3lt; flowa; flowb; pres; flag(0/1))--> all raw data, without filters
% O3_Level1_5min (dd-mmm-yyyy HH:MM; sec_since_1970; mean_o3; std_o3; median_o3; p25_o3; p75_o3; max_o3; min_o3; num_o3)
% O3_Level1_30min
% Time stamps in the output files refer to the start of the period (note that Thermo gas monitors time stamps refer to the end of the period)
% Criteria for Level1: (Keep records when) Flow A and B >= 0.55 LPM; Intensity A and B >= 75000 Hz (and the subsequent point)
% -50<ConcO3<1000 ppb
% Criteria for Level2: remove calibration and maintanace periods and eventually correct the baseline
%
% Obs: trace gas Thermo monitors report mixing ratios normalized to 1013.25 hPa and 0 oC
% Para converter de unix time para matlabtime: datestr(level1_5min(1,1)/(24*60*60)+datenum(1970,1,1))
% ---------------- Parametros -----------------------
caminho='F:\O3\Raw_Data_RB\';
caminho1='C:\Users\milen\OneDrive\PosDoc\RioBranco\O3\';
fuso=5; % fuso horario para calcular hora local
tit='RioBranco-2018';
level2=false; % logical to produce level2 data
ano=2018;
%
% ---------------- Fim dos parametros --------------------
%
data=zeros(1,12);
%
files=dir(caminho);
disp('-----Reading raw data files')
for ifile=3:size(files,1)
nome=getfield(files, {ifile,1}, 'name');
disp(nome)
fid=fopen(strcat(caminho,nome));
textscan(fid,'%s',1);
textscan(fid,'%s %s %s %s %s %s %s %s %s',1);
textscan(fid,'%s',1);
textscan(fid,'%s %s %s %s %s %s %s',1);
textscan(fid,'%s',1);
label=textscan(fid,'%s %s %s %s %s %s %s %s %s %s %s %s %s',1);
C=textscan(fid,'%s %s %s %f %f %f %f %f %f %f %f %f %f');
for icol=1:size(C,2)-3
aux(:,icol)=C{1,icol+3};
end
aux1=strcat((C{1,2}(:,1)),(C{1,1}(:,1)));
aux2=datevec(aux1,'mm-dd-yyHH:MM'); % date vector
tmat=datenum(aux2); % date number in Matlab
aux3=find(strcmpi('--------',C{1,3})==0); % look for flaged data
nflag=zeros(length(tmat),1);
nflag(aux3)=1;
data1=cat(2,tmat,aux,nflag);
data=cat(1,data,data1);
flag=cat(1,flag,C{1,3});
clear aux aux1 aux2 aux3 nflag C data1 nome icol
fclose(fid);
end
outdata=sortrows(data,1); % sort data based on day of year
outdata(1,:)=[];
clear data
outdata=cat(2,outdata(:,1),datevec(outdata(:,1)),outdata(:,2:12));
%
disp('----Checking for repeated data')
[newdata,erro,linha]=frepeated(outdata,8);
if (strcmp(erro,'none')~=1)
disp(erro)
disp(linha)
else
clear outdata
outdata=newdata;
clear newdata erro
end
% Corrigindo o time stamp do monitor para reportar o início do período
aux=diff(outdata(1:2,1));
outdata(2:length(outdata),1)=outdata(2:length(outdata),1)-diff(outdata(:,1));
outdata(1,1)=outdata(1,1)-aux;
clear aux
outdata(outdata(:,1)<datenum(ano,1,1),:)=[];
save([caminho1,tit,'_O3_Level0.txt'],'outdata','-ASCII')
%
orig=outdata;
%
figure
subplot(2,3,1),plot(outdata(:,1),outdata(:,10),'b.')
d1=outdata(1,1); d2=outdata(length(outdata),1);
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Intensity Cell A (Hz)','FontSize',12),grid on, axis tight
subplot(2,3,4),plot(outdata(:,1),outdata(:,11),'c.')
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Intensity Cell B (Hz)','FontSize',12),grid on, axis tight
subplot(2,3,2),plot(outdata(:,1),outdata(:,15),'b.')
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
title(tit,'FontSize',12)
ylabel('Flow Cell A','FontSize',12),grid on, axis tight
subplot(2,3,5),plot(outdata(:,1),outdata(:,16),'m.')
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Flow Cell B','FontSize',12),grid on, axis tight
subplot(2,3,3),plot(outdata(:,1),outdata(:,17),'g.')
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Pressure','FontSize',12),grid on, axis tight
subplot(2,3,6),plot(outdata(:,1),outdata(:,18),'r.'), axis tight
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Flag','FontSize',12),grid on
saveas(gcf,strcat(caminho1,'O3_diag'),'png')
%
figure
plot(orig(:,1),orig(:,8),'r.'),grid on, hold on,dynamicDateTicks
%
% ------ Filters for Level1 data ---------------------------
origsize=size(outdata,1);
aux=find(outdata(:,15)<0.55 | outdata(:,16)<0.55 | outdata(:,10)<60000 | outdata(:,11)<60000);
aux1=aux(diff(aux)>1)+1; % subsequent point
aux2=find(outdata(:,8)<=-50 | outdata(:,8)>=1000);
rem=cat(1,aux,aux1,aux2);
outdata(rem,:)=[];
clear aux aux1 aux2 rem
%
% ------- Additional periods to remove (specific for Rio Branco) ---------------
rem=[datenum(2018,08,24,19,41,00) datenum(2018,08,24,20,19,00);... % zero check
datenum(2018,09,02,19,27,00) datenum(2018,09,02,20,15,00);...
datenum(2018,09,10,19,53,00) datenum(2018,09,10,20,44,00);... % inicializando
datenum(2018,09,11,19,47,00) datenum(2018,09,15,23,21,00);...
datenum(2018,09,15,23,05,00) datenum(2018,09,15,23,21,00);...
datenum(2018,09,19,21,36,00) datenum(2018,09,19,22,40,00);...
datenum(2018,10,10,21,35,00) datenum(2018,10,15,19,37,00);...
datenum(2018,10,30,15,23,00) datenum(2018,10,30,16,19,00);...
datenum(2018,10,30,21,15,00) datenum(2018,11,05,23,36,00);];

for i=1:length(rem)
outdata(outdata(:,1)>=rem(i,1) & outdata(:,1)<=rem(i,2),:)=[];
end
%
figure
plot(orig(:,1),orig(:,8),'r.'),grid on, hold on,dynamicDateTicks
plot(outdata(:,1),outdata(:,8),'b.')
plot([outdata(1,1) outdata(length(outdata),1)],[0 0],'-y')
legend('removed','keeped');
saveas(gcf,strcat(caminho1,'O3_removed_points'),'png')

%
% ----------------- Save Level 1 data -----------------------------------------
Instrument='Thermo';
Model='49i';
serial_num='1406960636 30.192';
Site=tit;
Val_level='1';
Proc_vers='1.0';
labels{1}='O3';
units{1}='(ppb)';
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,outdata(:,1),outdata(:,8))
%
%%
%----------------- Level 2 processing -----------------------------------------
if(level2)
% Lê os arquivos Level 1 5 min
% no PC LFA:
% C:\Users\milen\OneDrive\PosDoc\RioBranco\0_Ascii-files\Level1\
% file: RioBranco-2018_Thermo_49i_Level1_5min.csv
aux=importdata('C:\Users\milen\OneDrive\PosDoc\RioBranco\0_Ascii-files\Level1\RioBranco-2018_Thermo_49i_Level1_5min.csv');
level1_5min=aux.data;
clear aux
level1_5min(:,1)=level1_5min(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
origsize=sum(~isnan(level1_5min(:,2)));
orig=level1_5min;
%
% Detrending baseline accoring to zero checks
% manut(:,1)=manut(:,1)-mod(manut(:,1),5/60/24); % arredondando horários de manutenção para 5 minutos inteiros
% manut(:,2)=manut(:,2)-mod(manut(:,2),5/60/24);
% for i=1:size(manut,1)
% aux=find(level1_5min(:,1)>=manut(i,1) & level1_5min(:,1)<=manut(i,2));
% if(manut(i,3)==0 && ~isempty(aux))
% level1_5min(aux,2:size(level1_5min,2))=NaN; % remove maintanance periods from dataset
% elseif(~isempty(aux))
% % detrend baseline (ozone_detrend.m)
% if(manut(i,4)==0)
% level1_5min(aux,2)=level1_5min(aux,2)-manut(i,5);
% else
% level1_5min(aux,2)=level1_5min(aux,2)+(polyval(manut(i,4:5),level1_5min(aux(1),1)))+(-polyval(manut(i,4:5),level1_5min(aux,1)));
% end
% end
% clear aux aux1
% end
figure
plot(orig(:,1),orig(:,2),'r.'),grid on, hold on
plot(level1_5min(:,1),level1_5min(:,2),'b.'),dynamicDateTicks
%
% Save Level 2 data
Instrument='Thermo';
Model='49i';
serial_num='1406960636 30.192';
Site=tit;
Val_level='2';
Proc_vers='1.0';
labels{1}='O3';
units{1}='(ppb)';
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level1_5min(:,1),level1_5min(:,2))
end
%
%