% Compare TMI OI daily and NOAA Reynolds weekly SST.
%
% Simon de Szoeke
%

% for TMI data:
lonc =   0.125:0.25:359.875; % center of the patches
latc = -89.875:0.25:89.875;
lonrange = [0 360];
latrange = [-40 40];
lon = lonc(lonc>=lonrange(1) & lonc<=lonrange(2));
lat = latc(latc>=latrange(1) & latc<=latrange(2));
scale_factor = 0.15;
add_offset   = -3;
%averaging matrices for .25->1.0 degree
j = 1:length(lon)';
i = ceil(j/4);
s = ones(1,1440)/4;
Alon = sparse(i,j,s,360,1440);
j = 1:length(lat)';
i = ceil(j/4);
s = ones(1,320)/4;
Alat = sparse(i,j,s,80,320);

% for Reynolds data
% see http://www.cdc.noaa.gov/cdc/data.noaa.oisst.v2.html
nc = netcdf('/home/sol/deszoeke/data/oi_sst/Reynolds/sst.wkmean.1990-present.nc');
time = nc{'time'}+datenum('1-1-1800'); % Sundays (beginning of avg interval)
timec = time+3; % Wednesdays (center of interval)
latr = nc{'lat'}(:); % descending! 89.5:-1:-89.5 (180)
lonr = nc{'lon'}(:); % .5:1:359.5 (360)

start = datenum('1-may-2001')-datenum('0-jan-2001');
endt  = datenum('31-oct-2001')-datenum('0-jan-2001'); % 1may->31oct!
%start = datenum('1-aug-2001')-datenum('0-jan-2001');
%endt  = datenum('30-sep-2001')-datenum('0-jan-2001'); % aug-sep
dt = 7;
sm    = zeros(80,360);
Rey1  = sm;
TMI1  = sm;
vari1 = sm;

switch 'load' % read and *process* the data again (slow), or *load* it from a file
 case 'process'
  % Start: read the previous and next weekly SST for Wednesdays surrounding
  % the first day, start.
  mday = start+datenum('0-jan-2001');
  mdayp = mday - mod(mday-timec(1),dt);
  indp = (mdayp-timec(1))/dt + 1; % f(1:n) index
  % read the current week sst (indp)
  sstp = nc{'sst'}(indp,:,:)*nc{'sst'}.scale_factor(:) + nc{'sst'}.add_offset(:);
  sstp(sstp<nc{'sst'}.valid_range(1) | sstp>nc{'sst'}.valid_range(2)) = NaN;
  % read the next week sst (indp+1)
  sstn = nc{'sst'}(indp+1,:,:)*nc{'sst'}.scale_factor(:) + nc{'sst'}.add_offset(:);
  sstn(sstn<nc{'sst'}.valid_range(1) | sstn>nc{'sst'}.valid_range(2)) = NaN;
  % loop over desired days for jday
  for jday = start:endt
    mday  = jday+datenum('0-jan-2001');
    mdayp = mday - mod(mday-timec(1),dt);
    if mday==mdayp % mday is a Wednesday, increment weekly data
      indp = (mdayp-timec(1))/dt + 1;
      sstp = sstn; % first timestep sstn undefined!
		   % read the next week sst (indp+1)
		   sstn = nc{'sst'}(indp+1,:,:)*nc{'sst'}.scale_factor(:) + nc{'sst'}.add_offset(:);
		   sstn(sstn<nc{'sst'}.valid_range(1) | sstn>nc{'sst'}.valid_range(2)) = NaN;
    end
    % read the daily TMI sst -> t
    fname = ['/home/sol/deszoeke/data/oi_sst/tmi/2001/tmi.fusion.2001.',sprintf('%03d',jday),'.v01'];
    fid = fopen(fname);
    temp = reshape(fread(fid,1440*720,'char'),1440,720);
    fid = fclose(fid);
    temp(temp>250) = NaN;
    t = temp(:,latc>=latrange(1) & latc<=latrange(2))*scale_factor+add_offset;
    t1x1 = Alat*(Alon*t)';
    
    flat = flipud(find(abs(latr)<40));
    ti = sstp(flat,:) + (sstn(flat,:)-sstp(flat,:))*mod(mday-timec(1),dt)/dt;
    mask = isfinite(t1x1) & isfinite(ti);
    TMI1(mask)  = TMI1(mask)  + t1x1(mask);
    Rey1(mask)  = Rey1(mask)  + ti(mask);
    vari1(mask) = vari1(mask) + (t1x1(mask)-ti(mask)).^2;
    sm = sm + mask;
  end
  TMI = TMI1./sm;
  Rey = Rey1./sm;
  bias = TMI-Rey;
  vari = vari1./sm - bias.^2;
  
  mask = sm>(2*dt);
  TMI(~mask)  = NaN;
  Rey(~mask)  = NaN;
  vari(~mask) = NaN;
  bias(~mask) = NaN;

