免费分享的Matlab滑动EOF代码和文献资料

很多小伙伴们应该都用过EOF分析，提取主导模态的一种有效方法。滑动啥啥啥也是常用的统计方法，比如滑动平均，滑动相关，滑动方差。可以拿来看一些信号的年代际变化。何超师兄有篇CD，讲中国东部降水主导模态的年代际变化，从三极子型变成偶极子型(He et al. 2017 CD)。其实之前也有很多工作讨论过中国东部降水主导模态的年代际变化，例如Ye et al.(2012 AAS)；或者讨论主导模态的年代际变化，比如Kwon et al.(2005 GRL)。但是他们都是对于两个时段，分别做EOF，然后来说主导模态有所不同。滑动EOF就是滑动着做EOF，能够看得出来啥时候发现了模态的转变。蛮有意思的一种方法，但用的人可能比较少？正好有朋友问我这个的代码怎么写，就重复了师兄文章里的图。代码和文献一并附上。

+ s) T8 u# N4 R- S, g, T

[C] 纯文本查看 复制代码

clc,clear,close all;fclose all;
%% 读入原始的GPCP数据
ncdisp('F:/RMatlab/data20211219/precip.mon.mean.nc');
x=ncread('F:/RMatlab/data20211219/precip.mon.mean.nc','lon');
x=double(x);
y=ncread('F:/RMatlab/data20211219/precip.mon.mean.nc','lat');
y=double(y);
gpcp=ncread('F:/RMatlab/data20211219/precip.mon.mean.nc','precip');
%%
gpcp(:,:,517:end)=[];                           % % 把末尾数据（2021年）给删掉
gpcp=reshape(gpcp,length(x),length(y),12,[]);   % % 经度*纬度*月份*年份
gpcp=gpcp(x>=107 & x<=122 , y>=20 & y<=43, : , :);  % % 107-122E, 20-43N
gpcp=squeeze(mean(mean(gpcp(:,:,6:8,:),3),1));      % % JJA平均
for ii=1:9
    gpcp(ii,:)=zscore(squeeze(gpcp(ii,:)));         % % 标准化变量
end
runEOF1=nan(9,42);           % % 最后画图的数据
for ii=8:35
    XX=gpcp(:,ii-7:ii+7);    % % 纬度*15年滑动
    [v,d]=eig(XX*XX');
    lam=real(diag(d));
    [B,idx]=sort(lam,'descend');
    VV=v(:,idx(1:5));
    PC=VV'*XX;
    PC=PC';                  % % EOF的计算过程
    PC1=PC(:,1);             % % PC1
    EOF1=VV(:,1)*std(PC1);   % % EOF1
    runEOF1(:,ii)=EOF1;      % % 滑动EOF1
end
%% 但是这里有个问题，滑动EOF的符号可能突变，画图就不好看了
for ii=8:34
    if dot(runEOF1(:,ii),runEOF1(:,ii+1))<0
        runEOF1(:,ii+1)=-1*runEOF1(:,ii+1);
        % % 如果后一年符号有突变，就乘以个-1，换回来
    end
end
%% 画图看看
figure,hold on;
contourf(1979:2020,21.25:2.5:41.25,runEOF1,...
    -1.5:0.01:1.5,'linestyle','none');
colormap('jet');colorbar;caxis([-1 1]);
axis([1986 2013 21.25 41.25]);
c=colorbar('eastoutside','ticklength',0);
caxis([-1,1])
% 调用ncl色带
load('colorbar-mat/BkBlAqGrYeOrReViWh200_r.mat');
colormap(BkBlAqGrYeOrReViWh200_r);
ax = gca;
axpos = ax.Position;
c.Position(3) = 0.5*c.Position(3);
ax.Position = axpos;
cbarrow;

9 N* E5 M* l, {- f

  O4 m2 j  B! a2 A4 J+ Z( b

更多内容：

游客，如果您要查看本帖隐藏内容请回复

, T- m- c8 W& G% T+ K0 _
+ O3 j: ?3 B% W& Y8 T9 I$ m: J* X

非常感谢楼主的分享，本科小弟表示很受用