% Human-in-the-loop QC procedure to generate final QC'd OCSD data set.

% All the 'processing' (applying cal's, deriving of engineering units, etc.)
%  is now done without human input in python scripts, either in
%  real-time (makeOAwebQC_OA1.py), or post-cruise (loadOAdata_OA1.py).
% This Matlab script is used only to identify and remove bad data.

% Outputs final .csv files in the /QC directory
% A final .nc file (merged to a common 10-min time stamp) needs to be
%  created after this script is run.

% Updated on 03/21/19 by Devon Northcott to add in the ability to remove
% NaN's that were incorectly applied, and to streamline most of the
% operations into sub-funtions, hopefully making the script more modular,
% and making it easier to build QC scripts for other platforms

%Adapted OA code to OCSD 6/2020, by Diego Sancho
%Needs subfunctions of folder QCCODE in Diatom/users/dsancho/MATLAB/QCCode

%depl info
 veh = 'OCSD';
 dep = '20190918';  

%directories for post-cruise and QC'd data
 atlas = '//atlas.shore.mbari.org/';
 mydirPC = [atlas 'Diatom/socal moorings/deployment_data/' dep '/realTime/'];                       
 mydir = [atlas 'Diatom/socal moorings/deployment_data/' dep]; 
 mydirpH= [atlas 'Diatom/socal moorings/deployment_data/' dep '/Durafet pH recorder data/']; 
     
 %add mapping and qc code
 addpath(genpath([atlas 'Diatom/users/dsancho/MATLAB/utility'])) 
 addpath(genpath([atlas 'Diatom/users/dsancho/MATLAB/QCCode']))
 clear veh dep;
 %%
 %Figure info
 set(groot,'defaultFigureCreateFcn',@(fig,~)addToolbarExplorationButtons(fig))
 set(groot,'defaultAxesCreateFcn',@(ax,~)set(ax.Toolbar,'Visible','off'))
 
 %data often include more than just one full deployment, keep dates of
 %interest only
 tbeg=datenum(2019,09,18,19,30,00);
 tend=datenum(2019,12,02,22,00,00);
% read-in postCruise .csv files.  NaN zeros, convert time stamp
% .csv files must be in the format indicated otherwise the QC functions 
% will fail.
 dn = datenum([1970 1 1]);
 %Position [time,lon,lat]
  Position = Load_CSV_for_QC([mydirPC 'Position.csv']);
  f=find(Position(:,1)<tbeg | Position(:,1) > tend);
  Position(f,:)=[];
 %CO2 [time, pCO2water, pCO2air]
  CO2 = Load_CSV_for_QC([mydirPC 'CO2.csv']);
  f=find(CO2(:,1)<tbeg | CO2(:,1) > tend);
  CO2(f,:)=[];
 %CO2equil [time, temp-ON, temp-OFF, Press-ON, Press-OFF, xCO2-ON,
 %xCO2-OFF]
  CO2equil=Load_CSV_for_QC([mydirPC 'CO2equil.csv']);
  f=find(CO2equil(:,1)<tbeg | CO2equil(:,1) > tend);
  CO2equil(f,:)=[];
  %CO2stnd [time, temp-ON, temp-OFF, Press-ON, Press-OFF, xCO2-ON,
 %xCO2-OFF]
  CO2stnd=Load_CSV_for_QC([mydirPC 'CO2stnd.csv']);
  f=find(CO2stnd(:,1)<tbeg | CO2stnd(:,1) > tend);
  CO2stnd(f,:)=[];
  %CO2zero [time, temp-ON, temp-OFF, Press-ON, Press-OFF, xCO2-ON,
 %xCO2-OFF]
  CO2zero=Load_CSV_for_QC([mydirPC 'CO2zero.csv']);
  f=find(CO2zero(:,1)<tbeg | CO2zero(:,1) > tend);
  CO2zero(f,:)=[];
  
  %CO2air [time, temp-ON, temp-OFF, Press-ON, Press-OFF, xCO2-ON,
 %xCO2-OFF]
  CO2air=Load_CSV_for_QC([mydirPC 'CO2air.csv']);
 f=find(CO2air(:,1)<tbeg | CO2air(:,1) > tend);
  CO2air(f,:)=[];
 %Temp [time,T 5m, T 35m, T 60m]
  Temp = Load_CSV_for_QC([mydirPC 'Temperature.csv']);
  f=find(Temp(:,1)<tbeg | Temp(:,1) > tend);
  Temp(f,:)=[];
 %Density [time, Dens 5m, Dens 35m, Dens 60m]
  Dens = Load_CSV_for_QC([mydirPC 'Density.csv']);
  f=find(Dens(:,1)<tbeg | Dens(:,1) > tend);
  Dens(f,:)=[];
 %Fluor [time,fluor 35m]
  Fluor = Load_CSV_for_QC([mydirPC 'Fluor.csv']);
  f=find(Fluor(:,1)<tbeg | Fluor(:,1) > tend);
  Fluor(f,:)=[];
 %MET [time,airPress, airTemp]
  MET = Load_CSV_for_QC([mydirPC 'MET.csv']);
  f=find(MET(:,1)<tbeg | MET(:,1) > tend);
  MET(f,:)=[];
 %MET wind [time,windDir, windSpd]
  METw = Load_CSV_for_QC([mydirPC 'METwind.csv']);
  f=find(METw(:,1)<tbeg | METw(:,1) > tend);
  METw(f,:)=[];
 %Salinity [time, Sal 5m, Sal 35m, Sal 60m]
  Sal = Load_CSV_for_QC([mydirPC 'Salinity.csv']);
  f=find(Sal(:,1)<tbeg | Sal(:,1) > tend);
  Sal(f,:)=[];
 %O2 [time,O2 5m, O2 35m, O2 60m]
  O2=Load_CSV_for_QC([mydirPC 'Oxygen.csv']);
  f=find(O2(:,1)<tbeg | O2(:,1) > tend);
  O2(f,:)=[];
  
  Press=Load_CSV_for_QC([mydirPC 'Pressure.csv']);
  f=find(Press(:,1)<tbeg | Press(:,1) > tend);
  Press(f,:)=[];
  
  %pH --processed by convOCSD_internal to convert from volts to pH
  CTD_ph=load([mydirpH 'CTD_pH_conv.mat']);
  %%
 f=find(CTD_ph.pH79int(:,1)>=tbeg & CTD_ph.pH79int(:,1)<=tend);
 ph1=CTD_ph.pH79int(f,3);
  f=find(CTD_ph.pH80int(:,1)>=tbeg & CTD_ph.pH80int(:,1)<=tend);
 ph2=CTD_ph.pH80int(f,3);
 ph=[CTD_ph.pH80int(f,1) ph1 ph2];
  
  %%
 CO2hum=nan(size(CO2air));
%Plot all sensor data in one fig
% This is to compare different data sets for simultaneous 'events' -
%  i.e. if there is a spike in O2 and pCO2w at the same time, they 
%  are both likely to be 'real.' This is to aid removal of data which
%  is implemented in the individual sensor break-out figures.
 qcOCSD_plotAll(Position,CO2,ph,Temp,Sal,Press,Dens,O2,Fluor,MET,METw);

%%
%m_map of position
m_proj('Mercator','lon',([-128 -116]),'lat',([32 40]))
m_scatter(-Position(:,2),Position(:,3),15,'filled')
hold on
m_gshhs_f('patch',[.5 .5 .5],'edgecolor','k');

 %%
 %Plot CO2 sys diagnostics
  ha=tight_subplot(7,1,0.05,0.1,0.1);
 g = get(ha,'Position');
 set(ha,'Box','on','NextPlot','add');
 axes(ha(1)); %press
   plot(CO2equil(:,1),CO2equil(:,5),'.-','Color',[.8 0 0]);
   plot(CO2zero(:,1),CO2zero(:,5),'.-','Color',[0 0 0]);
   plot(CO2air(:,1),CO2air(:,5),'.-','Color',[0 0 .8]);
   f = find(~isnan(CO2stnd(:,7)));
   plot(CO2stnd(f,1),CO2stnd(f,5),'.-','Color',[0 .8 0]);
   ylabel('[kPa]');
   title('Pressure');
 axes(ha(2)); %temp
   plot(CO2equil(:,1),CO2equil(:,3),'.-','Color',[.8 0 0]);
   plot(CO2zero(:,1),CO2zero(:,3),'.-','Color',[0 0 0]);
   plot(CO2air(:,1),CO2air(:,3),'.-','Color',[0 0 .8]);
   plot(CO2stnd(f,1),CO2stnd(f,3),'.-','Color',[0 .8 0]);
   ylabel('[C]');
   title('Temperature');
 axes(ha(3)); %humidity
   plot(CO2hum(:,1),CO2hum(:,2),'.-','Color',[.8 0 0]);
   plot(CO2hum(:,1),CO2hum(:,3),'.-','Color',[0 0 .8]);
   plot(CO2hum(f,1),CO2hum(f,4),'.-','Color',[0 .8 0]);
   ylabel('[V]');
   title('Humidity');
 axes(ha(4)); %xCO2
   plot(CO2equil(:,1),CO2equil(:,7),'.-','Color',[.8 0 0]);
   plot(CO2zero(:,1),CO2zero(:,7),'.-','Color',[0 0 0]);
   plot(CO2air(:,1),CO2air(:,7),'.-','Color',[0 0 .8]);
   plot(CO2stnd(f,1),CO2stnd(f,7),'.-','Color',[0 .8 0]);
   ylabel('[ppm]');
   title('xCO2');
   
 axes(ha(5)); %pCO2w
   f = find(~isnan(CO2(:,3)));
   plot(CO2(f,1),CO2(f,3),'.-');
   ylabel('[ppm]');
   title('pCO2w');

 axes(ha(6)); %pCO2a
   f = find(~isnan(CO2(:,4)));
   plot(CO2(f,1),CO2(f,4),'.-');
   ylabel('[ppm]');
   title('pCO2a');
  
 axes(ha(7)); %pH
  %pH
   plot(ph(:,1),ph(:,2),'.-');
   hold on
   plot(ph(:,1),ph(:,3),'.-');
   if size(ph,2)>3
       hold on
       plot(ph(:,1),ph(:,4),'.-');
       legend('5m','35m','60m');
   end
   legend('35m','60m')
   ylabel('pH');
   title('pH')
   
 set(ha,'XLim',[min(CO2(:,1)) max(CO2(:,1))] + [-1 1]*0.1);
 set(ha(1:end-1),'XTickLabel',[]);
 linkaxes(ha,'x');
 datetick('x','mm/dd','keeplimits')
 %%
 %manually nan bad data 
 CO2=nan; %all CO2 data bad for this deployment
 %nan out temperature from Wx
 MET(:,3)=nan; 
 %remove one spike in wind speed
 METw(METw(:,3)>40)=nan;
%%

%Final Check:
 qcOCSD_plotAll(Position,CO2,ph,Temp,Sal,Press,Dens,O2,Fluor,MET,METw);

 %%
 %save QC'd matfile
 cd(mydir)
 save('20190918_QC.mat','CO2','CO2air','CO2stnd','CO2equil','CO2zero','Fluor',...
     'MET','METw','O2','ph','Sal','Temp','Press','Dens','Position'); 
