matlab EOF程序＋数据. 可使用eof和svd两种方法计算

quinpool

kongfeng0824 发表于 2013-1-18 23:05

                               
登录/注册后可看大图

好像没有经纬度信息呢？

作图时经纬度的信息我没加，主要原因是需要使用m_map，matlab中没有带这个软件包。

quinpool

kongfeng0824 发表于 2013-1-19 02:15

                               
登录/注册后可看大图

上面的这段程序，如果把第二行和第三行的%都去掉，得到的结果和不去掉得到的结果有很大的差别。不知道是怎么 ...

需要求原始资料的协方差，只用第二行或第三行就可以了。

quinpool

tintin 发表于 2012-11-27 15:29

                               
登录/注册后可看大图

请问下载下来就好了吗 还是要安装一下？

下载后，把目录名（例如：放在D:\m_map和D:\m_namebox）加入matlab路径中。具体做法：

建一个addpath.m文件，放到你画图的目录中，
文件内容是：

path(path,'D:\m_map')
path(path,'D:\m_namebox\')
画图程序的开头加上
addpath
就可以了。

quinpool

Heqiwei 发表于 2012-12-6 11:30

                               
登录/注册后可看大图

time=1856+0.5/12:1/12:2003+11.5/12;
帮忙说明下你的数据资料的 时间和格点 谢谢

资料为ZC模式在热带区域逐月的SST距平场，
时间1856年1月-->2003年12月
共1776个月

空间格点数m=25*18;

quinpool

新版本；计算没变，仅加上用m_map画地图的代码。

没有m_map的请自己去下载。然后在计算画图的路径下（例如：D:\）新建一个addpath.m，内容是：   
path(path,'D:\m_map')
path(path,'D:\m_namebox\')

然后把addpath加到eof代码的开头。代码附后。
－－－－－－－－－－－－－－－－－－－






clear; close all;clc;addpath;
      sst= load('m_dt.dat')';
      [m,it]=size(sst);
       ix=25; iy=18;  
      
       time=1856+0.5/12:1/12:2003+11.5/12;
      
       lon=135:5.625:270;
       lat=-17:2:17;
      
    %   for i=1:m
   %   sst(i,:)=deseason(sst(i,:));
    %   end


% (a) get the mean for each month:
        for i=1:m % m is the number of spatial grids (m=25*18=450)
        for k=1:12  % average for each month
        aa(i,k)=mean(sst(i,k:12:it)); % n is the length of sample
        end
        end
%  (b) remove seasonal cycle

        for i=1:m
        for j=1:it
        k=mod(j-1,12)+1;
        sst(i,j)=sst(i,j)-aa(i,k);
        end
        end
        %}
% (4)        calculate covariance matrix
     %s=cov(sst',1);  
     %sst = bsxfun(@minus,sst,sum(sst,1)/n);
   
     %option='eof'
     
     option='svd'
     
     if(option == 'eof' )  
     
     
     
     C=(sst*sst')/it; %remove mean by column (i.e. time dimension, =1776)
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
[E,g]=eig(C);
  E=fliplr(E);   g=rot90(g,2);      gd=diag(g);
  PCA=E'*sst;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  e1=E(:,1);    f1=reshape(e1,ix,iy);
  e2=E(:,2);    f2=reshape(e2,ix,iy);
  e3=E(:,3);    f3=reshape(e3,ix,iy);
   PC1=PCA(1,:);PC2=PCA(2,:);PC3=PCA(3,:);
   

   gd(1:3)
   total=sum(gd);
   first_eof=gd(1)/total*100; %
   second_eof=gd(2)/total*100; %
   third_eof=gd(3)/total*100; %
     %%%%%%%%%%%%%%%%%      
    figure(1);clf;
    subplot(3,2,1)
    m_proj('Equidistant','lat',[-20 20],'lon',[135 270]);   
[c, h] =m_contour(lon,lat,(f1)');
clabel(c, h, 'FontSize', 7);
m_grid('linestyle','none','tickdir','out','linewidth',2);
m_coast('patch','white','edgecolor','k');  
    title (['The First mode of EOF  '  num2str( round(first_eof*10)/10 )   '%'])
  
    subplot(3,2,3)
    m_proj('Equidistant','lat',[-20 20],'lon',[135 270]);   
[c, h] =m_contour(lon,lat,(f2)');
clabel(c, h, 'FontSize', 7);
m_grid('linestyle','none','tickdir','out','linewidth',2);
m_coast('patch','white','edgecolor','k');  
    title (['The second mode of EOF  '   num2str( round(second_eof*10)/10 )   '%'])

    subplot(3,2,5)
   m_proj('Equidistant','lat',[-20 20],'lon',[135 270]);   
[c, h] =m_contour(lon,lat,(f3)');
clabel(c, h, 'FontSize', 7);
m_grid('linestyle','none','tickdir','out','linewidth',2);
m_coast('patch','white','edgecolor','k');  
    title (['The third mode of EOF  '   num2str( round(third_eof*10)/10 )   '%'])
   
   
    subplot(3,2,2)
    plot(time,PC1)
    grid
    title ('PC1 of EOF')
    xlim([1850 2005])

    subplot(3,2,4)
    plot(time,PC2)
    grid
    title ('PC2 of EOF')
    xlim([1850 2005])
   
    subplot(3,2,6)
    plot(time,PC3)
    grid
    title ('PC3 of EOF')
    xlim([1850 2005])
     
     else
         
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  SVD

[u1,v,v1]=svd(sst);
   whos u1 v1
   gd_svd=diag(v);    sigma=gd_svd.^2/it; total_svd=sum(sigma);  
   sigma(1:3)
   first_svd=sigma(1)/total_svd*100;
   second_svd=sigma(2)/total_svd*100;
   third_svd=sigma(3)/total_svd*100;
   
    left1=u1(:,1);left2=u1(:,2);left3=u1(:,3);
    svd1=reshape(left1,ix,iy);
    svd2=(1)*reshape(left2,ix,iy);
    svd3=(1)*reshape(left3,ix,iy);
    %PCA_svd=sst'*u1;
    PCA_svd=v*v1';
    PC1_svd=PCA_svd(1,:);
    PC2_svd=PCA_svd(2,:);
    PC3_svd=PCA_svd(3,:);
%%%%%%%%%%%%%%%%%
    figure(2);clf;
    subplot(3,2,1)
     
m_proj('Equidistant','lat',[-20 20],'lon',[135 270]);   
[c, h] =m_contour(lon,lat,(svd1)');
clabel(c, h, 'FontSize', 7);
m_grid('linestyle','none','tickdir','out','linewidth',2);
%colorbar('position',[0.8 0.59 0.018 0.335],'FontSize',8);
m_coast('patch','white','edgecolor','k');  
   
    title (['The First mode of SVD  '  num2str( round(  first_svd*10)/10 )   '%'])
   

    subplot(3,2,2)
    plot(time,PC1_svd)
    grid
    title ('PC1 of SVD')
    xlim([1850 2005])

   
    subplot(3,2,3)
  m_proj('Equidistant','lat',[-20 20],'lon',[135 270]);   
[c, h] =m_contour(lon,lat,(svd2)');
clabel(c, h, 'FontSize', 7);
m_grid('linestyle','none','tickdir','out','linewidth',2);
%colorbar('position',[0.8 0.59 0.018 0.335],'FontSize',8);
m_coast('patch','white','edgecolor','k');  
    title (['The second mode of SVD  '   num2str( round(  second_svd*10)/10 )   '%'])
   
   
    subplot(3,2,4)
    plot(time,PC2_svd)
    grid
    title ('PC2 of SVD')
    xlabel('year')
    xlim([1850 2005])
   
    subplot(3,2,5)
m_proj('Equidistant','lat',[-20 20],'lon',[135 270]);   
[c, h] =m_contour(lon,lat,(svd3)');
clabel(c, h, 'FontSize', 7);
m_grid('linestyle','none','tickdir','out','linewidth',2);
%colorbar('position',[0.8 0.59 0.018 0.335],'FontSize',8);
m_coast('patch','white','edgecolor','k');  
    title (['The third mode of SVD  '   num2str( round(  third_svd*10)/10 )   '%'])
   
   
    subplot(3,2,6)
    plot(time,PC3_svd)
    grid
    title ('PC3 of SVD')
    xlabel('year')
    xlim([1850 2005])
   
   
     end

quinpool

liang

多谢楼主分享，先收藏备着

陌小默的

楼主好样的~~顶起啦

liujianbao

谢谢楼主，学习学习

willpower-pan

好东西 先留着