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

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

% QC's pH only, for now. Later should update to do pCO2 at the same time,
%  possibly other M1 data. Based on 'QC_wgTiny.m'

% Outputs a final .csv file in the /M1 pH data directory

%depl info
 buoy = 'm1-slave1';
 dep = '201507';

%directories for post-cruise and QC'd data
 atlas = '//atlas.shore.mbari.org/';
 mydirQC = [atlas 'jsevadjian/M-Buoy/final pH data/'];

%paths to M1 code
 addpath([atlas 'jsevadjian/code/M1/QC code/']);
 addpath([atlas 'jsevadjian/M-Buoy/code/']);

%read M1 ctd data
 ctd = load_M1_CTD('m1',dep);

%read & convert M1 pH data
 pHcal = [ ...
   datenum([2015 7 30 19 0 0; 2015 10 23 18 50 0])', -0.378917, 284.99; ...
   datenum([2015 10 23 19 0 0; 2016 8 29 0 0 0])', -0.384977, 279.15];

 pH = load_M1_pH(buoy,dep,ctd,pHcal);

%load already-QC'd CO2 data
 load([atlas 'jsevadjian/M-Buoy/final CO2 data/' dep '.mat']);

%QC pH:

% (1) Time-series plot
% This is used to manually NaN data, for example, if the pH battery voltage
%  was low, there were issues with the controller, or a dummy pH sensor was used.
 fig_pH = figure('Name','pH','Units','inches','Position',[0 0 11 8.5]);
 ha = tight_subplot(2,1,0.075,0.15,0.1);
 g = get(ha,'Position');
 set(ha,'Box','on','NextPlot','add');
 axes(ha(1));
   plot(ctd.time,ctd.sal,'.-');
   ylabel('');
   title('Salinity');
 axes(ha(2));
   plot(pH.time,pH.pH,'.-');
   ylabel('[total scale]');
   title('pH');
   uicontrol('Units','normalized','Style','pushbutton','String','NaN',...
     'Position',[g{2}(1)+g{2}(3)-0.08 g{2}(2)+g{2}(4)+0.003 .08 .03],...
     'Callback','pH.pH = qcM1_NaNcb(pH.time,pH.pH,ha(2));');
   set(ha,'XLim',[min([ctd.time; pH.time]) max([ctd.time; pH.time])] + [-1 1]*0.1);
 linkaxes(ha,'x');

% (2) Sliding variance filter on pH
% This is used to identify intermittent spikes in pH similar to most
%  salinity time series.
 win = 30; %window
 nst = 3; %number of stdev's above median
 sm = slidingmed(pH.pH, win);
 ss = slidingstdev(pH.pH, win);
 fs = find(abs(pH.pH - sm) > nst*ss);
 clear win nst sm ss;
%check
 axes(ha(2));
 plot(pH.time(fs),pH.pH(fs),'m.');
%NaN bad data
 pH.pH(fs) = NaN;
 clear fs;

% (3) T vs pH
% This is used to identify bad pH data by looking for outliers in
%  water-mass properties. This only 'flags' the outliers; you still need to
%  zoom-in and manually click the 'NaN' button in the time-series plot to
%  NaN the data.
 f = [];
 x = interp1(ctd.time,ctd.temp,pH.time);
 y = pH.pH;
 figTpH = figure('Name','TpH','Units','inches','Position',[0 0 11 8.5]);
 plot(x,y,'.');
 hold on;
 xlabel(['T [' char(186) 'C]']);
 ylabel(['pH']);
 g = get(gca,'Position');
 uicontrol('Units','normalized','Style','pushbutton','String','FLAG',...
   'Position',[g(1)+g(3)-0.08 g(2)+g(4)+0.003 .08 .03],...
   'Callback',...
   ['gx=get(gca,''XLim'');' ...
    'gy=get(gca,''YLim'');' ...
    'f=[f; find(x>=gx(1) & x<=gx(2) & y>=gy(1) & y<=gy(2))];' ...
    'figure(figTpH); plot(x(f),y(f),''k.'');' ...
    'figure(fig_pH); axes(ha(2)); plot(pH.time(f),pH.pH(f),''k.'');']);
 uiwait(figTpH);
 clear f x y g gx gy;

% (4) pCO2w vs pH
% This is used to identify bad pCO2w and/or pH data by comparing to the
%  theoretical relation obtained from CO2SYS using an estimated alkalinity
%  derived from S. This only 'flags' the outliers; you still need to
%  zoom-in and manually click the 'NaN' button in the time-series plot to
%  NaN the data.
%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.5:0.01:8.5];
 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
 f = [];
 x = interp1(co2.time,co2.pCO2water,pH.time);
 y = pH.pH;
 figPP = figure('Name','PP','Units','inches','Position',[0 0 11 8.5]);
 plot(x,y,'.');
 hold on;
 plot(xPred(:,1),pH_th,'-','Color',[.4 .4 .4],'LineWidth',2);
 plot(xPred(:,2),pH_th,'k-','LineWidth',2);
 xlabel('pCO2 [ppm]');
 ylabel('pH');
 g = get(gca,'Position');
 uicontrol('Units','normalized','Style','pushbutton','String','FLAG',...
   'Position',[g(1)+g(3)-0.08 g(2)+g(4)+0.003 .08 .03],...
   'Callback',...
   ['gx=get(gca,''XLim'');' ...
    'gy=get(gca,''YLim'');' ...
    'f=[f; find(x>=gx(1) & x<=gx(2) & y>=gy(1) & y<=gy(2))];' ...
    'figure(figPP); plot(x(f),y(f),''k.'');' ...
    'figure(fig_pH); axes(ha(2)); plot(pH.time(f),pH.pH(f),''k.'');']);
 uiwait(figPP);
 clear x y g gx gy xPred pH_th;

%save QC'd .csv file
 dn = datenum([1970 1 1]);
 time = (pH.time - dn) * 86400; %*1000;
 fid = fopen([mydirQC 'M1_pH_' dep '_QC.csv'],'w');
 fprintf(fid, ['Time[GMT/Unix],pH[TotalScale]']);
 for n=1:length(time)
   fprintf(fid, '\n%.0f,%.3f', [time(n) pH.pH(n)]);
 end
 fclose(fid);
 clear time fid n;

%------------%

%Final Check:
 figFinal = figure('Name','Final','Units','inches','Position',[0 0 11 8.5]);
 ha = tight_subplot(1,1,0,0.3,0.1);
 set(ha,'Box','on','NextPlot','add');
 plot(pH.time,pH.pH,'.-');
 ylabel('[total scale]');
 title('pH');
 set(ha,'XLim',[min(pH.time) max(pH.time)] + [-1 1]*0.1);
 
