% loads float0015Surface, modis_chl, and seawifs_chl, calculates average
   % of monthly readings, and then compares float and satellite data to check
   % for biases due to the larger volume of float data.  Note: have to run this file in parts.

close all
clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part I SEAWIFS

load float0015Surface
% float0015Surface contains:
  % float0015SurfaceChl: [date, chl_corr]

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

% get average 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

% do the same with seawifs data
load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\APEX\0015\Seawifs\seawifs_chl'
% seawifs_chl contains was generated from gen_seawifs_data.m:
  % seawifsChl [date, chl] (1033x170)
  % seawifs_chl_ts (10336x169) just chlorophyll THIS IS A 13X13 GRID
  % distgrid_seawifs (1033x169) grid of distance around float
  % clat (1x1033)
  % clon (1x1033)seawifs_date=seawifsChl(:,1);

seawifs_date=seawifsChl(:,1);
seawifs_chl=seawifsChl(:,2);

[uu vv]=size(seawifs_chl_ts); % 1276x169

% get all points within 7.5km of float
I=find(distgrid_seawifs(1,:)*4<7.5);

for i=1:uu
    dist_A(i)=nanmean(double(seawifs_chl_ts(i,I)),2);
    dist_B(i)=nanmedian(double(seawifs_chl_ts(i,I)),2);
end

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

% get average chl concentration per month    
[yy,month,dd] = datevec(double(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_B = []; %initializes array CHL_APEX_B
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_B = [CHL_APEX_B; 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);
B=datenum(CHL_APEX_B(:,1),CHL_APEX_B(:,2),15);

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

figure(1);
[AX,H1,H2] = plotyy(A,CHL_APEX(:,6)./CHL_APEX(:,3),B,CHL_APEX_B(:,5),'plot');
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',16);
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',[datenum('01 jun 2004') datenum('01 jun 2008')]);
set(AX(2),'xlim',[datenum('01 jun 2004') datenum('01 jun 2008')]);
set(AX(1),'XTick',[732099 732282 732464 732647 732829 733012 733194])
set(AX(2),'XTick',[]);
datetick(AX(1),'x', 'mmmyy','keeplimits','keepticks');
title({'Standard Deviation in Float',' and SeaWIFS Coverage vs Time'});

figure(2);
% subplot(2,1,1);
plot(A,CHL_APEX(:,5),'*b',B,CHL_APEX_B(:,5),'*r-') % plots number of float and satellite points read per month
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',16);
set(gca,'XTick',[732464 732829 733194 733560]);
datetick(gca,'x','mmmyy','keepticks');
ylabel('Number of Data Points');
legend('Float','SeaWIFS');
% subplot(2,1,2);
% plot(CHL_APEX_B(:,5),100*(CHL_APEX(2:27,4)-CHL_APEX_B(:,4))./CHL_APEX(2:27,4),'*b'); % mean
% hold on
% plot(CHL_APEX_B(:,5),100*(CHL_APEX(2:27,3)-CHL_APEX_B(:,4))./CHL_APEX(2:27,3),'*r'); % median
% hold off
% xlabel('Number of Seawifs Satellite Data Points in a Given Month');
% ylabel('% Deviation in Float Chlorophyll Due to Sampling');
% legend('Mean','Median');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part II MODIS

close all
clear all

load float0015Surface
% float0015Surface contains:
  % float0015SurfaceChl: [date, chl_corr]

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

% get average 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

% do the same with modis data
load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\APEX\0015\MODIS\modis_chl'
% modis_chl contains:
  % modis_chl_ts (1105x2601) - chl data only from 08/20/05 - 04/09/07 THIS IS A 51X51 GRID
  % modisChl0015 [date, chl] (1105x2602)
  % distgrid (1105x2601) grid of distance around float

modis_date=modisChl0015(:,1);

%%% get mean and median modis values within specified distances of float by way of Emmanuel's method
[uu vv]=size(modis_chl_ts); % 2116x2601

% get all points within 7.5km of float
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(3); % need to run this plot separately from Figure 4 or you lose axes labels
% 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',16);
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',[datenum('01 jun 2004') datenum('01 jun 2008')]);
set(AX(2),'xlim',[datenum('01 jun 2004') datenum('01 jun 2008')]);
set(AX(1),'XTick',[732099 732282 732464 732647 732829 733012 733194]);
set(AX(2),'XTick',[]);
datetick(AX(1),'x', 'mmmyy','keeplimits','keepticks');
title({'Standard Deviation in Float',' and MODIS Coverage vs Time'});

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

figure(4);
% subplot(2,1,1);
plot(A,CHL_APEX(:,5),'*b',C,CHL_APEX_C(:,5),'*g-'); % plots number of float and satellite points read per month
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',16);
set(gca,'XTick',[732464 732829 733194 733560]);
datetick(gca,'x','mmmyy','keepticks');
% xlabel('date');
% datetick;
ylabel('Number of Data Points');
legend('Float','Modis');
% 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');
