MATLAB做REOF

伽蓝鸟

有朋友问我怎么用MATLAB做REOF，今天尝试了下。之前算EOF的时候，都是自己码的代码，发现内置函数pca就够了。而旋转也有内置函数，rotatefactors。重复了一下Hannachi et al. (2007 IJC)里的图，看着还行。

代码如下
clc,clear,close all;fclose all;
%%
x=0:2.5:357.5;                          % % 经度范围
y=20:2.5:87.5;                          % % 纬度范围
fid=fopen('D:\winter1920\slp.dat');     % % 用GrADS写的1948-2000年NCEP/NCAR冬季的海平面气压数据
slp=fread(fid,'single');fclose(fid);
slp=reshape(slp,length(x),length(y),[]);    % % 读入数据
%%
% % 扣去平均值，计算距平
for ii=1:length(x)
    for jj=1:length(y)
        slp(ii,jj,:)=slp(ii,jj,:)-mean(squeeze(slp(ii,jj,:)));
    end
end
% % 纬度加权
for jj=1:length(y)
    slp(:,jj,:)=slp(:,jj,:)*sqrt(cosd(y(jj)));
end
%%
XX=reshape(slp,length(x)*length(y),[]);             % % 把数据重新整理成二维的
[EOF,PC,lam] = pca(XX');                                   % % MATLAB自带的函数做EOF
%%
% % 用rotatefactors做旋转
[REOF,G]=rotatefactors(EOF(:,1:6),'method','varimax','Normalize','on');
RPC=PC(:,1:6)*G;
%%
fid=fopen('D:\winter1920\REOF.dat','w');
fwrite(fid,reshape(REOF,[],1),'single');fclose(fid);
fid=fopen('D:\winter1920\REOF.ctl','w');
fprintf(fid,'%s\n','dset D:\winter1920\REOF.dat');
fprintf(fid,'%s\n','title variable');
fprintf(fid,'%s\n','undef -9.99e+33');
fprintf(fid,'%s\n','xdef 144 linear 0 2.5');
fprintf(fid,'%s\n','ydef 28 linear 20 2.5');
fprintf(fid,'%s\n','zdef 1 levels 1000 1');
fprintf(fid,'%s\n','tdef 6 linear 00Z01Jan1979 1yr');
fprintf(fid,'%s\n','vars 1');
fprintf(fid,'%s\n','REOF=>REOF 1 t,z,y,x variable');
fprintf(fid,'%s\n','endvars');
fclose(fid);

文献是这篇：Hannachi A, Jolliffe I T, Stephenson D B. Empirical orthogonal functions and related techniques in atmospheric science: A review[J]. International Journal of Climatology, 2007, 27(9): 1119-1152.

师姐刚刚给了我一个网址，里头做REOF的代码给出的结果也是类似的。网址是websites.pmc.ucsc.edu/~dmk/notes/EOFs/

更大华

{:5_213:}{:5_213:}{:5_213:}

LiShuangshuang

作者棒棒哒

rceclx

谢谢分享！
请问 'Normalize','on' 的作用是什么？

伽蓝鸟

rceclx 发表于 2020-2-11 18:01
谢谢分享！
请问 'Normalize','on' 的作用是什么？

关于标准化，我也不是很懂诶。Wilks那本书上写了好多。见他书上第12.5.3，正交旋转对初始特征向量尺度的敏感性。

rceclx

伽蓝鸟 发表于 2020-2-11 20:27
关于标准化，我也不是很懂诶。Wilks那本书上写了好多。见他书上第12.5.3，正交旋转对初始特征向量尺度的 ...

好的，我再查一下这个函数的参数。

nkuhyl

感谢分享，学到了

shishuilingshan

请教楼主啊，如果手上是站点资料 需要先插值到格点再继续做reof吗？还是可以直接做啊 希望得到回复 谢谢~

伽蓝鸟

shishuilingshan 发表于 2020-2-20 09:48
请教楼主啊，如果手上是站点资料 需要先插值到格点再继续做reof吗？还是可以直接做啊 希望得到回复 谢谢~

我建议是直接算REOF，画图的时候再插值（个人意见仅供参考啊）

shishuilingshan

伽蓝鸟 发表于 2020-2-21 08:08
我建议是直接算REOF，画图的时候再插值（个人意见仅供参考啊）

感谢楼主回复哈~这几天一直在站里看reof的帖子，似乎都是先用griddata插值然后算reof，想继续问下怎么用站点资料直接做reof啊？