EOF MATLAB

liuxiuqing

1. clc;clear;close all
   
2. % 参考书籍
   
3. % 黄嘉佑《气象统计分析与预报方法》（第四版）
   
4. % 魏凤英《现代气候统计诊断与预测技术》（第2版）
   
5. % matlab eig 函数：https://ww2.mathworks.cn/help/matlab/ref/eig.html
   
6. % matlab M_Map 工具箱：https://www.eoas.ubc.ca/~rich/map.html
   
7. % matlab M_Map 工具箱工具箱的安装：http://blog.sina.com.cn/s/blog_8fc890a20102v6pm.html
   
8. %% 数据读取
   
9. % ncdisp('sst-195001-201712.nc');             % NOAA Extended Reconstructed SST V4
   
10. sst=ncread('sst-195001-201712.nc','sst');     % 海温
    
11. lon=ncread('sst-195001-201712.nc','lon');     % 经度，分辨率：2°
    
12. lat=ncread('sst-195001-201712.nc','lat');     % 纬度，分辨率：2°
    
13. time=ncread('sst-195001-201712.nc','time');   % 时间，月平均，68年，12个月
    
14. %% 选区域
    
15. % 框选热带太平洋地区(Tropical Pacific：120°E-70°W；30°S-30°N）
    
16. sst_TP=sst(61:146,30:60,:);
    
17. % 把数据展开为四维（经度、纬度、月、年）
    
18. sst_TP_12=reshape(sst_TP,86,31,12,68);
    
19. [lon_TP,lat_TP,month_TP,year_TP]=size(sst_TP_12);
    
20. %% 资料预处理
    
21. % 求平均(每个格点不同月份的年平均)
    
22. average=squeeze(nanmean(sst_TP_12(:,:,:,1:year_TP),4));  % 68年
    
23. % 算距平（减去对应月份的平均值）
    
24. sst_TP_anomaly_12=sst_TP_12-repmat(average,[1,1,1,year_TP]);
    
25. % 降维，排列成m×n的形式（m：站点；n：时间）
    
26. sst_TP_anomaly=reshape(sst_TP_anomaly_12,lon_TP,lat_TP,month_TP*year_TP);
    
27. X_NaN=reshape(sst_TP_anomaly,lon_TP*lat_TP,month_TP*year_TP);
    
28. % 去除NaN数据（NaN为陆地，没有海温）
    
29. X_logical=isnan(X_NaN);  % 逻辑判断
    
30. i=0;
    
31. for j=1:lon_TP*lat_TP
    
32.     if sum(X_logical(j,:))==0
    
33.        i=i+1;
    
34.        X(i,:)=X_NaN(j,:);
    
35.     end
    
36. end
    
37. X=X/sqrt(month_TP*year_TP);  % 《气象统计分析与预报方法》142页
    
38. %% EOF运算
    
39. [X_x,X_y]=size(X);
    
40. if X_x<=X_y             % 不考虑时空转换
    
41.     S=X*X';             % 计算协方差
    
42.     [v,d]=eig(S);       % 使用eig函数求特征向量v和特征值d
    
43.     % 使用sort函数将特征值、特征向量按升序排序
    
44.     [~,ind]=sort(diag(d));
    
45.     Ds=d(ind,ind);
    
46.     Vs=v(:,ind);
    
47.     V=fliplr(Vs);  % 排序后的特征向量V
    
48.     D=rot90(Ds,2); % 排序后的特征值D
    
49. else                     % 考虑时空转换
    
50.     S=X'*X;              % 计算协方差
    
51.     [v,d]=eig(S);        % 使用eig函数求特征向量v和特征值d
    
52.     % 使用sort函数将特征值、特征向量按升序排序
    
53.     [~,ind]=sort(diag(d));
    
54.     Ds=d(ind,ind);
    
55.     Vs=v(:,ind);
    
56.     VR=fliplr(Vs);  % 排序后的特征向量V
    
57.     D=rot90(Ds,2);  % 排序后的特征值D
    
58.     % 求X*X'的特征值（《现代气候统计诊断与预测技术》109页）
    
59.     VR=X*VR;        
    
60.     for i=1:length(D)
    
61.         V(:,i)=VR(:,i)/sqrt(D(i,i));
    
62.     end
    
63. end
    
64. % 时间系数
    
65. T=V'*X*sqrt(month_TP*year_TP);
    
66. % 特征向量的方差贡献
    
67. diagonal=diag(D);
    
68. total=0;
    
69. for i=1:length(diagonal)
    
70.     R(i)=diagonal(i)/sum(diagonal);  % 方差贡献率
    
71.     total=total+diagonal(i);
    
72.     G(i)=total/sum(diagonal);        % 累积方差贡献率
    
73. end
    
74. % 标准化处理
    
75. for i=1:length(diagonal)
    
76.     V(:,i)=V(:,i)*sqrt(D(i,i));    % 特征向量乘以特征值的平方根
    
77.     T(i,:)=T(i,:)/sqrt(D(i,i));    % 时间系数除以特征值的平方根
    
78. end
    
79. %% 把删去的NaN数据补回
    
80. [a,b]=size(V);
    
81. 
    
82. V_NaN=NaN(lon_TP*lat_TP,b);
    
83. j=1;
    
84. for i=1:lon_TP*lat_TP
    
85.     if sum(X_logical(i,:))==0
    
86.        V_NaN(i,:)=V(j,:);
    
87.        j=j+1;
    
88.     end
    
89. end
    
90. %% 升维
    
91. V_NaN=reshape(V_NaN,lon_TP,lat_TP,b);
    
92. %% 显著性检验
    
93. figure
    
94. err=sqrt(2/(month_TP*year_TP))*diagonal/sum(diagonal); % 特征值的取样误差（North et al.,1982）
    
95. dot=1:4;
    
96. R_1_4=R(1:4)*100;
    
97. err=err(1:4)'*100;
    
98. plot(dot,R_1_4,'.','Markersize',12);
    
99. hold on
    
100. errorbar(dot,R_1_4,err,'LineStyle','none');
     
101. axis([0 5 0 50]);
     
102. set(gca,'YTick',0:10:50,'tickdir','out','yminortick','off');
     
103. set(gca,'YTicklabel',{'0','10','20','30','40','50'});
     
104. set(gca,'XTick',0:1:5,'tickdir','out','xminortick','off');
     
105. set(gca,'XTicklabel',{'0','1','2','3','4','5'});
     
106. set(gca,'linewidth',1,'fontsize',15,'fontname','Times');
     
107. set(gcf,'color','white','Position',[100,100,1000,500]);
     
108. title('Percentage variance (%)','fontsize',15,'color','b','position',[2.5,52]);
     
109. %% 第一模态画图
     
110. load colormap.mat
     
111. figure ('colormap',colormap)
     
112. % 空间分布
     
113. subplot(2,1,1);
     
114. V_1=V_NaN(:,:,1);  
     
115. m_proj('miller','lon',[120 290],'lat',[-30 30]);
     
116. [X1,Y1]=meshgrid(120:2:290,-30:2:30);
     
117. m_contourf(X1,Y1,rot90(V_1),'linestyle','none');
     
118. hold on
     
119. m_contour(X1,Y1,rot90(V_1),[-0.1 0.1 0.4 0.8],'linecolor','k','ShowText','on');
     
120. m_coast('patch',[.85 .85 .85 ],'edgecolor',[0.6 0.6 0.6]);
     
121. m_grid('linestyle','none','xtick',[140 180 220 260],'ytick',[30 15 0 -15 -30],'fontsize',14,'tickdir','out');
     
122. caxis([-1,1])
     
123. ax=colorbar;
     
124. set(ax,'linewidth',0.5,'fontsize',12,'edgecolor','k','YTick',[-1,-0.5,0,0.5,1],'YTickLabel',{'-1','-0.5','0','0.5','1'})
     
125. set(gcf,'color','white','Position',[100,100,1000,500]);
     
126. title(['(a) EOF1 = ' num2str(R(1:1)*100,4) '%'],'fontsize',15,'color','b');
     
127. % 时间系数
     
128. subplot(2,1,2);
     
129. T_1=T(1,:);
     
130. t=[1:month_TP*year_TP];
     
131. plot(t,T_1,'k','linewidth',1);
     
132. hold on;
     
133. T_0=zeros(1,month_TP*year_TP);
     
134. plot(t,T_0,'k--','linewidth',1);
     
135. axis([0 840 -4 4]);
     
136. set(gca,'YTick',-4:1:4,'tickdir','out','yminortick','on');
     
137. set(gca,'YTicklabel',{'-4','-3','-2','-1','0','1','2','3','4'});
     
138. set(gca,'XTick',0:60:840,'tickdir','out','xminortick','off');
     
139. set(gca,'XTicklabel',{'1950','1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010','2015','2020'});
     
140. set(gca,'linewidth',1,'fontsize',15,'fontname','Times');
     
141. set(gcf,'color','white','Position',[100,100,1000,500]);
     
142. title('(b) PC1','fontsize',15,'color','b');
     
143. %% 第二模态画图
     
144. load colormap.mat
     
145. figure ('colormap',colormap)
     
146. % 空间分布
     
147. subplot(2,1,1);
     
148. V_2=V_NaN(:,:,2);  
     
149. m_proj('miller','lon',[120 290],'lat',[-30 30]);
     
150. [X1,Y1]=meshgrid(120:2:290,-30:2:30);
     
151. m_contourf(X1,Y1,rot90(V_2),'linestyle','none');
     
152. hold on
     
153. m_contour(X1,Y1,rot90(V_2),[-0.1 0.1 0.4 0.8],'linecolor','k','ShowText','on');
     
154. m_coast('patch',[.85 .85 .85 ],'edgecolor',[0.6 0.6 0.6]);
     
155. m_grid('linestyle','none','xtick',[140 180 220 260],'ytick',[30 15 0 -15 -30],'fontsize',14,'tickdir','out');
     
156. caxis([-1,1])
     
157. ax=colorbar;
     
158. set(ax,'linewidth',0.5,'fontsize',12,'edgecolor','k','YTick',[-1,-0.5,0,0.5,1],'YTickLabel',{'-1','-0.5','0','0.5','1'})
     
159. set(gcf,'color','white','Position',[100,100,1000,500]);
     
160. title(['(a) EOF2 = ' num2str(R(2:2)*100,4) '%'],'fontsize',15,'color','b');
     
161. % 时间系数
     
162. subplot(2,1,2);
     
163. T_2=T(2,:);
     
164. t=[1:month_TP*year_TP];
     
165. plot(t,T_2,'k','linewidth',1);
     
166. hold on;
     
167. T_0=zeros(1,month_TP*year_TP);
     
168. plot(t,T_0,'k--','linewidth',1);
     
169. axis([0 840 -4 4]);
     
170. set(gca,'YTick',-4:1:4,'tickdir','out','yminortick','on');
     
171. set(gca,'YTicklabel',{'-4','-3','-2','-1','0','1','2','3','4'});
     
172. set(gca,'XTick',0:60:840,'tickdir','out','xminortick','off');
     
173. set(gca,'XTicklabel',{'1950','1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010','2015','2020'});
     
174. set(gca,'linewidth',1,'fontsize',15,'fontname','Times');
     
175. set(gcf,'color','white','Position',[100,100,1000,500]);
     
176. title('(b) PC2','fontsize',15,'color','b');

复制代码

沐然丶

你这复制粘贴别人代码还不标出处不太合适吧

吃肉的菜鸡

您好，我想请问一下，做EOF分析是不是一定要要求数据是正态分布的，如果不是正态分布的数据可以做EOF嘛？