%  This file computes the regression and correlation maps of SST and SLP
%  onto the monthly, or winter-averaged CTI.  This is used in HW3 problem
%  3 (Spring 2008).  All data are defined as November through March
%  (either monthly or averaged), from 1949 through 2002 (year corresponds
%  to January for each year.  Hence, November 1948 - March 1949
%  corresponds to 1949).
%

%  Start by clearing the workspace.
clear


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Next, load data for the monthly-average regression maps
load CTI_monthly_data.mat

%  For winter-averaged, uncomment the following four lines.
% load CTI_winter_data.mat
% cti_mon = cti_win;
% skt_mon = skt_win;
% slp_mon = slp_win;

lims = [0 360 -90 90];  %  [minlon maxlon minlat maxlat].  This is the
                        %  region over which data is defined (I just
                        %  happen to know this).

%  Take a look at what variables are in the workspace by typing 'who', or
%  'whos': 

whos   %  Long list of variables in the workspace

%  Note that the variable 'description' is in the workspace.  This is a
%  string variable that contains a description of all the other
%  variables.  To see the variable 'description', type 'description' on
%  the command line.  (Note that I created this variable when I generated
%  the .mat file in the first place - it's not built in or anything).

description

%  Let's start by generating regression maps of SST onto the CTI.  Start
%  by standardizing the CTI (place a semi-colon on the end of the line to
%  avoid seeing the results of the operation):

cti = (cti_mon - mean(cti_mon)) / std(cti_mon);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Make regression maps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Now, let's get the regression map of of SST onto the CTI.  First, we
%  need to get the SST data into matrix format.  Currently, it's a
%  3-dimensional array, where the first, second and third dimensions are
%  time, latitude, and longitude, respecitvely (in general, NetCDF
%  climate data have dimension time, level, latitude, longitude):

[ntim, nlat, nlon] = size(skt_mon);

%  Now, reshape into a two-dimensional matrix (time X space):

skt_mon = reshape(skt_mon, ntim, nlat*nlon);
size(skt_mon)  %  Just to double-check

%  I've already removed the monthly climatology (hence the mean) of the
%  data, so we can go ahead and get the regression maps as the covariance
%  of the cti with each column of skt_mon.  Note that the variance of
%  'cti' is one, so we don't need to divide by this in the regression
%  equation:  a1_map = cti' * skt_mon / (cti' * cti);

a1_map = cti' * skt_mon / (ntim - 1);

%  The correlation coefficient involves another step:  First, we need to
%  divide each column of skt_mon by its standard deviation, then get the
%  covariance.  To divide each element in skt_mon, use the element
%  division (which has a '.' in front of the '/':  ./):

corr_map = cti' * (skt_mon ./ (ones(ntim,1) * std(skt_mon))) / (ntim - 1);

%  Finally, reshape a1_map and corr_map so that it has dimensions (lat X
%  lon):

a1_map = reshape(a1_map, nlat, nlon);
corr_map = reshape(corr_map, nlat, nlon);

%  Now, let's calculate the same thing for SLP.
[ntim, nlat, nlon] = size(slp_mon);
slp_mon = reshape(slp_mon, ntim, nlat*nlon);

a1_map_slp = cti' * slp_mon / (ntim-1);
corr_map_slp = cti' * (slp_mon ./ (ones(ntim,1) * std(slp_mon))) / (ntim - 1);

a1_map_slp = reshape(a1_map_slp, nlat, nlon);
corr_map_slp = reshape(corr_map_slp, nlat, nlon);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Make Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Now, we're ready to generate a contour plot of this.  First, open a
%  figure window - I'm going to use the 'tall' orientation, and place two
%  maps in each window (using the subplot command):

figure(1);
orient tall;
wysiwyg;
subplot(2,1,1);
cla;

%  Set up a map axis - we'll do the contouring using the mapping
%  toolbox.

axesm('eqdcylin', 'maplatlimit', [-90 90], 'maplonlimit', [0 360], ...
      'origin', [0 180]);

%  First, define a 'graticule mesh' that has the latitude and longitude
%  points for the SST data we'll use.
[lat2d, lon2d] = meshgrat(lat_skt, lon_skt);  

%  Now, make a map surface of the SST regression map.
hh1s = surfacem(lat2d, lon2d, a1_map, -1*ones(size(a1_map)));
caxis([-1 1]);

%  Contour the SLP regression over the top.  Make positive contours
%  solid, and dashed contours negative.

hold on;  %  Don't erase what's already in the plot
[cc1, hh1] = contourm(lat_slp, lon_slp, a1_map_slp, [0.5:0.5:2.5], '-k');
[cc2, hh2] = contourm(lat_slp, lon_slp, a1_map_slp, [-2.5:0.5:-0.5], '--k');
[cc3, hh3] = contourm(lat_slp, lon_slp, a1_map_slp, [0 0], '-k');
set(hh3, 'linewidth', 2);  %  Make the 0 contour bold
hold off;

%  Make it pretty
lm = worldlo('POpatch');
for i = 1:length(lm);
  lm(i).altitude = -0.5;
  lm(i).otherproperty = {'facecolor', 0.7*[1 1 1], 'edgecolor', 'none'};
end

displaym(lm);  %  Add land map
gridm on;      %  Add grid lines
tightmap;      %  'tighten' image to limits of the subplot
title('\bfDJF monthly:  SST and SLP regression maps', 'fontsize', 14);
xlabel(['SLP contour:  0.5 mb (std dev)^-^1 ;  SST shading:  0.2 ^oC (std' ...
        ' dev) ^-^1'], 'fontsize', 12);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Now, plot the correlation map below this.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,1,2);
cla;

%  Set up map axis and make surface / contour maps.
axesm('eqdcylin', 'maplatlimit', [-90 90], 'maplonlimit', [0 360], ...
      'origin', [0 180]);

[lat2d, lon2d] = meshgrat(lat_skt, lon_skt);  
hh1s = surfacem(lat2d, lon2d, corr_map, -1*ones(size(a1_map)));
caxis([-1 1]);

hold on;  
[cc1, hh1] = contourm(lat_slp, lon_slp, corr_map_slp, [0.2:0.2:1], '-k');
[cc2, hh2] = contourm(lat_slp, lon_slp, corr_map_slp, [-1:0.2:-0.2], '--k');
[cc3, hh3] = contourm(lat_slp, lon_slp, corr_map_slp, [0 0], '-k');
set(hh3, 'linewidth', 2);  %  Make the 0 contour bold
hold off;

%  Make it pretty
lm = worldlo('POpatch');
for i = 1:length(lm);
  lm(i).altitude = -0.5;
  lm(i).otherproperty = {'facecolor', 0.7*[1 1 1], 'edgecolor', 'none'};
end

displaym(lm);  %  Add land map
gridm on;      %  Add grid lines
tightmap;      %  'tighten' image to limits of the subplot
title('\bfDJF monthly:  SST and SLP regression maps', 'fontsize', 14);
xlabel(['SLP contour:  0.5 mb (std dev)^-^1 ;  SST shading:  0.2 ^oC (std' ...
        ' dev) ^-^1'], 'fontsize', 12);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Repeat the above set of commands for the winter-averaged data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%