% calculates average of monthly chl readings
% compares float and satellite data to check for biases due to the larger
  % volume of float data
  
% loads:
  % float0005Surface.mat (generated by applyCalSurface.m)  
  % modis_chl, containing satellite chl data (downloaded 01/17/08 by Peter Brickley)
  
% saves cloudChl, A_chl, std_chl, med_chl, n_chl, C_chl

% written by E. Boss, UMaine, 2007
% email: emmanuel.boss@maine.edu

close all
clear all

load float0005Surface
% float0005Surface contains date, lat, lon, P, Chl_b_corr, B_LSS_b_corr, stat

% get first reading of the day
I=find(diff(float0005Surface(:,1))>4);
time_surf=[float0005Surface(1,1);float0005Surface(I+1,1)];
chl_surf=[float0005Surface(1,5);float0005Surface(I+1,5)];

% get average float chl concentration per month
[yy,month,dd] = datevec(time_surf); %converts date string to date vector
mmyy = [yy month chl_surf]; %finds each year (once) in dataset and stores in array uyear
uyear = unique(yy); %finds each year (once) in dataset and stores in array uyear
CHL_APEX = []; %initializes array CHL_APEX
for ky = 1:length(uyear) %for each year in dataset
    umonth = unique(month(yy==uyear(ky))); %finds each month (once) in each year
    for km = 1:length(umonth) %for each month in each year in each dataset
        midx = yy==uyear(ky) & month==umonth(km); %where year AND month exist (create logical array)
        Chl_b_midx = chl_surf(midx); %create array chl_b_midx of chl data at midx = true
        N=length(find(midx==1)); %get number of float points in the month
        thisrow = [uyear(ky) umonth(km) nan(1,2)]; %make array thisrow of yr, mo, NaNs for each chl_surf
        thisrow(3) = nanmedian(Chl_b_midx); %add median chl concentration for month to thisrow column 3
        thisrow(4) = nanmean(Chl_b_midx);  %add mean chl concentration for month
        thisrow(5) = N; %add number of points to column 5
        thisrow (6) = std(Chl_b_midx);
        CHL_APEX = [CHL_APEX; thisrow]; %append new row to big array CHL_APEX
    end
end

% get average modis chl concentration per month
load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\0005\MODIS\modis_chl'
% modis_chl contains:
  % modis_chl_ts (2116x2601) - chl data only from 06/13/04 - 01/01/07
  % modisChl [date, chl] (2116x2602)
  % distgrid (2116x2601) grid of distance around float
modis_date=modisChl(:,1);

% get mean and median modis values within 7.5km of float
[uu vv]=size(modis_chl_ts); % 2116x2601
I=find(distgrid(1,:)<7.5);
for i=1:uu
    dist_A(i)=nanmean(modis_chl_ts(i,I),2);
    dist_B(i)=nanmedian(modis_chl_ts(i,I),2);
end

% get all float points at time of satellite
Ikk=find(finite(dist_B));
float_sat_chl=interp1(time_surf,chl_surf,modis_date(Ikk),'linear'); % float at time of modis
I=finite(float_sat_chl);
interp_chl=float_sat_chl(I);
B=modis_date(Ikk);
time_interp=B(I);

% get average chl concentration per month    
[yy,month,dd] = datevec(time_interp); %converts date string to date vector
mmyy = [yy month interp_chl]; %finds each year (once) in dataset and stores in array uyear
uyear = unique(yy); %finds each year (once) in dataset and stores in array uyear
CHL_APEX_C = []; %initializes array CHL_APEX_C
for ky = 1:length(uyear) %for each year in dataset
    umonth = unique(month(yy==uyear(ky))); %finds each month (once) in each year
    for km = 1:length(umonth) %for each month in each year in each dataset
        midx = yy==uyear(ky) & month==umonth(km); %where year AND month exist (create logical array)
        Chl_b_midx = interp_chl(midx); %create array chl_b_midx of chl data at midx = true
        N=length(find(midx==1)); % get number of satellite points in the month
        thisrow = [uyear(ky) umonth(km) nan(1,2)]; %make array thisrow of yr, mo, NaNs for each chl_surf
        thisrow(3) = nanmedian(Chl_b_midx); 
        thisrow(4) = nanmean(Chl_b_midx); 
        thisrow(5) = N; % add number of points to column 5
        CHL_APEX_C = [CHL_APEX_C; thisrow]; %append new row to big array CHL_APEX (see plot.osd.m for CHL_APEX contents)
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTS

A=datenum(CHL_APEX(:,1),CHL_APEX(:,2),15);  
C=datenum(CHL_APEX_C(:,1),CHL_APEX_C(:,2),15);  

% plot(A(2:27),CHL_APEX(2:27,4),'*b',C,CHL_APEX_C(:,4),'*k')
% plot(A(2:27),(CHL_APEX(2:27,4)-CHL_APEX_C(:,4))./CHL_APEX(2:27,4),'*b')

