Global Index (all files) (short | long) | Local Index (files in subdir) (short | long)
Demo for calculating the first EOF of NH monthly SLP (the annular mode), and displaying regression maps of SLP and surface temperature onto the leading principal component time series. This assumes that the NCEP reanalysis monthly mean data is somewhere that is accessible.
| This script calls | |
|---|---|
% Start by clearing the workspace. I always do this.
clear
% Start with loading SLP data from a NetCDF file
% Define limits for loading SLP data
lims = [0 360 20 90]; % [minlon maxlon minlat maxlat];
tim = get_time(1950, 1999, 1948); % Month indices for 1950-1999. Input
% is yr1, yr2, start year for the
% entire data set.
lev = 1;
filename = '/Volumes/Data/NCEP/monthly/slp.mon.mean.nc';
varname = 'slp';
% Load data
[slp, lat, lon] = getnc2(filename, varname, lims, lev, tim);
% Now, remove the annual cycle from SLP:
[slp, slpclim] = remove_anncyc(slp);
% Before performing EOF/PC analysis, we need to weight SLP by
% sqrt(cos(lat)), so that the covariance matrix is weighted by
% cos(lat).
slp = cosweight(slp, lat);
% Now, get EOFs.
[lam, lds, pcs, per] = eof_dan(slp);
% Let's plot the leading two EOF's side by side on a sheet of paper, in
% a figure that will fit into a Journal of Climate format (6.25"
% across).
% Set up the figure window, and define size of the axes (plots):
figure_tall(1); % Sets up figure with orientation tall
wysiwyg; % What you see is what you get
% Now, define subplot size. Let's make it 3" by 3" with a 0.25" space
% in between. We'll do side by side. Also, a 1.5" margin from the
% top, and allow 0.75" between these two plots and the next two below
% them.
%
% The following routine sets some global variables that are used by the
% subplot2 command. This just defines the size of the subplots we'd
% like.
global_axes(3, 3, 0.25, 0.75, 1.5); % hsz, vsz, midh, midv, tmarg. Type
% 'help global_axes' for more info.
% Now, plot away!
clf; % clear current figure
% Open subplot 1:
subplot2(2,1); % There are two columns, this is the first.
cla; % Clear this subplot (in case it's full).
% Now, the loadings (EOF's) from above are in column format, so let's
% reshape them into maps, and get them ready for plotting.
nlat = length(lat); nlon = length(lon);
pat = reshape(lds(:,1), nlat, nlon);
pat = uncosweight(pat, lat);
pat = pat * sqrt(lam(1));
% Ok, now plot this. Let's start with a contour plot. Start by
% defining the global variables for the plot (lat, lon, frame limits).
global_latlon(lat, lon, [0 360 20 90]);
% Define map axis:
map_axis('stereo', [90 0]); % Origin is the pole, and Grenwich
% And, contour, using a contour interval of 0.25 mb
map_contour_pn(pat, 0.25, 'nozero');
% Add land, grid lines, etc:
fill_landmap('under');
gridm on;
tightmap;
framem on;
set(gca, 'visible', 'off');
t1 = title('{\bf NAM:} SLP reg. onto PC1', 'fontsize', 14, 'visible', 'on');
x1 = xlabel('Contour: 0.25 mb (std dev)^-^1', 'fontsize', 12, ...
'visible', 'on');
% There you have it! Now, let's do the same for the second EOF, the
% Aleutian low:
% Open subplot 2:
subplot2(2,2); % There are two columns, this is the second.
cla; % Clear this subplot (in case it's full).
nlat = length(lat); nlon = length(lon);
pat = reshape(lds(:,2), nlat, nlon);
pat = uncosweight(pat, lat);
pat = pat * sqrt(lam(2));
% We do not need to redefine the global lat and lon variables.
map_axis('stereo', [90 0]); % Origin is the pole, and Grenwich
map_contour_pn(pat, 0.25, 'nozero');
fill_landmap('under');
gridm on;
tightmap;
framem on;
set(gca, 'visible', 'off');
t1 = title('{\bf Aleutian Low:} SLP reg. onto PC2', ...
'fontsize', 14, 'visible', 'on');
x1 = xlabel('Contour: 0.25 mb (std dev)^-^1', 'fontsize', 12, ...
'visible', 'on');
% OK, now that we have all this, let's plot the correlation map of
% surface temperature onto the NAM and onto the Aleutian low.
filename2 = '/Volumes/Data/NCEP/monthly/skt.mon.mean.nc';
varname2 = 'skt';
[skt, lat2, lon2] = getnc2(filename2, varname2, lims, 1, tim);
[skt, sktclim] = remove_anncyc(skt);
% Now, get regressions:
[pat1, c1] = regress_map(skt, pcs(:,1));
[pat2, c2] = regress_map(skt, pcs(:,2));
global_latlon(lat2, lon2, [0 360 20 90]);
subplot2(2,3); cla;
map_axis('stereo', [90 0]);
map_surface_interp(pat1, -1);
caxis([-1 1]*1.5);
% Now, mask out the ocean:
fill_oceanmap(-0.5, [1 1 1]);
% And, draw in land:
draw_landmap(0.3*[1 1 1]);
gridm on;
tightmap;
framem on;
set(gca, 'visible', 'off');
t1 = title('{\bf NAM:} SKT reg. onto PC1', ...
'fontsize', 14, 'visible', 'on');
% And the same for the Aleutian Low
subplot2(2,4); cla;
map_axis('stereo', [90 0]);
map_surface_interp(pat2, -1);
caxis([-1 1]*1.5);
% Now, mask out the ocean:
fill_oceanmap(-0.5, [1 1 1]);
% And, draw in land:
draw_landmap(0.3*[1 1 1]);
gridm on;
tightmap;
framem on;
set(gca, 'visible', 'off');
t1 = title('{\bf Aleutian Low:} SKT reg. onto PC1', ...
'fontsize', 14, 'visible', 'on');
% Finally, add a colorbar to the bottom of this figure. We can control
% the size through the global_axes command:
global_axes(5, .15, .25, 0.75, 1.5+3+0.75+3+0.15);
subplot2(1,1); % Note, we just made the top margin big enough that it
% places the subplot below all the existing graphs.
cb = colorbar(gca);
set(gca, 'XTick', 1:11, 'XTickLabel', -1.5:.3:1.5);
x1 = xlabel('^oC (std dev)^-^1', 'fontsize', 12);
% Print the results somewhere.
pushd /Users/dvimont/matlab/Examples
print -depsc2 -painters NAM_Aleutian_regressions.eps
popd
% Last, let's perform spectral analysis on this time series. I feel
% comfortable with tem degrees of freedom, which will require that we
% chop the data set into 5 segments (for fifty years, this means the
% lowest allowable frequency will be 10 yrs.
timeseries = pcs(:,1);
ntim = length(timeseries);
nwind = ntim/5; % This has to be an integer
dof = 2*ntim/nwind; % Degrees of freedom
nover = nwind/2; % This also has to be an integer
[p, f] = spectrum(timeseries, nwind, nover, hanning(nwind), 12);
% Now, define best fit AR1 process for this
rho1 = corr(timeseries, timeseries, 1);
pAR = (1-rho1^2)./(1+rho1^2-2*rho1*cos(2*pi*[0:(nwind/2)]/nwind));
% Normalize pAR so it represents equal variance as p.
pAR = pAR * sum(p(2:end,1))/sum(pAR(2:end));
% Get confidence interval:
f95 = finv(0.95, dof, 100000);
f05 = finv(0.05, 100000, dof);
% OK! let's plot this four different ways.
figure_tall(2); clf;
global_axes(3, 1.5, .5, .75, 1.5);
% Start with power vs. frequency - this way variance is conserved.
subplot2(1,1); cla;
pp = plot((f), p(:,1), 'k', ...
(f), pAR*f95, '--k', ...
(f), pAR*f05, '--k');
set(pp(1), 'linewidth', 2);
set(gca, 'XTick', 0:6);
axis([0 6 0 2500]);
yl = ylabel('Power', 'fontsize', 10);
xl = xlabel('Frequency (cycles yr^-^1)', 'fontsize', 10);
t1 = title2('{\bfPC1 Spectrum:} Power vs. Frequency', 'fontsize', 12);
grid on;
% Now, power vs. log(f) - this stretches out the low frequencies, but
% variance is not proportional to area under the curve. Note that
% log(f) is proportional to period.
subplot2(1,2); cla;
pp = plot(log(f), p(:,1), 'k', ...
log(f), pAR*f95, '--k', ...
log(f), pAR*f05, '--k');
set(pp(1), 'linewidth', 2);
set(gca, 'XTick', log(2.^[-5:2]), 'XTickLabel', 2.^[5:-1:-2]);
axis([log(1/10)-.1 log(6)+.1 0 2500]);
yl = ylabel('Power', 'fontsize', 10);
xl = xlabel('Period (yr)', 'fontsize', 10);
t1 = title2('{\bfPC1 Spectrum:} Power vs. Log(f)', 'fontsize', 12);
grid on;
% Now, log(power) vs. log(f) - this stretches out the low frequencies,
% and compresses the variability (visually) in the estimates, but
% variance is not proportional to area under the curve. Note that
% log(f) is proportional to period.
subplot2(1,3); cla;
pp = plot(log(f), log(p(:,1)), 'k', ...
log(f), log(pAR*f95), '--k', ...
log(f), log(pAR*f05), '--k');
set(pp(1), 'linewidth', 2);
set(gca, 'XTick', log(2.^[-5:2]), 'XTickLabel', 2.^[5:-1:-2]);
axis([log(1/10)-.1 log(6)+.1 4.5 8.2]);
yl = ylabel('log(Power)', 'fontsize', 10);
xl = xlabel('Period (yr)', 'fontsize', 10);
t1 = title2('{\bfPC1 Spectrum:} log(Power) vs. Log(f)', 'fontsize', 12);
grid on;
% Now, f.*power vs. log(f) - this stretches out the low frequencies, and
% variance is proportional to area under the curve. It's a little
% unorthodox - some people have a tough time interpreting it. But, it
% also produces a visual maximum at 2*pi times the damping rate (for an
% AR1 process) and at the principal oscillation frequency for an AR2
% process. This is better for strongly autocorrelated time series.
subplot2(1,4); cla;
pp = plot(log(f), f.*p(:,1), 'k', ...
log(f), f'.*pAR*f95, '--k', ...
log(f), f'.*pAR*f05, '--k');
set(pp(1), 'linewidth', 2);
set(gca, 'XTick', log(2.^[-5:2]), 'XTickLabel', 2.^[5:-1:-2]);
axis([log(1/10)-.1 log(6)+.1 0 4000]);
yl = ylabel('f * Power', 'fontsize', 10);
xl = xlabel('Period (yr)', 'fontsize', 10);
t1 = title2('{\bfPC1 Spectrum:} f * Power vs. Log(f)', 'fontsize', 12);
grid on;