%  Tutorial for reconstructing decadal ENSO-like pattern, using only
%  interannual EOF's (as in Vimont, 2005).
%

%  Start by clearing the workspace.
clear


%  Step 1:  Load data

%  Start by defining data set and parameters (you could do this on the
%  command line, I like to keep track of these if I save data into a .mat
%  file).
filin = 'skt.mon.mean.nc';
varn = 'skt';
lims = [110 290 -30 60];
lev = 1;
tim = get_time(1948, 2002, 1948);

%  Now, load the data
cd ~/matlab/Data
[sst, lat, lon] = getnc2(filin, varn, lims, lev, tim);

%  Get the size of the sst data set - we'll use this later
[ntim, nlat, nlon] = size(sst);

%  Load land mask
lm = getnc2('land.sfc.gauss.nc', 'land', lims, 1, 1);
land = find(lm == 1);

%  Now, remove points that are ever covered with ice, and add to the land
%  mask
icept = find(min(sst) < 0);
land = sort(union(land, icept));

%  From here on, we can speed up calculations by setting land values to
%  NaN's, and just dealing with SST data
sst(:,land) = NaN;

%  Redefine ocean points here
ocept = find(~isnan(sst(1,:,:)));
sst2 = sst(:,ocept);


%  Step 2:  remove means, and apply consecutive running mean filters.
%  This latter step is much quicker if we define the filter weights, and
%  convolve the filter weights with the data.

%  Get rid of annual cycle
sst2 = remove_anncyc(sst2);

%  Apply 6-yr lowpass filter to SST
%  First, construct lowpass filter by convolving two running mean filters
%  (25 month and 27 month):
fweights = conv(ones(25,1)/25, ones(37,1)/37);

%  nf is the length of the filter.  This helps define the impulse
%  response function.  As this is a non-recursive filter, the first and
%  last nf/2-1 data points are affected by edge effects.
nf = length(fweights);

%  Next, initialize new data set, sst_lp (this makes calculations
%  faster - you should usually do this when defining a new large matrix
%  row-by-row).
[ntim, npts] = size(sst2);
sst_lp = repmat(NaN, [ntim, npts]);

%  Finally, apply the filter.  Note that the convolution has length
%  ntim+nf-1 (it adds values on to each end of the time series).  The
%  first (nf-1)/2 points in the filtered time series are affected by the
%  ends, so we have to trim 'nf-1' points from the beginning and end of
%  the convolved time series.

kpt = (nf-1)/2 + [1:(ntim-nf+1)];  %  indices of points that are
                                   %  unaffected by the edges 
for i = 1:npts;
  tem = conv(fweights, sst2(:,i));
  sst_lp(kpt,i) = tem(kpt + (nf-1)/2);
end

%  Define highpass-filtered SST by removing lowpass-filtered SST from
%  original SST data
sst_hp = sst2 - sst_lp;


%  Step 3:  Apply EOF analysis to the data

%  First, weight data by the cosine of latitude.  Do this by swapping
%  sst_lp or sst_hp with a 3D array called 'tem', so we can use the
%  'cosweight' routine.  (Recall, we only have ocean points, so 'tem' has
%  all points).
tem = repmat(NaN, [ntim nlat nlon]);

tem(:,ocept) = sst_lp;
tem = cosweight(tem, lat);
sst_lp = tem(:,ocept);

tem(:,ocept) = sst_hp;
tem = cosweight(tem, lat);
sst_hp = tem(:,ocept);

%  OK!  Now, let's get the EOF's.  'eof_dan' is my own routine, which
%  I'll post on the web.
[laml, ldsl, pcsl, perl] = eof_dan(sst_lp(kpt,:));
[lamh, ldsh, pcsh, perh] = eof_dan(sst_hp(kpt,:));


%  Step 4:  Obtain reconstructed EOFs

%  Now, project highpass-filtered patterns onto lowpass-filtered data.
%  This produces a set of low-frequency time series that describe how the
%  eof's of high-frequency data vary on low-frequency time scales.  Note
%  that these are _not_ equal to the pcsl above, so they're not
%  orthogonal or anything.
pcsh_lp = sst_lp * ldsh;

%  Now, we'd like to try to reconstruct the low-frequency time series /
%  spatial pattern using only the patterns from the high-frequency
%  patterns.  So, perform EOF analysis on pcsh_lp.  BUT, don't use all of
%  them, or they will span the space of the low-frequency data, and we'll
%  exactly reproduce ldsl for purely spurious reasons.
kp_mode = 1:4;