%  save Reynolds_TMI_vari_augsep.mat Rey TMI vari sm
  save Reynolds_TMI_vari.mat Rey TMI vari sm
  
case 'load'  % load the data
  flat = flipud(find(abs(latr)<40));
  load Reynolds_TMI_vari
  bias = TMI-Rey;
end

% plots
% mean TMI daily and Reynolds weekly
subplot(2,1,1)
contourf(lonr,latr(flat),TMI,10:32);
caxis([10 32]);colorbar;
title({'TMI daily (C)'})
ylabel('degrees latitude')
subplot(2,1,2)
contourf(lonr,latr(flat),Rey,10:32);
caxis([10 32]);colorbar;
title('Reynolds weekly (C)')
xlabel('degrees E longitude')
orient landscape
print -depsc TMI_vs_Reynolds_mean.eps
% zoom in on E. Pacific
%flonEP = find(lonr>=230 & lonr<=280);
%flatEP = find(latr>=-15 & latr<=15);
%flEP = flipud(find(latr(flat)>=-15 & latr(flat)<=15));
flonEP = find(lonr>=360-140 & lonr<=360-40);
flatEP = find(latr>=-30 & latr<=30);
flEP = flipud(find(latr(flat)>=-30 & latr(flat)<=30));

subplot(2,1,1)
contourf(lonr(flonEP),latr(flatEP),TMI(flEP,flonEP),17:29);
caxis([17 29]);colorbar;
title({'TMI daily (C)'})
ylabel('degrees latitude')
subplot(2,1,2)
contourf(lonr(flonEP),latr(flatEP),Rey(flEP,flonEP),17:29);
caxis([17 29]);colorbar;
title('Reynolds weekly (C)')
xlabel('degrees E longitude')
orient landscape
print -depsc TMI_vs_Reynolds_EPacific_mean.eps

% look at the differences
subplot(2,1,1)
contourf(lonr,latr(flat),bias,-2:.25:2);
caxis([-2 2]); colorbar;
title({'SST_{TMI daily} - SST_{Reynolds weekly} bias (K)'})
ylabel('degrees latitude')
subplot(2,1,2)
contourf(lonr,latr(flat),sqrt(vari),0:.25:2);
caxis([0 1.5]); colorbar;
title('RMS difference (K)')
xlabel('degrees E longitude')
orient landscape
print -depsc TMI_vs_Reynolds_diff.eps
% zoom in on E. Pacific
subplot(2,1,1)
contourf(lonr(flonEP),latr(flatEP),bias(flEP,flonEP),-2:.25:2);
caxis([-1 1]); colorbar;
title({'SST_{TMI daily} - SST_{Reynolds weekly} mean bias (K)'})
ylabel('degrees latitude')
subplot(2,1,2)
contourf(lonr(flonEP),latr(flatEP),sqrt(vari(flEP,flonEP)),0:.2:2);
caxis([0 1.2]); colorbar;
title('RMS difference (K)')
xlabel('degrees E longitude')
orient tall
print -depsc TMI_vs_Reynolds_EPacific_diff.eps