% OPS-TSI data processing
%
% O programa assume que os Dp bins são os mesmos do 1o arquivo.
% Do manual, entendi que os arquivos csv fornecem a contagem de partículas
% em cada bin. As contagens no bin1 correspondem ao número de partículas
% com bin1<Dp<bin2. As contagens no último bin referem-se a todas as
% partículas com Dp acima de 10um.
% Não ficou claro se essas contagens já estão em #/cm3 ou em #.
% Pág 5.4: entendi que nos arquivos csv as contagens estão em #, ou seja, precisamos dividir pelo fluxo.
% (as dimensões apresentadas na fórmula do manual estão inconsistentes)
% raw(#/cm3) = N(#)/(Q(cm3/s)*(ts-DTC*td))
% onde Q=16.67 cm3/s, ts(sample time), DTC(deadtimecorrectionfactor)=1, td(deadtime)
%
% Outputs: (time stamp in the processed files refer to the start of sampling)
% (first row contain Dp in um)
% raw (1-matdate,2-elapsedtime(s),3:19-Counts,20-Deadtime(s),21-Temp(C),22-RH(%),23-AmbientPress(kPa),24-Alarms)
% conc (1-matdate,2:17-conc(#/cm3),18-Temp(C),19-RH(%),20-AmbientPress(kPa),21-Alarms)
% numsd (1-matdate,2:17-dN/dlogDp in cm-3)
% volsd (1-matdate,2:17-dV/dlogDp in um3/cm3)
% total (1-matdate, 2-totalconc(cm-3), 3-volPM2.5(um3/cm3), 4-volPM10(um3/cm3))
% Contourpots of dV/dlogDp
%
% --------------- Parameters--------------------
caminho='C:\A\Diadema2016\data_trailer_2016\OPS_TSI3330163604\';
caminho1='C:\A\Diadema2016\processados\OPS\';
caminho2='C:\A\Diadema2016\processados\OPS\plots\';
tit='Diadema-2016';
filtmanut=true; % logical para filtrar dados em períodos de manutenção ou com problemas
basicplots=true;
contourplots=false;
lev1=false;
lev2=false;
%
% inicio e termino de periodos de manutencao (horario em utc): removidos das analises
if(filtmanut)
C1=xlsread('C:\A\Diadema2016\logbook_trailer_Diadema2016.xlsx','todos','A4:M135');
C2=xlsread('C:\A\Diadema2016\logbook_trailer_Diadema2016.xlsx','ops','A4:M4');
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));
end
%
% Read raw data files -----------------------------------
raw=zeros(1,24);
files=dir(caminho);
for ifile=3:37%size(files,1)
nome=getfield(files, {ifile,1}, 'name');
disp('-----Reading raw data file')
disp(nome)
fid=fopen(strcat(caminho,nome));
C1=textscan(fid,'%s %s','Delimiter',',');
start=datenum(strcat(C1{1,2}(8),C1{1,2}(7)),'yyyy/mm/ddHH:MM:SS');
N=str2double(C1{1,2}(36));
if (ifile==3)
Dpcut=str2double(C1{1,2}(12:28))';
end
fclose(fid);
fid=fopen(strcat(caminho,nome));
%Elapsed Time [s],Bin 1,Bin 2,Bin 3,Bin 4,Bin 5,Bin 6,Bin 7,Bin 8,Bin 9,Bin 10,Bin 11,Bin 12,Bin 13,Bin 14,Bin 15,Bin 16,Bin 17,Deadtime (s),Temperature (C),Humidity (%),Ambient Pressure (kPa),Alarms,Errors,
C1=textscan(fid,'%f%f%f%f%f%f%f%f%f%f %f%f%f%f%f%f%f%f%f%f %f%f%f%s','Delimiter',',','Headerlines',38);
aux=start+(C1{1,1}/60/60/24)-(C1{1,1}(1)/60/60/24);% time stamp in the end of the sampling period
raw=cat(1,raw,[aux cell2mat(C1(1:23))]);
fclose(fid);
clear C1 aux
end
raw(1,3:19)=Dpcut;
raw=sortrows(raw,1); % sort data based on day of year
save([caminho1,tit,'_OPS_raw_level0.mat'],'raw')
%
% Correct counts and normalize size distributions -----------------------------------
Dpavg=zeros(1,length(Dpcut)-1); dlogDp=zeros(1,length(Dpcut)-1);
for j=2:size(dlogDp,2)+1
Dpavg(1,j-1)=(Dpcut(j)+Dpcut(j-1))/2; % Diameter midpoint
dlogDp(1,j-1)=log10(Dpcut(j)/Dpcut(j-1)); % dlogDp
end
%
% Correct counts for coincidence error based on deadtime
% conc(#/cm3) = N(#)/(Q(cm3/s)*(ts-DTC*td))
% onde Q=16.67 cm3/s, ts(sample time in sec), DTC(deadtimecorrectionfactor)=1, td(deadtime)
conc=zeros(size(raw,1),length(Dpavg));
for j=1:length(Dpavg)
conc(2:length(raw),j)=raw(2:length(raw),j+2)./(16.67*(60-1.0*raw(2:length(raw),20))); % #/cm3
end
conc=[raw(:,1) conc raw(:,21:24)];
conc(1,2:17)=Dpavg;
save([caminho1,tit,'_OPS_conc_level0.mat'],'conc')
%
% Calculate normalized number size distributions from concentrations
numsd=zeros(size(raw,1),length(Dpavg));
for j=1:length(Dpavg)
numsd(2:length(raw),j)=conc(2:length(raw),j+1)/log10(Dpcut(j+1)/Dpcut(j)); % cm-3
end
numsd=[raw(:,1) numsd];
numsd(1,2:17)=Dpavg;
save([caminho1,tit,'_OPS_numsd_level0.mat'],'numsd')
%
% Calculate normalized volume size distributions
volsd=zeros(size(numsd));
for j=1:length(Dpavg)
volsd(2:length(raw),j)=(pi/6)*numsd(2:length(raw),j+1)*(Dpavg(j)^3); % um3/cm3
end
volsd=[raw(:,1) volsd];
volsd(1,2:17)=Dpavg;
save([caminho1,tit,'_OPS_volsd_level0.mat'],'volsd')
%
% Calculate total number concentration, volume 0.25-2.5um (PM2.5), volume 0.25-10um (PM10)
total=zeros(size(numsd,1),4);
total(:,1)=raw(:,1);
for i=2:size(numsd,1)
total(i,2)=numsd(i,2:17)*dlogDp(1,:)'; % total conc cm-3
total(i,3)=volsd(i,2:11)*dlogDp(1,1:10)'; % vol PM2.5 um3/cm3
total(i,4)=volsd(i,2:17)*dlogDp(1,1:16)'; % vol PM10 um3/cm3
end
total(1,:)=[];
save([caminho1,tit,'_OPS_totconc_level0.mat'],'total')
%
% Basic Plots -----------------------------------
if(basicplots)
aux=2:length(raw);
figure
subplot(2,2,1),plot(raw(aux,1),raw(aux,21),'r.'),grid on,ylabel('Temp(C)'),dynamicDateTicks
title(tit)
subplot(2,2,2),plot(raw(aux,1),raw(aux,22),'c.'),grid on,ylabel('RH(%)'),dynamicDateTicks
subplot(2,2,3),plot(raw(aux,1),raw(aux,23),'g.'),grid on,ylabel('Press(kPa)'),dynamicDateTicks
subplot(2,2,4),plot(raw(aux,1),raw(aux,24),'k*'),grid on,ylabel('Alarms'),dynamicDateTicks
saveas(gcf,[caminho2,tit,'_OPS_param_level0'],'jpg')
%
figure,plot(raw(aux,1),raw(aux,3),'r.'),grid on,ylabel('Counts Bin1'),dynamicDateTicks
%
figure
subplot(3,1,1),plot(total(:,1),total(:,2),'b.'),grid on,ylabel('TotalConc (cm^-^3)'),dynamicDateTicks
title(tit)
subplot(3,1,2),plot(total(:,1),total(:,3),'r.'),grid on,ylabel('Vol PM2.5 (um^3/cm^-^3)'),dynamicDateTicks
subplot(3,1,3),plot(total(:,1),total(:,4),'k.'),grid on,ylabel('Vol PM10 (um^3/cm^-^3)'),dynamicDateTicks
saveas(gcf,[caminho2,tit,'_OPS_totalconc_level0'],'jpg')
%
figure,plot(mod(total(:,1),1),total(:,2),'b.'),grid on
ylabel('TotalConc (cm^-^3)'),xlabel('UTC'),dynamicDateTicks,title(tit)
saveas(gcf,[caminho2,tit,'_OPS_totalconc_dcycle_level0'],'jpg')
end
%
% -------------------------------------------------------------------------
% Level 1
% Keep records when: 1< TotalConc< 1000 cm-3
% Remove mantainance periods according to the field logbook.
if(lev1)
C=load([caminho1,'Diadema-2016_OPS_raw_level0.mat']);
raw=C.raw;
C=load([caminho1,'Diadema-2016_OPS_conc_level0.mat']);
conc=C.conc;
C=load([caminho1,'Diadema-2016_OPS_numsd_level0.mat']);
numsd=C.numsd;
C=load([caminho1,'Diadema-2016_OPS_volsd_level0.mat']);
volsd=C.volsd;
C=load([caminho1,'Diadema-2016_OPS_totconc_level0.mat']);
total=C.total;
total_level0=total;
diam=raw(1,3:18);
tit='Diadema-2016';
origsize=size(raw,1)-1;
%
% Remove maintanance periods from dataset
for i=1:size(manut,1)
aux=find(raw(:,1)>=manut(i,1) & raw(:,1)<=manut(i,2));
raw(aux,:)=[];conc(aux,:)=[];
numsd(aux,:)=[];volsd(aux,:)=[];
aux1=find(total(:,1)>=manut(i,1) & total(:,1)<=manut(i,2));
total(aux1,:)=[];
clear aux aux1
end
%
aux=find(total(:,2)<1 | total(:,2)>1000); % remove lines with totconc<1cm-3 or >1000cm-3
total(aux,:)=[];
raw(aux+1,:)=[];conc(aux+1,:)=[];
numsd(aux+1,:)=[];volsd(aux+1,:)=[];
%
figure,plot(total_level0(:,1),total_level0(:,2),'r.'),grid on, hold on
plot(total(:,1),total(:,2),'b.')
dynamicDateTicks
%
% Save Level 1 data
Instrument='OPS_TSI';
Model='3330';
serial_num='3330163604';
Site=tit;
Val_level='1';
Proc_vers='1.0';
labels{1}='Total concentration';
labels{2}='PM2.5 Volume';
labels{3}='PM10 Volume';
units{1}='(cm-3)';
units{2}='(um3.cm-3)';
units{3}='(um3.cm-3)';
aux=2:length(numsd);
tempo=numsd(aux,1);
var=total(:,2:4);
sizedist=numsd(aux,2:size(numsd,2));
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,tempo,var,diam,sizedist);
end
%
% Contourplots -------------------------------------------------------
if(contourplots)
disp('Calculations for contourplots')
diam=volsd(1,2:17);
volsd(1,:)=[];
aux=find(diff(volsd(:,1))>=0.5);
volsd2(1,:)=volsd(1,:);
for i=2:length(aux) % insere linhas em branco em períodos com mais de 12h sem dados
volsd2=cat(1,volsd2,volsd(aux(i-1)+1:aux(i),:),nan(1,size(volsd,2)));
end
volsd2=cat(1,volsd2,volsd(aux(length(aux))+1:size(volsd,1),:));
aux=find(isnan(volsd2(:,1)));
volsd2(aux,1)=volsd2(aux-1,1)+0.4;
timesec=volsd2(:,1);
Z=volsd2(:,2:17);
[r,c]=find(Z<0); % replace negative values per 0
for k=1:length(r)
Z(r(k),c(k))=0;
end
clear r c
Znorm=Z;
for ii=1:size(Z,1)
Znorm(ii,:)=Z(ii,:)./max(Z(ii,:));
end
%
figure
contourf(timesec,diam,Z','LineStyle','none'),grid on
set(gca,'YScale','log')
set(gca,'YTick',[0.5 1 2.5 5 7.5])
colorbar
title([tit,'-OPS Vol Size Dist'],'FontSize',15),xlabel('UTC Time','FontSize',12),ylabel('Diameter (um)','FontSize',12)
set(gca,'XLim',[datenum(2016,11,22) datenum(2017,3,16)])
set(gca,'XTick',datenum(2016,11,22):15:datenum(2017,3,16))
set(gca,'XTickLabel', datestr(datenum(2016,11,22):15:datenum(2017,3,16),'mm/dd/yy'),'FontSize',12)
saveas(gcf,strcat(caminho2,tit,'_contour_tseries'),'jpg')
%
figure
contourf(timesec,diam,Znorm','LineStyle','none'),grid on
set(gca,'YScale','log')
set(gca,'YTick',[0.5 1 2.5 5 7.5])
legend('normalized')%,colorbar
title([tit,'-OPS Vol Size Dist'],'FontSize',15),xlabel('UTC Time','FontSize',12),ylabel('Diameter (um)','FontSize',12)
set(gca,'XLim',[datenum(2016,11,22) datenum(2017,3,16)])
set(gca,'XTick',datenum(2016,11,22):15:datenum(2017,3,16))
set(gca,'XTickLabel', datestr(datenum(2016,11,22):15:datenum(2017,3,16),'mm/dd/yy'),'FontSize',12)
saveas(gcf,strcat(caminho2,tit,'_contour_tseries_norm'),'jpg')
%
numsd(1,:)=[];
mediana_num=zeros(3,length(2:17));
mediana_vol=zeros(3,length(2:17));
for ii=2:17
mediana_num(1,ii-1)=prctile(numsd(:,ii),25);
mediana_num(2,ii-1)=prctile(numsd(:,ii),50);
mediana_num(3,ii-1)=prctile(numsd(:,ii),75);
mediana_vol(1,ii-1)=prctile(volsd(:,ii),25);
mediana_vol(2,ii-1)=prctile(volsd(:,ii),50);
mediana_vol(3,ii-1)=prctile(volsd(:,ii),75);
end
figure
subplot(1,2,1)
errorbar(diam,mediana_num(2,:),mediana_num(1,:),mediana_num(3,:),'-b.')
title(strcat(tit,'-Median Number Size Dist'),'FontSize',12,'fontweight','b'),xlabel('Dp (um)','FontSize',12,'fontweight','b'),ylabel('dN/dlogDp (cm^-^3)','FontSize',12,'fontweight','b')
set(gca,'XScale','log'), grid on
set(gca,'XTick',[0.5 1 2.5 5 7.5],'FontSize',12,'fontweight','b')
%set(gca,'XLim',[10 500])
subplot(1,2,2)
errorbar(diam,mediana_vol(2,:),mediana_vol(1,:),mediana_vol(3,:),'-r.')
title(strcat(tit,'-Median Volume Size Dist'),'FontSize',12,'fontweight','b'),xlabel('Dp (um)','FontSize',12,'fontweight','b'),ylabel('dV/dlogDp (um^3/cm^-^3)','FontSize',12,'fontweight','b')
set(gca,'XScale','log'), grid on
set(gca,'XTick',[0.5 1 2.5 5 7.5],'FontSize',12,'fontweight','b')
saveas(gcf,strcat(caminho2,tit,'_OPS_median_sizedist'),'jpg')
end
%
%
% ------------ Level 2 ------------------------------------------------
% Level2: STP correction (1013.25; 293.15).
if(lev2)
% Lê os arquivos Level 1 que ja salvou
aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_OPS_TSI_3330_Level1_5min.csv');
diam=zeros(1,16);
for j=24:39
diam(j-23)=str2num(char(aux.textdata(2,j)));
end
level1_5min=[aux.data(:,1:2) aux.data(:,9) aux.data(:,16) aux.data(:,23:38)];
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';
%
% Lê dados de meteorologia e calcula fator stp
C=importdata('C:\A\Diadema2016\meteorologia_IAG\Diadema-2016_Meteo_5min.csv');
met(:,1)=datenum(C.textdata(3:length(C.textdata)),'dd-mmm-yyyy HH:MM:SS ');
met(:,2:3)=C.data;
met(:,1)=met(:,1)+3/24; % local time to utc time
met(:,3)=met(:,3)+273.15; % oC to K
met(met(:,1)<level1_5min(1,1),:)=[];
met(met(:,1)>level1_5min(length(level1_5min),1),:)=[];
if(length(met)~=length(level1_5min))
disp('check length of meteorological data file')
stop
end
%
pref=1013.25; tref=293.15;
stp=(pref*met(:,3))./(tref*met(:,2));
for j=2:20
level1_5min(:,j)=level1_5min(:,j).*stp;
end
%
% Save Level 2 data
Instrument='OPS_TSI';
Model='3330';
serial_num='3330163604';
Site=tit;
Val_level='2';
Proc_vers='1.0';
labels{1}='Total concentration';
labels{2}='PM2.5 Volume';
labels{3}='PM10 Volume';
units{1}='(cm-3)';
units{2}='(um3.cm-3)';
units{3}='(um3.cm-3)';
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level1_5min(:,1),level1_5min(:,2:4),diam,level1_5min(:,5:20));
end
%
%