- 积分
- 276
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-8-5
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
clear
sst=zeros(360,180,240);
pre=zeros(144,72,240); %截取时间为1989年1月到2009年12月,共20年240个月
lon1=ncread('E:data\sst.mon.mean.nc','lon');
lat1=ncread('E:data\sst.mon.mean.nc','lat');
time1=ncread('E:data\sst.mon.mean.nc','time');
sstdata=ncread('E:data\sst.mon.mean.nc','sst');
lon2=ncread('E:data\precip.mon.mean (3).nc','lon');
lat2=ncread('E:data\precip.mon.mean (3).nc','lat');
time2=ncread('E:data\precip.mon.mean (3).nc','time');
pre=ncread('E:data\precip.mon.mean (3).nc','precip');
%%求sst的eof
for i=1:240
s=sstdata(:,:,i+100);
s(s>100)=nan;
sst(:,:,i)=s;
end
sst2=zeros(240,360*180); %%数组降维
for i=1:240
squ1=sst(:,:,i);
sst2(i,:)=reshape(squ1,1,64800);
end
pre_nan=isnan(sst2);
i=0;
sst3=[];
for j=1:64800 %对于每一个空间点
if sum(pre_nan(:,j))==0; %时间上求和等于0
i=i+1; %时间求和为0时,i++
sst3(:,i)=sst2(:,j); %筛选掉一些空间点,这些空间点是陆地,无值
end
end
%下面开始基于月尺度的EOF分解
Y1=zeros(size(sst3));
Y1=sst3-repmat(mean(sst3,1),size(sst3,1),1); %%Y1'为预处理后的数据距平矩阵
X1=Y1*Y1'/43798;
[V1,D1]=eig(X1);
V1=fliplr(V1); %矩阵作左右翻转
D1=rot90(D1,2); %矩阵上下翻转后再左右翻转(查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr(flipud(d))=rot90(d,2))
diagonal1=diag(D1); %d是特征值的对角矩阵,v是特征向量
%%得到D是特征值的对角矩阵,V是特征向量
g=[];
sum_d1=sum(diagonal1);
for i=1:240;
g(i)=diagonal1(i)/sum_d1; %方差贡献率矩阵
end
spacef1=Y1'*V1; %主成分
for i=1:240;
spacef1(:,i)=spacef1(:,i)/sqrt(diagonal1(i)); % 空间本征函数,每一列都是特征向量
end
timef1=Y1*spacef1; %时间本征函数
%将删去的填充值补回
sst4=zeros(240,64800);
sst4(:,:)=nan;
i=1;
spacef2=spacef1';
for j=1:64800 %对于每一个空间点
if sum(pre_nan(:,j))==0;
sst4(:,j)=spacef2(:,i);
i=i+1;
end
end
sst5=sst4';
%%——————————————————————————开始作图
figure(1); %第一模态时间序列图
x=1:1:240;
plot(x,timef1(:,1),'m','linewidth',2);
%xlabel('1961-2010年50年','fontsize',15,'fontname','隶书')
ylabel('INDEX','fontsize',12,'fontname','黑体')
title('年平均sst第1模态时间序列', 'color', 'k','fontsize',15,'fontname','幼圆')
grid on
grid minor;
hold off
saveas(gcf,'年平均sst第1模态时间序列','jpg');
%————————————————————————EOF图
fig=zeros(360,180);
fig=reshape(sst5(:,1),360,180);
figure(2)
[LON1,LAT1]=meshgrid(lon1,lat1);
[cs,h]=contourf(LON1,LAT1,fig',20);
z = clabel(cs, h, 'FontSize', 10, 'Color', 'r', 'Rotation', 0);
fig1=zeros(360,180);
fig1=reshape(sst5(:,2),360,180);
figure(5)
[LON1,LAT1]=meshgrid(lon1,lat1);
[cs,h]=contourf(LON1,LAT1,fig1',20);
z = clabel(cs, h, 'FontSize', 10, 'Color', 'r', 'Rotation', 0);
for i=1:240
s=pre(:,:,i+120);
s(s<0)=nan;
pre(:,:,i)=s;
end
%----------------------------------------------------
pre2=zeros(240,144*72); %%数组降维
for i=1:240
squ=pre(:,:,i);
pre2(i,:)=reshape(squ,1,144*72);
end
pre_nan=isnan(pre2);
i=0;
pre3=[];
for j=1:144*72 %对于每一个空间点
if sum(pre_nan(:,j))==0; %时间上求和等于0
i=i+1; %时间求和为0时,i++
pre3(:,i)=pre2(:,j); %筛选掉一些空间点,这些空间点是海洋,无值
end
end
%下面开始基于年尺度的EOF分解
Y2=zeros(size(pre3));
Y2=pre3-repmat(mean(pre3,1),size(pre3,1),1);
X2=Y2*Y2'/9556; %这里做了时空转换,因为空间序列数远远大于时间序列数
[V2,D2]=eig(X2);
V2=fliplr(V2); %矩阵作左右翻转
D2=rot90(D2,2); %矩阵作左右 %矩阵上下翻转后再左右翻转(查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr(flipud(d))=rot90(d,2))
diagonal2=diag(D2); %d是特征值的对角矩阵,v是特征向量
%%得到D是特征值的对角矩阵,V是特征向量
g2=[];
sum_d=sum(diagonal2);
for i=1:240
g2(i)=diagonal2(i)/sum_d; %方差贡献率矩阵
end
spacef_1=Y2'*V2; %主成分
for i=1:240;
spacef_1(:,i)=spacef_1(:,i)/sqrt(diagonal2(i)); % 空间本征函数,每一列都是特征向量
end
timef_1=spacef_1'*Y2'; %时间本征函数
%将删去的填充值补回
pre4=zeros(240,144*72);
pre4(:,:)=nan;
i=1;
spacef_2=spacef_1';
for j=1:144*72 %对于每一个空间点
if sum(pre_nan(:,j))==0;
pre4(:,j)=spacef_2(:,i);
i=i+1;
end
end
pre5=pre4';
%%开始作图
figure(3); %%第一模态图
x=1:1:240;
plot(x,timef_1(1,:),'m','linewidth',2);
%xlabel('1961-2010年50年','fontsize',15,'fontname','隶书')
ylabel('INDEX','fontsize',12,'fontname','黑体')
set(gca ,'xtick',(1:3:240))
set(gca,'xticklabel',{'1989','','1994','','1999','','2004','','2009',''})
title('年平均precip第1模态时间序列', 'color', 'k','fontsize',15,'fontname','幼圆')
grid on
grid minor;
hold off
saveas(gcf,'年平均precip第1模态时间序列','jpg');
%----------eof
pres=zeros(144,72);
pres=reshape(pre5(:,1),144,72);
figure(4)
contourf(flipud(pres'));
hold off
|
|