eof

xueqixiang

家园的各位老师,用HADISST海温数据作了EOF处理,结果如下,第一模态时间序列出错了.请老师们指点迷津.谢谢
sst = ncread('HadISST_sst.nc','sst');                         % 读取nc格式数据
%**************************************************************************
ran = 372;
sst1=sst(1:110,40:101,1405:1776);                        % 选取所需要区域的数据
sst2=sst(241:360,40:101,1405:1776);
sst3=zeros(230,62,ran);
sst3(110:-1:1,1:62,1:ran)=sst1;
sst3(230:-1:111,1:62,1:ran)=sst2;
sst=sst3;                                          
%**************************************************************************
sst_area1=zeros(ran,14260);                             % zeros全零数组
for i=1:ran;
squ=squeeze(sst(:,:,i));                              % 执行该指令后sst数据转换为二维数组                        
sst_area1(i,:)=reshape(squ,1,14260);                   % 将数据转变为二维
end
%**************************************************************************
% 剔除与其它站点相关系数小的站点的数据~简单的认为剔除陆地和冬季结冰点的数据
sst_area1(sst_area1<-10000)=NaN;                   % 陆地和冰点的填充值为-1.00000001504747e+30~将此值定义为NaN


    sst_nan=isnan(sst_area1);

i=0;
for j=1:14260
    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:ran-12+i,:)=sst_region(i:12:end,:) - repmat( mean(sst_region(i:12:end,:),1) , size(sst_region(i:12:end,:),1), 1);
end


%**************************************************************************


R=X*X';                                          % 协方差矩阵R=X*X'是14260*14260的方阵~现定义矩阵R=X'*X是ran*ran的矩阵
[v,d]=eig(R);                                          % 进行EOF分解~因为X'*X与X*X'的秩相同所以特征值相同~d为x的特征值组成的对角阵~v为X*X'的特征向量~
v=fliplr(v);                                           % 矩阵作左右翻转
d=rot90(d,2);                                          % 矩阵上下翻转后再左右翻转(查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr(flipud(d))=rot90(d,2))
diagonal=diag(d);
spacef=X'*v;                                       
for i=1:ran;
    spacef(:,i)=spacef(:,i)/sqrt(diagonal(i));         % 空间本征函数
end
timef=X*spacef;                                        % 时间本征函数
sum_d=sum(diagonal);
count=0;
for i=1:ran;
    count=count+diagonal(i);           
    G1(i)=count/sum_d;                                 % G1(i)是累积方差贡献率
end
for i=1:ran;
    G2(i)=diagonal(i)/sum_d;                           % G2(i)是方差贡献率
end
%**************************************************************************
                                                       % 将删去的陆地与冰点的填充值补回
sst_area2=zeros(ran,14260);
sst_area2(:,:)=NaN;
i=1;
spacef2=spacef';
for j=1:14260
   if sum(sst_nan(:,j))==0;
      sst_area2(:,j)=spacef2(:,i);
      i=i+1;
   end
end
sst_area3=sst_area2';
%**************************************************************************                                                     
% 画图


subplot(2,2,2)          % 将绘图窗口划分为2*1个子窗口，并在第2个子窗口中绘图
x=1:ran;
plot(x,timef(:,1),'b');

set(gca,'xtick',(1:60:ran));
set(gca,'xticklabel',{'1987','1992','1997','2002','2007','2012','2017'});
title('北太平洋第1模态1987至2017年SST时间序列', 'color', 'k','fontsize',12,'fontname','幼圆');
grid on;
hold off

subplot(2,2,1)


lon=[60.5:289.5];
lat=[-10.5:50.5];
m_proj('Equidistant cylindrical','long',[60.5 289.5],'lat',[-10.5 50.5]);
m_contourf(lon,lat,rot90(reshape(sst_area3(:,1),230,62)',2));
m_coast('patch',[1 .85 .7]);
m_elev('contourf',[500:500:6000]);
m_grid('box','fancy','tickdir','in');
colormap(flipud(copper));

title(['北太平洋第1模态填色图'],'fontsize',12,'fontname','幼圆');

subplot(2,2,3)
lon=[60.5:289.5];
lat=[-10.5:50.5];
m_proj('Equidistant cylindrical','long',[60.5 289.5],'lat',[-10.5 50.5]);
m_contourf(lon,lat,rot90(reshape(sst_area3(:,2),230,62)',2));
% colorbar;
m_coast('patch',[1 .85 .7]);
m_elev('contourf',[500:500:6000]);
m_grid('box','fancy','tickdir','in');
colormap(flipud(copper));

title(['北太平洋第2模态填色图'],'fontsize',12,'fontname','幼圆');

subplot(2,2,4)          % 将绘图窗口划分为2*1个子窗口，并在第2个子窗口中绘图
x=1:ran;
plot(x,timef(:,2),'b');

set(gca,'xtick',(1:60:ran));
set(gca,'xticklabel',{'1987','1992','1997','2002','2007','2012','2017'});
title('北太平洋第2模态1987至2017年SST时间序列', 'color', 'k','fontsize',12,'fontname','幼圆');
grid on

xueqixiang

又跑了一遍,各位老师能指点下吗?是数据还是程序问题,搞不出来了
为什么第二模态时间序列是正常的

膘膘

范围画小一点

xueqixiang

膘膘 发表于 2019-10-23 19:26
范围画小一点

第二模态是正常的，

膘膘

xueqixiang 发表于 2019-10-24 16:11
第二模态是正常的，

SST这种东西，你设nan值太多，矩阵不一定可逆。这个第二模态正不正常没有关系。还有就是你数据距平或者标准化没有

xueqixiang

膘膘 发表于 2019-10-25 13:49
SST这种东西，你设nan值太多，矩阵不一定可逆。这个第二模态正不正常没有关系。还有就是你数据距平或者标 ...

谢谢，改成中低纬度了，可能数据本身在高纬度异常值太多。可以了，谢谢

xueqixiang

膘膘 发表于 2019-10-25 13:49
SST这种东西，你设nan值太多，矩阵不一定可逆。这个第二模态正不正常没有关系。还有就是你数据距平或者标 ...

谢谢，改成中低纬度了，可能数据本身在高纬度异常值太多。可以了，谢谢

yijianmingliu

学习一下，我也要相同问题

希诺漾

楼主，新手刚刚学EOF，不能读懂每个含义，想先试着运行一下，能把数据发我一下吗?或者加下我的QQ，指教一下我，急，有偿感谢   

希诺漾

楼主，新手刚刚学EOF，不能读懂每个含义，想先试着运行一下，能把数据发我一下吗?或者加下我的QQ，指教一下我，急，有偿感谢3208945419@qq.com  会EOF的大佬欢迎来!急，有偿