%  Now, run the EOF routine on the reconstructed time series.  Note that
%  'eof_dan' assumes you only want the first 10 modes, so indicate that
%  you only want 4 (there aren't 10 time series).
[lamr, ldsr, pcsr, perr] = eof_dan(pcsh_lp(kpt,kp_mode), 4);

%  Finally, get the reconstructed spatial patterns.
ldsr_pat = ldsh(:,kp_mode)*ldsr;


%  Step 5:  Plot results

%  Start by defining a figure window, and the size of the subplots we'll
%  be using (these are defined to fit in AMS column format:  3.125 in
%  wide.  The height is defined by trial and error):
figure_tall(1); clf;
global_axes(3.1, 1.8, .5, .2, 1.5);

%  Define global variables 'XAX', 'YAX', and 'FRAME' which will be used
%  in future routines:
global_latlon(lat, lon, lims);

%  Start with plotting first EOF of highpass-filtered and
%  lowpass-filtered SST.  Note that we need to add land points back to
%  the patterns, so initialize 'pat1_hp' as a matrix of NaN's:
pat1_hp = repmat(NaN, [nlat nlon]);

%  Replace ocean points with the first EOF of hp-filt SST (scaled).  Note
%  polarity of ldsh is switched for plotting purposes.
pat1_hp(ocept) = -1*ldsh(:,1)*sqrt(lamh(1));

%  Finally, reverse the sqrt(cos) weighting
pat1_hp = uncosweight(pat1_hp, lat);

%  Now, plot it up!  Define subplot
subplot2(1,1); cla

%  Make subplot a 'map' axis:
map_axis('eqdcylin');

%  Contour the map
map_contour_pn(pat1_hp, 0.1, 'zero');

%  And, make it pretty:
fill_landmap('over');
gridm on;
tightmap;

title2('{\bfEOF1:}  Highpass Filtered SST (32%)', 'fontsize', 12);

%  Do the exact same for the Lowpass-filtered EOF1
pat1_lp = repmat(NaN, [nlat nlon]);
pat1_lp(ocept) = -1*ldsl(:,1)*sqrt(laml(1));
pat1_lp = uncosweight(pat1_lp, lat);

subplot2(1,2); cla;
map_axis('eqdcylin');
map_contour_pn(pat1_lp, 0.1, 'zero');
fill_landmap('over');
gridm on;
tightmap;

title2('{\bfEOF1:}  Lowpass Filtered SST (41%)', 'fontsize', 12);

%  Finally, plot the reconstructed EOF
pat1_r = repmat(NaN, [nlat nlon]);
pat1_r(ocept) = 1*ldsr_pat(:,1)*sqrt(lamr(1));
pat1_r = uncosweight(pat1_r, lat);

subplot2(1,3); cla;
map_axis('eqdcylin');
map_contour_pn(pat1_r, 0.1, 'zero');
fill_landmap('over');
gridm on;
tightmap;

title2('{\bfEOF1:} Reconstructed LP\_SST (69%)', 'fontsize', 12);

%  Add units
xlabel('Contour:  0.1 ^oC', 'fontsize', 10);

%  Finally, plot time series

%  Redefine subplot axes so it fits at the very bottom.  Note, we just
%  change the top margin so it's below all the previous subplots, then
%  define our subplot as the first one.
global_axes(3.1, 1.5, .5, .2, 1.5+3*1.8+2*0.2+0.5);
subplot2(1,1); cla;

%  Now, plot PC's, the pcsh above, the pcsl and pcsr below.  Note that
%  we retain the scaling on the PC's, even though eventually we won't
%  know the units.  Visually, the variations are comparable.
yr = 1948:1/12:2003.999;
pp = plot(yr(kpt), -pcsh(:,1)+50, 'k', ...
	  yr(kpt), pcsr(:,1)+0, '--k', ...
	  yr(kpt), -pcsl(:,1)+0, '-k');

%  Define axis limits, and put a grid on
axis([1949 2001 -30 110]);
set(gca, 'YTick', [0 50], 'YTickLabel', [], 'XTick', 1950:10:2000);
grid on;

%  Change the line width:
set(pp, 'linewidth', 1);

%  Set the color of the reconstructed pc to grey, and the linewidth to 2,
%  so it's distinguishable:
set(pp(2), 'color', 0.6*[1 1 1], 'linewidth', 2);

%  Finally, add a title and xlabel
t2 = title2(['{\bfPC1:}  HP (top), LP (bot), LPR (bot)'], 'fontsize', 12);
xlabel('Year', 'fontsize', 10);

%  For kicks, let's add labels to the curves (this is an interactive
%  plotting function):
t1a = gtext('\bfHPC1', 'fontsize', 12);
t2a = gtext('\bfRPC1', 'fontsize', 12, 'color', 0.6*[1 1 1]);
t3a = gtext('\bfLPC1', 'fontsize', 12);


%  Step 6:  Print it up!  I prefer to print as an encapsulated
%  post-script, using the 'painters' algorithm.  It makes really nice
%  printer-quality plots with little memory, and you can convert these to
%  jpeg using 'gs' (ask me for a quick routine to do this if you'd
%  like).  Also, you can change the length of the dashes, which enhances
%  final print quality.

cd ~/matlab/Figs
print -depsc2 -painters tutorial_fig1.eps