% MAAP data processing - Diadema 2016
% Read raw data files and use LFA criteria to generate Level1, 2 and 3 datasets
% Functions used: frepeated; export_LFA_GoAmazon_jan2015
% Outputs: Level0 (mat_time(utc); bc(ug/m3); bc_mass; temp(oC); ?; press(mbar); ?; ?
% Level1_5min (dd-mmm-yyyy HH:MM; sec_since_1970; mean; std; median; p25; p75; max; min; num)
% Level1_30min
% Level 1: # of iterations must be <50 and >0; 1-min BC data must be >-1 ug.m-3 and <100 ug.m-3
% Level 2: Apply 5% correction (Müller et al.,2011). Apply STP corrections.
% Level 3: MAAP BC must be smaller than 50% of PM2.5 measured from TEOM (Thermo).
% If no TEOM data is available, use SMPS integrated mass and average density of 1.4 g cm-3
%
% ---------------- Parameters -----------------------
caminho1='C:\A\Diadema2016\flex_data_files\MAAP156\csv\'; % arquivos flex (Labview)
caminho2='C:\A\Diadema2016\data_trailer_2016\maap\'; % arquivos do datalogger do instrumento
caminho3='C:\A\Diadema2016\processados\maap\';
tit='Diadema-2016';
pref=1013.25; % Reference pressure for STP correction (mbar)
tref=273.15; % Reference temp for STP correction (K)
lerflex=true; % logical para ler arquivos flex (Labview)
lerdatalogger=true; % logical para ler arquivos do datalogger do instrumento
corrigirtempoflex=true;
corrigirtempodatalogger=true;
filtmanut=true; % logical para filtrar dados em períodos de manutenção ou com problemas
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','neph_maap','A4:M5');
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
%
% ---------------- Reading raw files -----------------------
if(lerdatalogger)
data1=zeros(1,2);
files=dir(caminho2);
disp('-----Reading raw data files')
for ifile=3:size(files,1)
nome=getfield(files, {ifile,1}, 'name');
disp(nome)
fid=fopen(strcat(caminho2,nome));
C=textscan(fid,'%s%s%*s%f');
fclose(fid);
aux(:,4)=C{1,3};
aux(:,1)=datenum(C{1,1},'yy-mm-dd');
aux(:,2)=mod(datenum(C{1,2},'HH:MM'),1);
aux(:,3)=aux(:,1)+aux(:,2);
data1=cat(1,data1,aux(:,3:4)); % matdate,bc(ug/m3)
clear aux C nome icol
end
data1(1,:)=[];
data1=[data1(:,1) nan(length(data1),5) data1(:,2) nan(length(data1),1)];
else
data1=[];
end
%
if(lerflex)
files=dir(caminho1);
data2=zeros(1,8);
disp('-----Reading flexdata files')
for ifile=3:size(files,1)
nome=getfield(files, {ifile,1}, 'name');
disp(nome)
fid=fopen(strcat(caminho1,nome));
C=cell(floor(getfield(files, {ifile,1}, 'bytes')/330),10);
i=1;
while (~feof(fid))
C(i,:)=textscan(fid,'%s%s%s%*s%*s%*s %*f%*f%*f%*f%*f%*f%*f %f%f%f%f%f%f%f',1,'Delimiter',',','MultipleDelimsAsOne',1);
%C(i,:)=textscan(fid,'%*s%s%s%s%*s%*s%*s %*f%*f%*f%*f%*f%*f%*f%*f %f%f%f%f%f%f%f%f',1,'Delimiter',',');
i=i+1;
end
fclose(fid);
%
% removendo linhas com defeito
rem=0;
for j=4:10
aux=find(cellfun(@isempty,C(:,j)));
rem=cat(1,rem,aux);
end
rem(1)=[];
C(rem,:)=[];
%
% removendo mais linhas com defeito
aux1=find(cellfun(@isnan,C(:,5)));
aux2=find(cellfun(@isnan,C(:,10)));
rem=cat(1,aux1,aux2);
rem=unique(rem);
C(rem,:)=[];
%
data_(:,2:8)=cell2mat(C(:,4:10));
for i=1:length(C)
try
data_(i,1)=datenum(strcat(C{i,1},C{i,2},C{i,3}),'mm/dd/yyyyHH:MM:SSAM');
catch
data_(i,1)=NaN;
end
end
data_(find(isnan(data_(:,1))),:)=[];
%
% Um erro que acontece às vezes, é aparecer o mês na última coluna
% de dados. Removendo estes casos aqui.
[~,mes,~,~,~,~]=datevec(data_(1,1));
data_(data_(:,8)==mes,8)=NaN;
data2=cat(1,data2,data_);
clear data_ C mes
end
data2(1,:)=[];
else
data2=[];
end
%
if(corrigirtempodatalogger) % tapando buracos do flexdata (períodos com computador travado, dados armazenados no datalogger do instrumento)
aux1=find(data1(:,1)>=datenum(2016,11,26,18,31,00) & data1(:,1)<=datenum(2016,11,28,20,03,00));
aux2=find(data1(:,1)>=datenum(2016,12,03,22,42,00) & data1(:,1)<=datenum(2016,12,05,19,03,00));
aux3=find(data1(:,1)>=datenum(2016,12,10,21,18,00) & data1(:,1)<=datenum(2016,12,11,20,00,00));
aux4=find(data1(:,1)>=datenum(2016,12,20,20,05,00) & data1(:,1)<=datenum(2016,12,22,14,15,00));
aux5=find(data1(:,1)>=datenum(2016,12,31,19,07,00) & data1(:,1)<=datenum(2017,01,01,19,08,00));
aux6=find(data1(:,1)>=datenum(2017,01,02,19,40,00) & data1(:,1)<=datenum(2017,01,05,16,26,00));
aux7=find(data1(:,1)>=datenum(2017,02,24,20,15,00) & data1(:,1)<=datenum(2017,02,25,21,15,00));
aux8=find(data1(:,1)>=datenum(2017,02,27,21,35,00) & data1(:,1)<=datenum(2017,02,28,15,05,00));
aux9=find(data1(:,1)>=datenum(2017,03,02,21,45,00) & data1(:,1)<=datenum(2017,03,02,22,10,00));
aux10=find(data1(:,1)>=datenum(2017,03,04,17,15,00) & data1(:,1)<=datenum(2017,03,06,11,59,00));
aux11=find(data1(:,1)>=datenum(2017,03,06,14,29,00) & data1(:,1)<=datenum(2017,03,07,16,52,00));
aux12=find(data1(:,1)>=datenum(2017,03,07,23,44,00) & data1(:,1)<=datenum(2017,03,09,21,43,00));
aux=data1; clear data1
data1=aux(cat(1,aux1,aux2,aux3,aux4,aux5,aux6,aux7,aux8,aux9,aux10,aux11,aux12),:);
clear aux1 aux2 aux
end
if(corrigirtempoflex) % correção segundo o logbook
%subtrair 1h dos horarios anteriores a 25nov2016 16:25
aux1=find(data2(:,1)==datenum(2016,11,25,17,24,36));
aux2=find(data2(:,1)<=datenum(2016,11,25,17,24,06));
if(length(aux1)~=1)
disp('Problema na correcao do tempo')
stop
end
data2(aux2(1):aux1,1)=data2(aux2(1):aux1,1)-datenum(0,0,0,1,0,0);
%somar 2h dos horarios entre 25nov2016 14:25 e 15:19:33
aux1=find(data2(:,1)==datenum(2016,11,25,14,25,10));
aux2=find(data2(:,1)==datenum(2016,11,25,15,19,06));
if(length(aux1)~=1 || length(aux2)~=1)
disp('Problema na correcao do tempo')
stop
end
data2(aux1:1:aux2,1)=data2(aux1:1:aux2,1)+datenum(0,0,0,2,0,0);
end
%
data=cat(1,data1,data2);
aux=sortrows(data,1); % sort data based on day of year
data=aux;
%
% Obs: é comum a repetição de horários nos arquivos flexdata. Nada grave,
% dura alguns segundos por dia, com valores de espalhamento ligeiramente diferentes.
% Tem a ver com o fato de que o LabView loga os dados a cada segundo,
% enquanto que o neph mede a cada 5 segundos. Assim, vou dispensar a verificação de
% dados repetidos aqui, e remover linhas com mesmo time stamp.
disp('----Checking for repeated data')
aux=find(diff(data(:,1))==0);
disp(['Número de linhas com timestamp repetido: ',num2str(length(aux))])
data(aux,:)=[];
clear aux
% [newdata,erro]=frepeated(data,2);
% if (strcmp(erro,'none')~=1)
% disp(erro)
% else
% data=newdata;
% end
% clear aux erro
%
save([caminho3,tit,'_maap_level0.mat'],'data','-mat')
clear data1 data2
%
% Remove períodos de manutencao/contaminacao de acordo com o logbook
if(filtmanut)
dataorig=data;
for i=1:size(manut,1)
aux=find(data(:,1)>=manut(i,1) & data(:,1)<=manut(i,2));
if(manut(i,3)==0)
data(aux,:)=[]; % remove maintanance periods from dataset
end
end
save([caminho3,tit,'_maap_clean.mat'],'data','-mat')
figure,plot(dataorig(:,1),dataorig(:,7),'r.'),hold on,grid on
plot(data(:,1),data(:,7),'b.')
ylabel('BC (ng/m^3)'),title(tit)
dynamicDateTicks
saveas(gcf,strcat(caminho3,tit,'_maap_level0'),'jpg')
else
figure
plot(data(:,1),data(:,7),'k.'),grid on
ylabel('BC (ng/m^3)'),title('tit')
dynamicDateTicks
saveas(gcf,strcat(caminho3,tit,'_maap_level0'),'jpg')
end
%
%
if(lev1)
% ------ Level1 ---------------------------
% Level 1: # of iterations must be <50 and >0;
% 1-min BC data must be >-5 ug.m-3 and <100 ug.m-3
C=load('C:\A\Diadema2016\processados\maap\Diadema-2016_maap_clean.mat');
data=C.data;
C=load('C:\A\Diadema2016\processados\maap\Diadema-2016_maap_level0.mat');
data_level0=C.data;
tit='Diadema-2016';
%
origsize=size(data_level0,1);
data(data(:,7)<-5,:)=[];
data(data(:,7)>100,:)=[];
data(data(:,3)>50 | data(:,3)<0,:)=[];
%
figure,plot(data_level0(:,1),data_level0(:,7),'r.'),grid on,hold on
plot(data(:,1),data(:,7),'k.')
dynamicDateTicks
%
% Save Level 1 data ----------
%data(:,7)=data(:,7)*1000; % ug/m3 to ng/m3
Instrument='MAAP';
Model='5012';
serial_num='516';
Site=tit;
Val_level='1';
Proc_vers='1';
labels={'BC'};
units={'(ug.m-3)'};
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,data(:,1),data(:,7))
end
%
if(lev2)
% ------ Level2 ---------------------------
% Level 2: Apply 5% correction (Müller et al.,2011). Apply STP corrections.
% Lê os arquivos Level 1 que ja salvou
aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_MAAP_5012_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
origsize=sum(~isnan(level1_5min(:,2)));
tit='Diadema-2016';
clear aux
%
% 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));
%
% STP and 5% correction (Muller et al., 2011)
level1_5min(:,2)=level1_5min(:,2).*stp*1.05;
%
% Save Level 2 data ----------------
level1_5min(:,3)=6.6*level1_5min(:,2);
Instrument='MAAP';
Model='5012';
serial_num='516';
Site=tit;
Val_level='2';
Proc_vers='1';
labels={'BC637';'Abs637'};
units={'(ug.m-3)';'(Mm-1)'};
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level1_5min(:,1),level1_5min(:,2:3))
%
end

% ------ Level3 ---------------------------
% Level 3: MAAP BC must be smaller than 50% of PM2.5 measured from TEOM (Thermo).
% If no TEOM data is available, use SMPS integrated mass and average density of 1.4 g cm-3.

%
% Verificação final
% aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_MAAP_5012_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_MAAP_5012_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.')
% figure,plot(level1_5min(:,2),level2_5min(:,2),'b.'),grid on