- 积分
- 26784
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-4-25
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 伽蓝鸟 于 2020-2-7 18:37 编辑
有朋友问我怎么用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/
|
评分
-
查看全部评分
|