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

%updated 12.2019 by Diego Sancho to compare prior to deployment using test data
%from parking lot.

% Air temp from the 200WX is biased due to short power-on time.
%  Here, we compare air temp from the 200WX and the nearby MLML
%  weather station, and derive a relationship to correct 200WX air temp.

%define buoy, deployment, data directory
 buoy = 'oa1';
 dep = '202010';

 atlas = '//atlas.shore.mbari.org/';
 mydirOA = [atlas 'oasis_coredata/deployments/' ...
    buoy '/' dep '/' buoy];
%%
%read-in mlml .csv files
 
 mlml_9 = csvread([atlas 'OA_Moorings/cal_data/MET/MLML_2020-09.csv'],1,2); %ignores first column w UTC time
 mlml_10 = csvread([atlas 'OA_Moorings/cal_data/MET/MLML_2020-10.csv'],1,2);
%  mlml_5 = csvread([atlas 'OA_Moorings/cal_data/MET/MLML_2020-05.csv'],1,2);
%  mlml_6 = csvread([atlas 'OA_Moorings/cal_data/MET/MLML_2020-06.csv'],1,2);
 mlml=[mlml_9;mlml_10];
 
 %read in text files and concatenate 
 
 %how many days back you want to go from last day of test data
%%
 numdays=10; %take number of days
 
 D=dir([atlas 'oasis_coredata/deployments/' buoy '/' dep '/' buoy '/pb200']);
 
 %loads each line of txt file into cell
 tlines=cell(0,1);
 for i=length(D)-numdays:length(D)
     filename=D(i).name;
     file=fopen([mydirOA '/pb200/' filename]);
     tline=fgetl(file);
     while ischar(tline)
         tlines{end+1,1}=tline;
         tline=fgetl(file);
     end
     fclose(file);
 end
  %%
 %clean up non-data lines
 for i=1:length(tlines)
     f(i)=strcmp('#',tlines{i,1}(1));
 end
 tlines(f,:)=[];
 
 %extract time and air T from char string
times=str2double(cellfun(@(tlines) tlines(1:10), cellstr(tlines),'un',0));
mettime=datenum(2020,1,1,0,0,0)+times; %wx200 samples every hour

%split by commas
splitline=split(tlines,',');
metairt=str2double(splitline(:,20));
        
 
%convert MLML time in Unix to Matlab datenum
 dn = datenum([1970 1 1]);
% met(:,1) = dn + met(:,1)/86400/1000;
 mlml(:,1) = dn + mlml(:,1)/86400;
 clear dn;
%%
%simple time-series plots
 figure; hold on;
 plot(mlml(:,1),mlml(:,6),'b.-');
 plot(mettime,metairt,'g.-');
 
 %%
 %clean mlml data
 f=find(mlml(:,6)<5 | mlml(:,6)>40);
 mlml(f,6)=nan; 
 
 
%%
%look at time-averaged (smoothed data)

%align times--round to nearest
d1=datetime(mettime,'ConvertFrom','datenum'); 
met_t=dateshift(d1,'start','minute');%removes trailing seconds

d=datetime(mlml(:,1),'ConvertFrom','datenum');
mlml_t=dateshift(d,'start','minute'); %round to nearest minute (down) :05

%find mlml indices in each hour
mlmltavg=nan(length(met_t)-1,1);
for i=1:length(met_t)-1
    f=find(mlml_t>=met_t(i) & mlml_t<met_t(i+1)); 
    mlmltavg(i)=mean(mlml(f,6),'omitnan');
end
%%

%NaN empty MET data
 %met(met(:,3)==0,3) = NaN;
 
%align
 for n=1:size(mettime,1)
   [m(n),i(n)] = min(abs(mettime(n) - mlml(:,1)));
 end
 f = find(m < 30/1440);
 x = metairt(f,1);
 y = mlml(i(f),6);
 time=mettime(f);
 
 clear n m i f;

%get rid of NaN values
 f = find(sum(isnan([x y]')) > 0);
 x(f) = [];
 y(f) = [];
 time(f)= [];
 clear f;

%scatter plot
 figure; hold on;
 plot(x,y,'k.');

%get regression
 pf = polyfit(x,y,1); % y = 0.4123x + 8.2120

%check
 pv = polyval(pf,[min(x) max(x)]);
 plot([min(x) max(x)],pv,'b-');
 
 %%
 %convert to new T, compare to old T and MLML t
 newT=x.*pf(1)+pf(2);
 figure
 plot(mettime,metairt,'k')
 hold on
 plot(time,newT,'r')
 hold on
 plot(time,y,'g')
 
 %%
 %compare calibrations (to fewer points--oldT)
 old_pf=[0.3393 9.2791];
 oldT=x.*old_pf(1)+old_pf(2);
 
 figure
 plot(time,newT,'g')
 hold on
 plot(time,oldT,'k')