% keeling_fit.m -- script to try and manage real pCO2 data from the
% Mauna Loa Observatory 
%
% Download the data from:
% (http://scrippsco2.ucsd.edu/data/atmospheric_co2)


% set up the inputs

data_file='monthly_mlo.csv'	% the monthly file is easier to parse
skip_header_lines=57; 		% choose a line with data after all the headers

% read in the input data

indata=csvread(data_file,skip_header_lines,0);
[ndata,ncolumns]=size(indata); 
ndata=ndata-4; 			% remove the padding at the end of the file
				% (in the data I downloaded, the last four lines
				% had no data, so I am excluding them from the 
				% count -- your data file may be different...)

dates_decimal=indata(1:ndata,4);
ts_start=dates_decimal(1);
pCO2=indata(1:ndata,5);

% and plot them!
% (the dates in decimal years are in column 4; the raw pCO2 values are in 
%  column 5)

figure;
plot(dates_decimal,pCO2); hold on
xlabel('date')
ylabel('pCO2 (ppm)')
title('Mauna Loa Observatory uncorrected CO2 data')

% without doing any additional data conditioning, it seems there are holes in the
% data set! indeed, if you look at the data file, there are multiple lines with 
% '-99.99' as their pC02 values -- these are data gaps that we might want to fill

pCO2_interp=zeros(ndata,1);

% loop through the data and look for holes!

for i=3:ndata             % start at row 3 since there are two holes before it!
                          
   if(pCO2(i)==-99.99)    % if we find a hole...

      disp(sprintf('hole at line %d',i));

      test_data=pCO2(i-3:i+3);		% grab some data either side
      test_dates=dates_decimal(i-3:i+3);	% (choose your window wisely...)
      interp_date = dates_decimal(i);

      test_data_cleaned=nonzeros((test_data~=-99.99).*test_data);   % remove the hole(s) from
      test_dates_cleaned=nonzeros((test_data~=-99.99).*test_dates); % these subsets

      interp_data=pchip(test_dates_cleaned,test_data_cleaned,interp_date); % and interpolate!
%     interp_data=spline(test_dates_cleaned,test_data_cleaned,interp_date); % and interpolate!
      pCO2_interp(i)=interp_data;       			           % (I think pchip works better here)

      disp(sprintf('estimated pCO2 value %f  at epoch %f',interp_data,interp_date)); disp(' ');

%      if you want to check the interpolation, uncomment these
%      figure;
%      plot(test_dates_cleaned,test_data_cleaned,'b.'); hold on
       plot(interp_date,interp_data,'r.')  % can leave this uncommented to plot interpolated values
%      title(sprintf('hole at %f',interp_date))

   else                  % if no hole

      pCO2_interp(i)=pCO2(i);

   end

end

% remove the first two rows, since there's no data there, recount the number of data points

dates_decimal(1:2,:)=[];
pCO2_interp(1:2,:)=[];
ndata=length(dates_decimal);

% plot your interpolated time series

hold off;
figure;
plot(dates_decimal,pCO2_interp,'r')
xlabel('date')
ylabel('pCO2 (ppm)')
title('Mauna Loa Observatory cleaned CO2 data')
hold on  % any subsequent plot commands will plot on top of this...

% and now they should look a lot better!

% ------

% next: let's fit some straight lines to the data

% extract the data and dates 
earlyyearsfilter=(dates_decimal>=1959).*(dates_decimal<=1965);
earlyyearsdates=dates_decimal(find(earlyyearsfilter))-ts_start;
earlyyearsdata=pCO2_interp(find(earlyyearsfilter));
% assemble your matrix and do the backslash thing
A=[earlyyearsdates ones(size(earlyyearsdates))];
mearly=A\earlyyearsdata;

% extract the data and dates 
lateyearsfilter=(dates_decimal>=2000).*(dates_decimal<=2015);
lateyearsdates=dates_decimal(find(lateyearsfilter))-ts_start;
lateyearsdata=pCO2_interp(find(lateyearsfilter));
% assemble your matrix and do the backslash thing
A=[lateyearsdates ones(size(lateyearsdates))];
mlate=A\lateyearsdata;

% and display the results

disp(sprintf('1959-1965 average rate: %f',mearly(1)));
disp(sprintf('2000-2015 average rate: %f',mlate(1)));

% plot yer straight lines...
predicted_dates=(1958:2040)';
predicted_dates_shifted=predicted_dates-ts_start;
earlyyearsline=[predicted_dates_shifted ones(size(predicted_dates_shifted))]*mearly;
plot(predicted_dates,earlyyearsline,'k');
lateyearsline=[predicted_dates_shifted ones(size(predicted_dates_shifted))]*mlate;
plot(predicted_dates,lateyearsline,'b');


% ------
% next, why not fit some polynomials, too!

% example: estimate a function and calculate its misfit
%polyparams2=polyfit(dates_decimal-ts_start,pCO2_interp,2);	% fit a parabola
%pCO2_bestfit2=polyval(polyparams2,dates_decimal-ts_start);	% evaluate said parabola
%misfit2=sqrt((pCO2_interp-pCO2_bestfit2)'*(pCO2_interp-pCO2_bestfit2));
%disp(sprintf('parabola misfit: %f ppm',misfit2));

% example: fit different polynomials, extrapolate and plot them
% predicted_dates=1958:1:2040;
% polyparams2=polyfit(dates_decimal-ts_start,pCO2_interp,2);	% fit a parabola
% pCO2_bestfit2=polyval(polyparams2,predicted_dates-ts_start);	% evaluate said parabola
% plot(predicted_dates,pCO2_bestfit2,'k')
% polyparams3=polyfit(dates_decimal-ts_start,pCO2_interp,3);	% fit a cubic
% pCO2_bestfit3=polyval(polyparams3,predicted_dates-ts_start);	% evaluate said cubic
% plot(predicted_dates,pCO2_bestfit3,'y')
% polyparams4=polyfit(dates_decimal-ts_start,pCO2_interp,4);	% fit a quartic
% pCO2_bestfit4=polyval(polyparams4,predicted_dates-ts_start);	% evaluate said quartic
% plot(predicted_dates,pCO2_bestfit4,'c')
% polyparams5=polyfit(dates_decimal-ts_start,pCO2_interp,5);	% fit a 5th order poly
% pCO2_bestfit5=polyval(polyparams5,predicted_dates-ts_start);	% evaluate it
% plot(predicted_dates,pCO2_bestfit5,'m')
% polyparams6=polyfit(dates_decimal-ts_start,pCO2_interp,6);	% fit a 6th order poly
% pCO2_bestfit6=polyval(polyparams6,predicted_dates-ts_start);	% evaluate it
% plot(predicted_dates,pCO2_bestfit6,'g')
