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

quinpool

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

option='eof'
用eof方法计算

option='svd'  －－>用svd方法计算。

两种方法结果完全一致。


文件名:  eof.rar
下载地址:  http://www.rayfile.com/files/14b103d7-7bf7-11e1-ab61-0015c55db73d/

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

附件看不到。重新贴上来。

clear; close all;clc;
      sst= load('m_dt.dat')';
      [m,it]=size(sst);
       ix=25; iy=18;  
       time=1856+0.5/12:1/12:2003+11.5/12;
      
      
      
    %   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)
    contour(f1')
    title (['The First mode of EOF  '  num2str( round(first_eof*10)/10 )   '%'])
  
    subplot(3,2,3)
    contour(f2')
    title (['The second mode of EOF  '   num2str( round(second_eof*10)/10 )   '%'])

    subplot(3,2,5)
    contour(f3')
    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)
    contour(svd1')
    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)
    contour(svd2')
    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)
    contour(svd3')
    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

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
就可以了。

言深深

不带这样的，楼主发两次，一个财神卡一个幸运星······

tingting

先备着  指不定就要用呢 谢谢分享

大熊

问下啊，你那个测试数据画出来标的是格点数啊，它的经纬度是什么啊？

。７０％淡

{:eb302:}好东西 先留着

quinpool

大熊 发表于 2012-4-21 14:14

                               
登录/注册后可看大图

问下啊，你那个测试数据画出来标的是格点数啊，它的经纬度是什么啊？

经纬度信息没画上去，主要考虑matlab画地图要调用地图map软件包。这个小软件包需要另外下载。
也许新版本的matlab可以直接画地图，我现在还是用调用老的地图包的方式。

自己画图时设上自己资料网格所对应的经纬度。

latesunshine

谢谢分享~
正在学习中。请问这个实验数据是什么数据？可以简要说明下吗？谢谢！

威震寰宇_佳

留着备用 多谢分享~~

陌小默的

哈哈给力啊