fid=fopen('','r');
Num=360;
data=zeros(360,180,Num);
for i=1:Num
aaa=fscanf(fid,'%s',7);
data(:,:,i)=fscanf(fid,'%f',[360,180]); end
sst1=data(1:90,11:70,1:Num); sst2=data(311:360,11:70,1:Num);
sst3=zeros(140,60,Num); sst3(90:-1:1,1:60,1:Num)=sst1; sst3(140:-1:91,1:60,1:Num)=sst2; sst=sst3;
for i=1:140
for j=1:60
for k=1:Num
if(sst(i,j,k)==-1000)||((sst(i,j,k)==-32768))
sst(i,j,k)=NaN;
end
end
end
end
end X=zeros(size(sst_region));
for i=1:12 X(i:12:Num-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 是 156*156 的矩阵 sst_area1=zeros(Num,8400); for i=1:Num; squ=squeeze(sst(:,:,i)); 换为二维数
组 sst_area1(i,:)=reshape(squ,1,8400);
end 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
% zeros 全零数组 执行该指令后 sst 数据转 将数据转变为二维
选取所需要区域的数据
求距平 ~注意季节的变换
学者给的程序
协方差矩阵 R=X*X' 是 8400*8400
[v,d]=eig(R); % 进行 EOF 分解-因为 X'*X 与X*X'的秩相同所以特征值相同 ~d 为x 的特征值组成的对角阵 ~v 为X*X'的特征向量-
[diagonal,I]=sort(diag(d),'descend'); v=v(:,I);
G=diagonal/sum(diagonal); Gs=0;
for i=1:length(G) Gs(i)=sum(G(1:i)); if Gs(i)> break; end end N=length(Gs)
%v=fliplr(v); %d=rot90(d,2);
翻转( 查看生成 的对 角 阵是由小到大排列的 ~此指令可使其由大到小排列
~fliplr(flipud(d))=rot90(d,2))
%diagonal=diag(d);
spacef=X'*v;
for i=1:Num;
spacef(:,i)=spacef(:,i)/sqrt(diagonal(i)); %
end
%**************************************************************************
将删去的陆地与冰点的填
充值补回
sst_area2=zeros(Num,8400); 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';
%************************************************************************** % 画图
%subplot(2,1,2) figure(1)
x=1:Num; plot(x,timef(:,1),'g'); %ylim([ -80 80 ]);
% xlabel('1980-1992 年 156 个月 ','fontsize',15,'fontname','
ylabel('INDEX','fontsize',12,'fontname',' 黑体 ') 矩阵作左右翻转 矩阵上下翻转后再左右 空间本征函数 timef=X*spacef;
sum_d=sum(diagonal);
count=0;
for i=1:Num;
count=count+diagonal(i);
G1(i)=count/sum_d;
end
for i=1:Num;
G2(i)=diagonal(i)/sum_d;
end
% G1(i) % G2(i) 时间本征函数 是累积方差贡献率 是方差贡献率 将绘图窗口划分为 2*1 个子窗口， 并在第 2 个子窗口中绘图
隶书 ')
set(gca,'xtick',(1:6:Num))
set(gca,'xticklabel',{'1980','','1981','','1982','','1983','','1984','','1985',
'','1986','','1987','','1988','','1989','','1990','','1991','','1992','','1993' }) title(' 北 太 平 洋 第 1 模 态 1980 至 1992 年 SST 时 间 序 列 ', 'color', 'k','fontsize',15,'fontname',' 幼圆 ') grid on hold off % %subplot(2,1,1)
% lon=[:];
% lat=[:];
% m_proj('Equidistant Cylindrical','lat',[ ],'lon',[ ]);
%
m_contourf(lon,lat,rot90(reshape(sst_area3(:,6),140,60)',2),30,'linestyle','non e');
% colorbar
% m_coast('patch',[.95 .95 .95]);
% m_coast('color','k');
% m_grid('linestyle','none','tickdir','out','linewidth',;
% % xlabel('longitude','fontsize',15,'fontname','comic sans ms');
% % ylabel('latitude','fontsize',15,'fontname','comic sans ms');
% title([' 北太平洋第 6 模态填色图 '],'fontsize',15,'fontname','
lon=[:];
lat=[:];
figure(2)
m_proj('lambert','lat',[ ],'lon',[ ]);
m_contourf(lon,lat,rot90(reshape(sst_area3(:,1),140,60)',2)); % colorbar; m_coast('patch',[1 .85 .7]); %m_elev('contourf',[500:500:6000]);
m_grid('box','fancy','tickdir','in');
%colormap(flipud(copper));
% xlabel('longitude','fontsize',15,'fontname',' 幼圆 '); %
ylabel('latitude','fontsize',15,'fontname',' 幼圆 ');
title([' 北太平洋第 1 模态填色图 '],'fontsize',15,'fontname','
colorbar 幼圆 ');
幼圆 ');