% -- 
% Jeff Sevadjian
% Monterey Bay Aquarium Research Institute
% 7700 Sandholdt Rd
% Moss Landing, CA 95039, USA
% ph: +1 831-775-1888

% Loads real-time .csv files output from makeOAwebQC_oa1.py,
%  does a fairly thorough manual removal of bad data,
%  saves as a .mat file for Diego to extend the analyses.

% This script was run during the 'middle' of the deployment and is not
%  intended as a 'final' QC'd .mat file.

% Output matrices are same as the older QC'd OA1, OA2, M1 data, e.g.:

%  CTD:   Time [Matlab serial datenum], Temperature [C], Salinity 
%  Fluor: Time, fluorescence [ug/l] 
%  MET:   Time, barometric pressure [bar], raw (uncorrected) air temp [C],
%         relative humidity [%], corrected air temp [C], wind direction [T],
%         wind speed [m/s]
%  O2:    Time, oxygen [umol/kg] 
%  pCO2:  Time, pCO2 in water, pCO2 in air [ppm] 
%  pH:    total scale 
%  VR2C:  Time, tag ID

atlas = '//atlas.shore.mbari.org/';
buoy = 'oa1';
dep = '201607';
mydir = [atlas 'OA_Moorings/deployment_data/' buoy '/' dep '/realTime/'];

dn = datenum([1970 1 1]);

%--%

%CTD
 CTD = csvread([mydir 'CTD.csv'],1,0);
 CTD(:,1) = dn + (CTD(:,1)/86400/1000);

%get rid of non-unique time stamps (?)
 [~,i] = unique(CTD(:,1));
 CTD = CTD(i,:);
 clear i;

%QC temp
 figure; hold on;
 plot(CTD(:,1),CTD(:,2),'.-');
 %no manual outliers

%QC salinity
 figure; hold on;
 plot(CTD(:,1),CTD(:,3),'.-');
 %fixed floor
  f = find(CTD(:,3) < 15);
  plot(CTD(f,1),CTD(f,3),'m.');
  CTD(f,3) = NaN;
  clear f;
 %variance filter
  win = 100; %window
  nst = 5; %number of stdev's above median
  sm = slidingmed(CTD(:,3), win);
  ss = slidingstdev(CTD(:,3), win);
  fs = find(abs(CTD(:,3) - sm) > nst*ss);
  plot(CTD(fs,1),CTD(fs,3),'r.');
  CTD(fs,3) = NaN;
  clear win nst sm ss fs;
 %manual outliers
  f = find( CTD(:,1)>datenum([2017 3 22 16 55 0]) & ...
    CTD(:,1)<datenum([2017 3 22 19 0 0]) );
  plot(CTD(f,1),CTD(f,3),'k.');
  CTD(f,2:end) = NaN;
  clear f;

%get rid of 4th column (pressure)
 CTD(:,4) = [];

%--%

%O2
 O2 = csvread([mydir 'Oxygen.csv'],1,0);
 O2(:,1) = dn + (O2(:,1)/86400/1000);

%get rid of non-unique time stamps (?)
 [~,i] = unique(O2(:,1));
 O2 = O2(i,:);
 clear i;

%QC O2
 figure; hold on;
 plot(O2(:,1),O2(:,2),'.-');
 %missing values
  f = find(O2(:,2) == 0);
  plot(O2(f,1),O2(f,2),'m.');
  O2(f,2) = NaN;
  clear f;
 %manual outliers
  f = find( O2(:,1)>datenum([2017 3 22 16 55 0]) & ...
    O2(:,1)<datenum([2017 3 22 19 0 0]) );
  plot(O2(f,1),O2(f,2),'k.');
  O2(f,2) = NaN;
  clear f;

%--%

%Fluor
 Fluor = csvread([mydir 'Fluor.csv'],1,0);
 Fluor(:,1) = dn + (Fluor(:,1)/86400/1000);

