EOF分解

EOF分解，用ncdisp命令可以看到HadISST_sst.nc是1870年1月至最近的全球SST数据，按比例截取北太平洋（20N-80N 130E-90W）1971年1月至1990年12月240个月份的SST数据进行经验正交函数（Empirical Othorgnal Function）分解,简记为EOF分解，得到该区域该时段的海温时空特征。

clear
%ncdisp HadISST_sst.nc
 %time = 1730  （UNLIMITED）
 %latitude  = 180
 %longitude = 360
 %nv        = 2
%units = ‘days since 1870-1-1 0:0:0‘

lon=ncread（‘HadISST_sst.nc‘‘longitude‘）;%提取经度数据
lat=ncread（‘HadISST_sst.nc‘‘latitude‘）; %提取纬度数据
sst=ncread（‘HadISST_sst.nc‘‘sst‘）;      %提取海表面sst数据% 读取nc格式数据

sst1=sst（1:9021:801201:1440）;          %选取所需要区域的数据
sst2=sst（311:36021:801201:1440）;       %130E-90W20N-80N
sst3=zeros（14060240）;
sst3（90:-1:11:601:240）=sst1;
sst3（140:-1:911:601:240）=sst2;
sst=sst3;                                          

sst_area1=zeros（2408400）;                %zeros全零数组
for i=1:240;
 squ=squeeze（sst（::i））;                 %执行该指令后sst数据转换为二维数组                         
 sst_area1（i:）=reshape（squ18400）;      %将数据转变为二维
end 

% 剔除与其它站点相关系数小的站点的数据
sst_area1（sst_area1==-1000）=NaN;% 简单的认为剔除陆地和冬季结冰点的数据
sst_nan=isnan（sst_area1）;
i=0;
for j=1:8400
    if sum（sst_nan（:j））==0;
        i=i+1;
        sst_region（:i）=sst_area1（:j）;
    end
end
%**************************************************************************
% 求距平~注意季节的变换
X=zeros（size（sst_region））;          % 学者给的程序
for i=1:12
X（i:12:228+i:）=sst_region（i:12:end:） - repmat（ mean（sst_region（i:12:end:）1）  size（sst_region（i:12:end:）1） 1）;
end

% 由于变量数（8400）>观测样本数~协方差矩阵的阶数较大，因此可用时空转换的方法提高计算速度
                                                       
R=X*X‘;         % 协方差矩阵R=X*X‘是8400*8400的方阵~现定义矩阵R=X‘*X是240*240的矩阵
[vd]=eig（R）;   % 进行EOF分解~因为X‘*X与X*X‘的秩相同所以特征值相同~d为x的特征值组成的对角阵~v为X*X‘的特征向量
v=fliplr（v）;    % 矩阵作左右翻转
d=rot90（d2）;   % 矩阵上下翻转后再左右翻转（查看生成的对角阵是由小到大排列的，用命令doc rot90可查询
diagonal=diag（d）;
spacef=X‘*v;                                       
for i=1:240;
    spacef（:i）=spacef（:i）/sqrt（diagonal（i））;     % 空间本征函数
end
timef=X*spacef;                                    % 时间本征函数
sum_d=sum（diagonal）;
count=0;
for i=1:240;
    count=count+diagonal（i）;           
    G1（i）=count/sum_d;                            % G1（i）是累积方差贡献率
end
for i=1:240;
    G2（i）=diagonal（i）/sum_d;                     % G2（i）是方差贡献率
end
                                                 % 将删去的陆地与冰点的填充值补回
sst_area2=zeros（2408400）;
sst_area2（::）=NaN;
i=1;
spacef2=spacef‘;
for j=1:8400
   if sum（sst_nan（:j））==0;
      sst_area2（:j）=spacef2（:i）;
      i=i+1;
   end
end
sst_area3=sst_area2‘;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                
%画图

x=1:240;
plot（xtimef（:6）‘b‘）;
ylim（[ -80 80 ]）;
% xlabel（‘1971-1990年240个月‘‘fontsize‘15‘fontname‘‘隶书‘） 
ylabel（‘INDEX‘‘fontsize‘12‘fontname‘‘黑体‘）
set（gca‘xtick‘（1:6:162））
set（gca‘xticklabel‘{‘1971‘‘‘‘1972‘‘‘‘1973‘‘‘‘1974‘‘‘‘1975‘‘‘‘1976‘‘‘‘1977‘‘‘‘1978‘‘‘‘1979‘‘‘‘1980‘‘‘‘1981‘‘‘‘1982‘‘‘‘1983‘‘‘‘1984‘‘‘‘1985‘‘‘‘1986‘‘‘‘1987‘‘‘‘1988‘‘‘‘1989‘‘‘‘1990‘‘‘‘1991‘}）
title（‘北太平洋第6模态1971至1990年SST时间序列‘ ‘color‘ ‘k‘‘fontsize‘15‘fontna

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件      337557  2014-04-24 23:32  010022011028蓝方宝.docx
     文件       14363  2014-04-24 20:19  010022011028lfb\图片\001.png
     文件       14334  2014-04-24 20:38  010022011028lfb\图片\002.png
     文件       14316  2014-04-24 20:38  010022011028lfb\图片\003.png
     文件       14305  2014-04-24 20:39  010022011028lfb\图片\004.png
     文件       14305  2014-04-24 20:40  010022011028lfb\图片\005.png
     文件       14310  2014-04-24 20:42  010022011028lfb\图片\006.png
     文件       61709  2014-04-24 19:34  010022011028lfb\图片\011.png
     文件       61474  2014-04-24 19:35  010022011028lfb\图片\022.png
     文件       65711  2014-04-24 19:40  010022011028lfb\图片\033.png
     文件       59179  2014-04-24 19:42  010022011028lfb\图片\044.png
     文件       69114  2014-04-24 19:44  010022011028lfb\图片\055.png
     文件       68440  2014-04-24 19:45  010022011028lfb\图片\066.png
     文件        4372  2014-04-24 20:02  010022011028lfb\程序\newEOF.m
     文件         918  2014-04-24 19:47  010022011028lfb\程序\new_draw.m
     文件         445  2014-04-24 20:41  010022011028lfb\程序\timef_draw.m