爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 2832|回复: 2

[程序设计] 画出的时间序列图有点奇怪,有人帮我看看吗

[复制链接]

新浪微博达人勋

发表于 2016-9-6 11:21:31 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-9-6 11:23:10 | 显示全部楼层
这是画出的sst的时间序列图
捕获.PNG
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-9-6 19:06:57 | 显示全部楼层
谢谢楼主分享~~~~~~~~~
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表