%get rid of non-unique time stamps (?)
 [~,i] = unique(Fluor(:,1));
 Fluor = Fluor(i,:);
 clear i;

%QC Fluor
 figure; hold on;
 plot(Fluor(:,1),Fluor(:,2),'.-');
 %no missing values
 %manual outliers
  f = find( Fluor(:,1)>datenum([2017 3 22 16 55 0]) & ...
    Fluor(:,1)<datenum([2017 3 22 19 0 0]) );
  plot(Fluor(f,1),Fluor(f,2),'k.');
  Fluor(f,2) = NaN;
  clear f;

%--%

%pH, pCO2
 data = csvread([mydir 'CO2.csv'],1,0);

%get rid of non-unique time stamps (?)
 [~,i] = unique(data(:,1));
 data = data(i,:);
 clear i;

%--%

%QC pH
 pH = [dn + (data(:,1)/86400/1000) data(:,2)];
 figure; hold on;
 plot(pH(:,1),pH(:,2),'.-');
 %missing values
  f = find(pH(:,2) == 0);
  plot(pH(f,1),pH(f,2),'m.');
  pH(f,2) = NaN;
  clear f;
 %manual outliers
  f = find( pH(:,1)>datenum([2017 3 18 0 37 0]) & ...
    pH(:,1)<datenum([2017 3 22 19 0 0]) );
  plot(pH(f,1),pH(f,2),'k.');
  pH(f,2) = NaN;
  clear f;

%--%

%QC pCO2w
 pCO2 = [dn + (data(:,1)/86400/1000) data(:,3:4)];
 clear data;
 %missing values
  f = find(pCO2(:,2) == 0);
  pCO2(f,:) = [];
  clear f;
 %plot
  figure('Units','inches','Position',[0 0 11 8.5]);
  ha = tight_subplot(3,1,0.02,0.05,0.05);
  axes(ha(1));
    hold on;
    plot(pH(:,1),pH(:,2),'.-');
  axes(ha(2));
    hold on;
    plot(pCO2(:,1),pCO2(:,2),'.-');
 %manual outliers
  f1 = find( pCO2(:,1)>datenum([2016 8 29 18 55 0]) & ...
    pCO2(:,1)<datenum([2016 8 29 19 5 0]) );
  f2 = find( pCO2(:,1)>datenum([2016 11 20 19 0 0]) & ...
    pCO2(:,1)<datenum([2016 11 20 19 5 0]) );
  plot(pCO2([f1 f2],1),pCO2([f1 f2],2),'k.');
  pCO2([f1 f2],2) = NaN;
  clear f1 f2;

 %compare pH and pCO2w
  %index pCO2 to pH
   dt = 15/1440;
   xt = round(pCO2(:,1)/dt)*dt;
   yt = round(pH(:,1)/dt)*dt;
   [ti,ia,ib] = intersect(xt,yt);
   x = pCO2(ia,2);
   y = pH(ib,2);
   clear dt xt yt ia ib;
  %calc theoretical relation
   sali = interp1(CTD(:,1),CTD(:,3),ti);
   tempi = interp1(CTD(:,1),CTD(:,2),ti);
   alk = 2150.0 + 44.0.*(sali-31.25); %(GF)
  %CO2SYS
   sil = 0;
   po4 = 1;
   A = CO2SYS(alk,y,1,3,sali,tempi,NaN,0,NaN,...
     sil,po4,1,4,1);
   xPred = A(:,4); %(uatm)
   clear sali tempi alk sil po4 n A;
  %plot diff
   axes(ha(1));
     plot(ti,y,'g.');
   axes(ha(2));
     plot(ti,x,'g.');
     plot(ti,xPred,'c.-');
   axes(ha(3));
     hold on;
     plot(ti,x-xPred,'.-');
     plot(ti([1 end]),[0 0],'k');
   linkaxes(ha,'x');
  clear ha xPred;

  %calc theoretical relation
   sali = [33 34];
   tempi = [8 15];
   alk = 2150.0 + 44.0.*(sali-31.25); %(GF)
  %CO2SYS
   sil = 0;
   po4 = 1;
   pH_th = [7.7:0.01:8.7];
   for n=1:length(sali)
     A = CO2SYS(alk(n),pH_th,1,3,sali(n),tempi(n),NaN,0,NaN,...
       sil,po4,1,4,1);
     xPred(:,n) = A(:,4); %(uatm)
   end
   clear sali tempi alk sil po4 n A;
  %plot scatter
   figure; hold on;
   plot(xPred(:,1),pH_th,'-','Color',[.4 .4 .4],'LineWidth',2);
   plot(xPred(:,2),pH_th,'k-','LineWidth',2);
   plot(x,y,'b.');

  % Definitely some clusters that look worse than others, appears to be
  %  associated with broad-scale cooling beginning mid-November. Could be
  %  a different water mass (Davidson Current?) with different TA.
  %  May be difficult to correct for. Timing is not related to the sensor
  %  swap 3/22/2017.

   clear x y ti pH_th xPred

 %QC pCO2a (no manual outliers)
  figure;
  plot(pCO2(:,1),pCO2(:,3),'.-');

