[数据处理] 免费分享的Matlab滑动EOF代码和文献资料

[复制链接]

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

$ a" v4 B  C1 n) i' }
[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;
0 i+ K8 u5 @  |" n6 }1 `/ R
( o2 q. y) N& X  [& A( I' ^0 k

f40b4145395494c2720cbc1934d5f3f2.png


  _" R2 u2 s+ M" l

更多内容:

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


/ ?- Z6 s; G- R2 `0 S
* |* \# e& [; k! [
# ?* U4 _: i" v* `! H3 q8 E
回复

举报 使用道具

相关帖子

全部回帖
非常感谢楼主的分享,本科小弟表示很受用2 H! Z8 A% Z* I" c  p, L
发表于 2022-8-8 12:48:09

举报 回复 使用道具

懒得打字?点击右侧快捷回复 【吾爱海洋论坛发文有奖】
您需要登录后才可以回帖 登录 | 立即注册
象山港
活跃在2024-12-1
快速回复 返回顶部 返回列表