% SO2 monitor (Thermo 43i) data processing
% Read raw data files, remove repetitions, calculate averages, make plots
%
% Function used: fmedia_qualquer_coisa2; frepeated; export_LFA_GoAmazon_jan2015
%
% Outputs: SO2_Level0 (mat_time(utc);y;m;d;h;m;s; so2; hiso2; intt; rctt; pgast; pres; smplfl; pmtv; flowb; lmpv; lmpi; flag(0/1))--> all raw data, without filters
% SO2_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)
% SO2_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 >=0.38LPM;
% rctt > 43oC;
% lmpi>=88% (and the subsequent point)
% -10<ConcSO2<100 ppb
%
% Criteria for Level2: correct the baseline if needed
%
% Obs: trace gas Thermo monitors report mixing ratios normalized to 1013.25 hPa and 0 oC
%
% ---------------- Parametros -----------------------
caminho='C:\A\Diadema2016\data_trailer_2016\so2_unifesp\';
caminho1='C:\A\Diadema2016\processados\SO2\';
fuso=0; % fuso horario para calcular hora local
tit='Diadema-2016';%'SO2-ZF2-2016'; %'SO2-UEA';%'SO2-TIWA';
filtflag=false; % logical para filtrar dados com flag
filtmanut=true; % logical para filtrar dados em períodos de manutenção ou com problemas
lertxt=false; % logical para ler arquivos txt com concentrações
basicplots=true;
lev1=false;
lev2=false;
%
if(filtmanut)
C1=xlsread('C:\A\Diadema2016\logbook_trailer_Diadema2016.xlsx','todos','A4:M135');
C2=xlsread('C:\A\Diadema2016\logbook_trailer_Diadema2016.xlsx','so2','A4:O37');
manut(:,1)=cat(1,datenum(C1(:,1:6)),datenum(C2(:,1:6)));
manut(:,2)=cat(1,datenum(C1(:,7:12)),datenum(C2(:,7:12)));
manut(:,3)=cat(1,C1(:,13),C2(:,13));
manut(:,4:5)=cat(1,zeros(length(C1),2),C2(:,14:15));
end
%
% inicio e termino de periodo de zoom para fazer grafico
zoom1=[datenum(2016,11,22,0,0,0) datenum(2017,3,15,0,0,0)];
zoom2=[datenum(2017,03,05,0,0,0) datenum(2017,3,15,0,0,0)];%
%
% ---------------- Fim dos parametros --------------------
%
data=zeros(1,18);
flag='--------';
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)));
vtempo=datevec(aux1,'mm-dd-yyHH:MM'); % date vector
mattime=datenum(vtempo);%-datenum(ano,1,1)+1; % matlab time
%doy=datenum(vtempo)-datenum(ano,1,1)+1; % day of year
aux2=find(strcmpi('--------',C{1,3})==0); % look for flaged data
nflag=zeros(length(mattime),1);
nflag(aux2)=1;
data1=cat(2,mattime,vtempo,aux,nflag);
data=cat(1,data,data1);
flag=cat(1,flag,C{1,3});
clear aux aux1 aux2 vtempo doy nflag C data1 nome icol mattime
fclose(fid);
end
%aux={data(:,1) flag}; % falta implementar flags --> dificil manter outdata
%e flags na mesma base temporal; talvez seja melhor transformar os diferentes flags em
%um codigo numerico
outdata=sortrows(data,1); % sort data based on day of year
outdata(1,:)=[];
clear data
%
if(lertxt)
caminhotxt='C:\A\GOAMAZON\T1\SO2_1406960668sn\raw_txt\processar\';
datatxt=zeros(1,2);
files=dir(caminhotxt);
disp('-----Reading txt data files')
for ifile=3:size(files,1)
nome=getfield(files, {ifile,1}, 'name');
disp(nome)
C=importdata(strcat(caminhotxt,nome));
aux1=strcat((C.textdata(:,1)),(C.textdata(:,2)));
datatxt1(:,1)=datenum(aux1,'yyyy-mm-ddHH:MM:SS');
datatxt1(:,2)=C.data;
datatxt=cat(1,datatxt,datatxt1);
clear aux aux1 C datatxt1 nome
end
datatxt(1,:)=[];
datatxt=sortrows(datatxt,1); % sort data based on day of year
first =floor(datatxt(1,1)); % first averaging time
last = ceil(datatxt(length(datatxt),1)); % last averaging time
delta = 5/60/24; % 5min
tol=0;
stat=fmedia_qualquer_coisa2(datatxt,first,last,delta,tol);
aux=find(stat(:,1)<=outdata(length(outdata),1)); % procura pontos que já estão em outdata (arquivos do iport)
stat(aux,:)=[];
stat=cat(2,stat,zeros(length(stat),9));
stat(:,3:18)=NaN;
stat(:,8)=stat(:,2);
stat(:,2:7)=datevec(stat(:,1));
outdata=cat(1,outdata,stat);
clear datatxt stat
end
%
if(zf2_2016)
outdata(outdata(:,1)<datenum(2016,5,3,14,40,0),:)=[]; % amostrando ar do container
end
%
disp('----Checking for repeated data')
[newdata,erro]=frepeated(outdata,8);
if (strcmp(erro,'none')~=1)
disp(erro)
aux=find(diff(outdata(:,1))==0);
datestr(outdata(aux(1),1))
end
outdata=newdata;
clear aux erro
%
% Remove maintanance periods from dataset
if(filtmanut)
outdataorig=outdata;
for i=1:size(manut,1)
aux=find(outdata(:,1)>=manut(i,1) & outdata(:,1)<=manut(i,2));
if(manut(i,3)==0)
outdata(aux,:)=[]; % remove maintanance periods from dataset
else
% detrend baseline (ozone_detrend.m)
if(manut(i,4)==0)
outdata(aux,8)=outdata(aux,8)+(-polyval(manut(i,4:5),outdata(aux,1)));
else
outdata(aux,8)=outdata(aux,8)+(polyval(manut(i,4:5),outdata(aux(1),1))-polyval(manut(i,4:5),outdata(aux,1)));
end
end
end
figure,plot(outdataorig(:,1),outdataorig(:,8),'r.'),hold on,grid on
plot(outdata(:,1),outdata(:,8),'b.')
dynamicDateTicks
clear aux
end
%
% Remove flagged periods from dataset
if(filtflag)
aux=find(outdata(:,18)==1);
outdata(aux,:)=[];
clear aux
end
%
% Corrigindo o time stamp do monitor para reportar o início do período
% O sampling time variou no tempo: 1min até 22/dez; 5min entre 22dez-11jan; e 1min depois disso.
% figure,plot(outdata(2:length(outdata),1),diff(outdata(:,1))*24*60,'b.'),grid on, dynamicDateTicks
deltat1=datenum(0,0,0,0,1,0);
deltat2=datenum(0,0,0,0,5,0);
aux=find(outdata(:,1)<datenum(2016,12,22,15,5,0));
outdata(aux,1)=outdata(aux,1)-deltat1;
aux=find(outdataorig(:,1)<datenum(2016,12,33,15,5,0));
outdataorig(aux,1)=outdataorig(aux,1)-deltat1;
aux=find(outdata(:,1)>=datenum(2016,12,22,15,4,0) & outdata(:,1)<datenum(2017,1,11,12,56,0));
outdata(aux,1)=outdata(aux,1)-deltat2;
aux=find(outdataorig(:,1)>=datenum(2016,12,22,15,4,0) & outdataorig(:,1)<datenum(2017,1,11,12,56,0));
outdataorig(aux,1)=outdataorig(aux,1)-deltat2;
aux=find(outdata(:,1)>datenum(2017,1,11,12,55,0));
outdata(aux,1)=outdata(aux,1)-deltat1;
aux=find(outdataorig(:,1)>datenum(2017,1,11,12,55,0));
outdataorig(aux,1)=outdataorig(aux,1)-deltat1;
%
% Saving output files (level0 and clean)
save([caminho1,tit,'_SO2_level0.mat'],'outdataorig','-mat')
if(filtflag || filtmanut)
save([caminho1,tit,'_SO2_clean.mat'],'outdata','-mat')
end
%
% -------------------------------------------------------------------------
% Basic Plots
if(basicplots)
% Diagnostic ---------------
figure
subplot(2,3,1),plot(outdata(:,1)+(fuso/24),outdata(:,10),'b.'),hold on
plot(outdata(:,1)+(fuso/24),outdata(:,11),'r.'),hold off
ylabel('Temp (oC)','FontSize',12),grid on, dynamicDateTicks
subplot(2,3,4),plot(outdata(:,1)+(fuso/24),outdata(:,17),'c.')
ylabel('Lamp Intensity (%)','FontSize',12),grid on, dynamicDateTicks
subplot(2,3,2),plot(outdata(:,1)+(fuso/24),outdata(:,14),'b.')
title(tit,'FontSize',12)
ylabel('Sample Flow','FontSize',12),grid on, dynamicDateTicks
subplot(2,3,5),plot(outdata(:,1)+(fuso/24),outdata(:,15),'m.')
ylabel('PMT Voltage (V)','FontSize',12),grid on, dynamicDateTicks
subplot(2,3,3),plot(outdata(:,1)+(fuso/24),outdata(:,13),'g.')
ylabel('Pressure','FontSize',12),grid on, dynamicDateTicks
subplot(2,3,6),plot(outdata(:,1)+(fuso/24),outdata(:,18),'r.'), dynamicDateTicks
ylabel('Flag','FontSize',12),grid on
saveas(gcf,strcat(caminho1,tit,'_SO2_diag'),'jpg')
%
% Time series and diurnal cycle --------------
figure
subplot(2,1,1),plot(outdata(:,1)+(fuso/24),outdata(:,8),'b.')
title(tit,'FontSize',14)
ylabel('SO2 [ppb]','FontSize',14),grid on, dynamicDateTicks
subplot(2,1,2),plot(mod(outdata(:,1)+(fuso/24),1),outdata(:,8),'b.')
ylabel('SO2 [ppb]','FontSize',14),grid on,xlabel('UTC Time','FontSize',14), axis tight
set(gca,'XTick',[0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0])
set(gca,'XTickLabel',{'3:00';'6:00';'9:00';'12:00';'15:00';'18:00';'21:00';'0:00'})
saveas(gcf,strcat(caminho1,tit,'_SO2_conc'),'jpg')
%
% Zoom --------------
aux=find(outdata(:,1)>=zoom1(1) & outdata(:,1)<=zoom1(2));
figure
subplot(2,1,1),plot(outdata(aux,1)+(fuso/24),outdata(aux,8),'b.')
title(tit,'FontSize',14)
ylabel('SO2 [ppb]','FontSize',14),grid on%, dynamicDateTicks
set(gca,'XTick',zoom1(1):2:zoom1(2))
set(gca,'XTickLabel',datestr(zoom1(1):2:zoom1(2),'dd/mm'))
set(gca,'XLim',[zoom1(1) zoom1(2)])
%set(gca,'YLim',[-0.5 2])
subplot(2,1,2),plot(mod(outdata(aux,1)+(fuso/24),1),outdata(aux,8),'b.')
ylabel('SO2 [ppb]','FontSize',14),grid on,xlabel('UTC Time','FontSize',14), axis tight
set(gca,'XTick',[0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0])
set(gca,'XTickLabel',{'3:00';'6:00';'9:00';'12:00';'15:00';'18:00';'21:00';'0:00'})
%set(gca,'YLim',[-0.5 2])
saveas(gcf,strcat(caminho1,tit,'_SO2_conc_zoom1'),'jpg')
%
% Zoom 2 --------------
aux=find(outdata(:,1)>=zoom2(1) & outdata(:,1)<=zoom2(2));
figure
subplot(2,1,1),plot(outdata(aux,1)+(fuso/24),outdata(aux,8),'b.')
title(tit,'FontSize',14)
ylabel('SO2 [ppb]','FontSize',14),grid on%, dynamicDateTicks
set(gca,'XTick',zoom2(1):1:zoom2(2))
set(gca,'XTickLabel',datestr(zoom2(1):1:zoom2(2),'dd/mm'))
set(gca,'XLim',[zoom2(1) zoom2(2)])
%set(gca,'YLim',[-0.5 6])
subplot(2,1,2),plot(mod(outdata(aux,1)+(fuso/24),1),outdata(aux,8),'b.')
ylabel('SO2 [ppb]','FontSize',14),grid on,xlabel('UTC Time','FontSize',14), axis tight
set(gca,'XTick',[0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0])
set(gca,'XTickLabel',{'3:00';'6:00';'9:00';'12:00';'15:00';'18:00';'21:00';'0:00'})
%set(gca,'YLim',[-0.5 6])
saveas(gcf,strcat(caminho1,tit,'_SO2_conc_zoom2'),'jpg')
%
% Ciclo diurno médio
first = 0; % first averaging time
last = 23.5/24; % last averaging time
delta = 30/60/24; % 30min
mdata=cat(2,mod(outdata(:,1),1),outdata(:,8));
tol=0;
stat=fmedia_qualquer_coisa2(mdata,first,last,delta,tol);
figure
errorbar(stat(:,1),stat(:,2),stat(:,3),'b.')
dynamicDateTicks
set(gca,'XLim',[0 1])
ylabel('SO_2 (ppb)','FontSize',14),grid on,xlabel('UTC Time','FontSize',14)
title('Average Diurnal Cycle','FontSize',14)
hgsave(gcf,strcat(caminho1,tit,'_SO2_dcycle_alldata'))
saveas(gcf,strcat(caminho1,tit,'_SO2_dcycle_alldata'),'jpg')
end
%
% -------------------------------------------------------------------------
% Level 1
% Keep records when:
% Flow >=0.38LPM; rctt > 43oC; lmpi>=88% (and the subsequent point); -10<ConcSO2<100 ppb
% Remove calibration periods according to the field logbook.
%
if(lev1)
C=load('C:\A\Diadema2016\processados\SO2\Diadema-2016_SO2_clean.mat');
outdata=C.outdata;
C=load('C:\A\Diadema2016\processados\SO2\Diadema-2016_SO2_level0.mat');
outdataorig=C.outdataorig;
tit='Diadema-2016';
%
origsize=size(outdataorig,1);
aux=find(outdata(:,14)<0.38 | outdata(:,11)<43 | outdata(:,17)<88);
aux1=aux(diff(aux)>1)+1; % subsequent point
aux2=find(outdata(:,8)<=-10 | outdata(:,8)>=100);
rem=cat(1,aux,aux1,aux2);
outdata(rem,:)=[];
clear aux aux1 aux2 rem
figure,plot(outdataorig(:,1),outdataorig(:,8),'r.'),grid on, hold on
plot(outdata(:,1),outdata(:,8),'b.'),dynamicDateTicks
%
% Save Level 1 data ----------
Instrument='Thermo';
Model='43iTLE';
serial_num='1333859828';
Site=tit;
Val_level='1';
Proc_vers='1.0';
labels{1}='SO2';
units{1}='(ppb)';
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,outdata(:,1),outdata(:,8))
end
%
%
% -------------------------------------------------------------------------
% Level 2
% Criteria for Level2: correct the baseline when needed
if(lev2)
% Lê os arquivos Level 1 que já salvou
aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_Thermo_43iTLE_Level1_5min.csv');
level1_5min=aux.data(:,1:2);
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)));
tit='Diadema-2016';
%
% Save Level 2 data
Instrument='Thermo';
Model='43iTLE';
serial_num='1333859828';
Site=tit;
Val_level='2';
Proc_vers='1.0';
labels{1}='SO2';
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
%
%
% Verificação final
% aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_Thermo_43iTLE_Level1_5min.csv');
% level1_5min=aux.data(:,1:2);
% level1_5min(:,1)=level1_5min(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
% aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_Thermo_43iTLE_Level2_5min.csv');
% level2_5min=aux.data(:,1:2);
% level2_5min(:,1)=level2_5min(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
% figure, plot(level1_5min(:,1),level1_5min(:,2),'b.'),grid on, hold on,dynamicDateTicks
% plot(level2_5min(:,1),level2_5min(:,2),'r.')