%--%

%MET
 MET = csvread([mydir 'MET.csv'],1,0);
 MET(:,1) = dn + (MET(:,1)/86400/1000);

%get rid of non-unique time stamps (?)
 [~,i] = unique(MET(:,1));
 MET = MET(i,:);
 clear i;

%QC atm press
 figure; hold on;
 plot(MET(:,1),MET(:,2),'.-');
 %fixed floor/ceiling
  f = find( MET(:,2)<980 | MET(:,2) >1040 );
  plot(MET(f,1),MET(f,2),'k.');
  MET(f,2) = NaN;
  clear f;

%QC air temp
 figure; hold on;
 plot(MET(:,1),MET(:,3),'.-');
 %fixed floor
  f = find( MET(:,3)<7 );
  plot(MET(f,1),MET(f,3),'k.');
  MET(f,3) = NaN;
  clear f;

%--%

%METwind
 METw = csvread([mydir 'METwind.csv'],1,0);
 METw(:,1) = dn + (METw(:,1)/86400/1000);

%get rid of non-unique time stamps (?)
 [~,i] = unique(METw(:,1));
 METw = METw(i,:);
 clear i;

%QC wind speed
 figure; hold on;
 plot(METw(:,1),METw(:,3),'.-');
 %missing values
  f = find(METw(:,3) == 0);
  plot(METw(f,1),METw(f,3),'r.');
  METw(f,2:3) = NaN;
  clear f;
 %fixed ceiling
  f = find( METw(:,3)>30 );
  plot(METw(f,1),METw(f,3),'k.');
  METw(f,2:3) = NaN;
  clear f;

%combine MET and METwind into same fmt as older QC'd data sets
% somehow have more values for METw than MET (?)
 MET = [MET NaN(size(MET,1),1) MET(:,3) ...
   interp1(METw(:,1),METw(:,2),MET(:,1),'nearest') ...
   interp1(METw(:,1),METw(:,3),MET(:,1),'nearest')];

%--%

%VR2C
 VR2C = csvread([mydir 'tagID.csv'],1,0);
 VR2C(:,1) = dn + (VR2C(:,1)/86400/1000);

%leave duplicate time stamps (multiple tags at same time)

%QC VR2C
 %fixed ceiling
  f = find( VR2C(:,2) > 2e5 );
  VR2C(f,2) = NaN;
  clear f;

%--%

%save to interns directory
 save([atlas 'jsevadjian/interns/Diego/mat/OA1_201607_QC.mat'],...
   'CTD', 'Fluor', 'MET', 'O2', 'pCO2', 'pH', 'VR2C');