% figure(1);
% subplot(2,1,1)
% % plotyy(A,CHL_APEX(:,6)./CHL_APEX(:,3),C,CHL_APEX_C(:,5)) % standard deviation of float and satellite coverage (normalized) vs time
% [AX,H1,H2] = plotyy(A,CHL_APEX(:,6)./CHL_APEX(:,3),C,CHL_APEX_C(:,5),'plot');
% set(gca,'FontName','Times New Roman');
% set(gca,'FontSize',12);
% set(get(AX(1),'Ylabel'),'String',{'Standard Deviation/Median' 'of Float Chlorophyll'});
% % set(get(AX(2),'Ylabel'),'String','Number of Matchups per Month');
% set(AX(1),'xlim',[min(time_surf-1) max(time_surf)]);
% set(AX(2),'xlim',[min(time_surf-1) max(time_surf)]);
% set(AX(1),'XTick',[732115 732480 732850 733215]);
% set(AX(2),'XTick',[]);
% % set(0,'DefaultAxesColorOrder',[0 0 0]);
% set(H1,'LineStyle','-');
% set(H2,'LineStyle','-');
% set(H1,'Marker','.');
% set(H2,'Marker','*');
% set(H1,'LineWidth',1);
% set(H2,'LineWidth',1);
% datetick(AX(1),'x', 'mmmyy','keeplimits','keepticks');
% set(gca,'XTick',[]);
% % title({'Standard Deviation in Float',' and MODIS Coverage vs Time'});

% figure(1);
% subplot(2,1,2);
% plot(A,CHL_APEX(:,5),'k.-',C,CHL_APEX_C(:,5),'*k-','LineWidth',1); % plots number of float and satellite points read per month
% set(gca,'FontName','Times New Roman');
% set(gca,'FontSize',12);
% set(gca,'XTick',[732115 732480 732850 733215]);
% datetick(gca,'x','mmmyy','keepticks');
% xlim([min(time_surf-1) max(time_surf)]);
% % xlabel('date');
% ylabel({'Number of Matchup' 'per Month'});
% legend('Float','Modis');

figure('Units','inches','Resize','off','position',[0 0 3.5 3],'paperposition',[0.5 4 3.5 3]);
axes1 = axes('Position',[0.187 0.167 0.65 0.74]);
% [AX,H1,H2] = plotyy(tim,totalchl',tim,totalbb','Parent',axes1);
% plotyy(A,CHL_APEX(:,6)./CHL_APEX(:,3),C,CHL_APEX_C(:,5)) % standard deviation of float and satellite coverage (normalized) vs time
H3=plot(A,CHL_APEX(:,5),'ko');
hold on
[AX,H1,H2] = plotyy(C,CHL_APEX_C(:,5),A,CHL_APEX(:,6)./CHL_APEX(:,4),'plot','Parent',axes1);
plot([732205 732254],[23.8 23.8],'k','LineWidth',1,'Marker','none'); % legend symbol
h=plot(732230,22.00,'Marker','o','MarkerSize',[8],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor','k'); % legend symbol
h=plot(732230,20.14,'Marker','s','MarkerSize',[8],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor','k'); % legend symbol
hold off
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
set(get(AX(1),'Ylabel'),'String',{'Number of samples per month'},'Color','k');
set(get(AX(2),'Ylabel'),'String',{'Coefficient of variation of float chlorophyll'},'Rotation',270.0,'VerticalAlignment','bottom','Color','k');
xlim([732090 733250]);
tick = [732099 732464 732829 733194 733560];
set(AX(1),'XTick',tick,'XTickLabels', datestr(tick,'dd mmm yy'),'YColor',[0 0 0]);
set(AX(2),'xlim',[min(time_surf-1) max(time_surf)],'XTick',[],'YColor',[0 0 0]);
set(H1,'Marker','s','MarkerFaceColor','none','MarkerEdgeColor','k','LineStyle','none');
set(H2,'LineWidth',2,'LineStyle','-','Color','k','Marker','none');
xlabel('Date');

% draw legend
annotation('textbox',[0.28 0.815 0.1 0.1],'String','CV of float chlorophyll','EdgeColor','none','FontName','Times New Roman','FitHeightToText','on');
annotation('textbox',[0.28 0.765 0.1 0.1],'String','float samples','EdgeColor','none','FontName','Times New Roman','FitHeightToText','on');
annotation('textbox',[0.28 0.710 0.1 0.1],'String','satellite samples','EdgeColor','none','FontName','Times New Roman','FitHeightToText','on');

corrcoef(CHL_APEX_C(1:35,5),CHL_APEX(2:36,6)./CHL_APEX(2:36,4));

% subplot(2,1,2);
% plot(CHL_APEX_C(:,5),100*(CHL_APEX(2:36,4)-CHL_APEX_C(:,4))./CHL_APEX(2:36,4),'*b'); % mean
% hold on
% plot(CHL_APEX_C(:,5),100*(CHL_APEX(2:36,3)-CHL_APEX_C(:,4))./CHL_APEX(2:36,3),'*g'); % median
% hold off
% xlabel('Number of Modis Satellite Data Points in a Given Month');
% ylabel('% Deviation in Float Chlorophyll Due to Sampling');
% legend('Mean','Median');

% plot all float data, not just averages:
% plotyy(time_surf,chl_surf,C,CHL_APEX_C(:,5)) % float chl and satellite coverage vs time

std_chl=CHL_APEX(:,6);
med_chl=CHL_APEX(:,3);
n_chl=CHL_APEX_C(:,5);
A_chl=datenum(CHL_APEX(:,1),CHL_APEX(:,2),15); 
C_chl=datenum(CHL_APEX_C(:,1),CHL_APEX_C(:,2),15); 

save cloudChl A_chl std_chl med_chl n_chl C_chl
