南信大气象统计实习8EOF分析

耿瑞满

%shix8 EOF分析
%有需要加QQ:3588430252或者发邮箱：3588430252@qq.com
%lat=-27.5:5:27.5N,12
%lon=120:10:290E,18
%time=i第i个时次
%先对18*12进行处理，再最后转置

clc;clear;
ncdisp("SSTPX.nc");
sst_ano = ncread('SSTPX.nc','s');

%为方便处理将NAN改成0并改三维成二维
%sst_X为二维距平-时间场
%sst_X=sst_V*sst_Z
%sst_V空间模态，sst_Z时间系数
%expvar方差贡献
sst_ano(isnan(sst_ano))=0;
[lon,lat,time]=size(sst_ano);
sst_X=zeros(lon*lat,time);
for j=1:time
    sst_X(:,j)=reshape(sst_ano(:,:,j),[216,1]);
end
sst_A=sst_X*sst_X';
[sst_V,sst_D]=sort_by_lambda(sst_A);
sst_Z=sst_V'*sst_X;
expvar=sst_D/sum(sst_D);


%North检验,主要是前四个模态通过North检验
%sig为该空间模态是否通过North检验
e=sst_D*(2/time).^0.5;
sst_D(end+1)=0;
sig=zeros(lon*lat,0);
for k=1:lon*lat
    if((sst_D(k)-sst_D(k+1))>=e(k))
        sig(k)=1;
    end
end
sig=sig';

%% 第一空间模态以及其对应时间系数
choice=2;
v1=sst_V(:,choice);
z1=sst_Z(choice,:);
v2=reshape(v1,[18,12]);
v2=v2';


figure(1);
subplot(2,1,1)
set(figure(1),'color','white')
addpath 'D:\Codes_Projects\MATLAB_codes\othercolor'
addpath 'D:\MATLAB\toolbox\m_map'

min_lon=120;
max_lon=290;
step_lon=10;
min_lat=-27.5;
max_lat=27.5;
step_lat=5;
lon2=min_lon:step_lon:max_lon;
lat2=min_lat:step_lat:max_lat;

m_proj('Equidistant Cylindrical','long',[min_lon max_lon],'lat',[min_lat max_lat]);%选中范围

% 使用 griddata 函数进行二维插值  
lon_grid = linspace(min_lon, max_lon, 200);  
lat_grid = linspace(min_lat, max_lat, 150);
[lon2,lat2]=meshgrid(lon2,lat2);
[lon_grid, lat_grid] = meshgrid(lon_grid, lat_grid );
v2_grid = griddata(lon2, lat2, v2,lon_grid, lat_grid,'v4');


%画填色图加通过显著性检验用阴影,
% x,y归化成图上坐标,Color(4)为透明度；
[c,h]=m_contourf(lon_grid,lat_grid,v2_grid);
set(h,'LineColor','none');

%色标
cbar=colorbar('eastoutside','ticklength',0);
t=get(cbar,'YTickLabel');
t=strcat(t,'K');
colormap(othercolor('RdBu10'))

%绘制范围内的地图
ma=shaperead('D:\Maps\world\World_countries.shp');
m_line([ma(:).X],[ma(:).Y],'color','k');
m_grid('linest','none','fontname','Times','fontsize',8,'linewidth',1);
title(sprintf("(a)eof%d  Variance=%.2f%%",choice,expvar(choice)*100),'FontSize',11);


subplot(2,1,2)
time=datetime(1948,1,1)+calmonths(0:515);
bar(time(z1>0), z1(z1>0),'b');
hold on
bar(time(z1<0), z1(z1<0),'r');
title(sprintf("(b)pc%d",choice))
saveas(gcf,sprintf('eof%d.png',choice))