matlab 计算AO指数并海洋气象指数分析都适用

% ao_code_saveclear;clc;close all% 下载的ao_index load('ao_index.mat')Ao_download = ao(:,3); % 自己计算的：load('ao_index_cla.mat')
% figuret=1:888;figureplot(t,AO,'-r')hold onplot(t,Ao_download,'-b')hold on;legend('self','download')
axis([0 888 -4 4]);set(gca,'YTick',-4:1:4,'tickdir','out','yminortick','on');set(gca,'YTicklabel',{'-4','-3','-2','-1','0','1','2','3','4'});set(gca,'XTick',0:60:900,'tickdir','out','xminortick','off');set(gca,'XTicklabel',{'1950','1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010','2015','2020','2025'});set(gca,'linewidth',1,'fontsize',15,'fontname','Times');set(gcf,'color','white','Position',[100,100,1000,500]);export_fig('AO_index123.jpg','-r600')

clc;clear;close all
%% 数据读取file = 'slp.mon.mean.nc';ncdisp('slp.mon.mean.nc');           slp=double(ncread(file,'slp'));    lon=double(ncread(file,'lon'));    lat=double(ncread(file,'lat'));   time=double(ncread(file,'time')); time1 = time./24 + datenum(1800,01,01) ;% 你看注释他从小时出发。time2 = datevec(time1);% 转化为天，因此要除以小时24ln1 = find(time2(:,1)>=1950&time2(:,1)<=2023);%% 选区域ln2 = find(lat >=20);lat = lat(ln2);slp_TP=slp(:,ln2,ln1);

slp_TP_12=reshape(slp_TP,144,29,12,74);[lon_TP,lat_TP,month_TP,year_TP]=size(slp_TP_12);

average=squeeze(nanmean(slp_TP_12(:,:,:,1:year_TP),4));  % 74年
% 算距平slp_TP_anomaly_12=slp_TP_12-repmat(average,[1,1,1,year_TP]);
% 降维，排列成m×n的形式（m：站点；n：时间）slp_TP_anomaly=reshape(slp_TP_anomaly_12,lon_TP,lat_TP,month_TP*year_TP);X_NaN=reshape(slp_TP_anomaly,lon_TP*lat_TP,month_TP*year_TP);
X_logical=isnan(X_NaN);  
i=0;for j=1:lon_TP*lat_TP    if sum(X_logical(j,:))==0        i=i+1;        X(i,:)=X_NaN(j,:);    endend
X=X/sqrt(month_TP*year_TP); 
%[X_x,X_y]=size(X);if X_x<=X_y             % 不考虑时空转换
    S=X*X';             % 计算协方差    [v,d]=eig(S);       % 使用eig函数求特征向量v和特征值d
    % 使用sort函数将特征值、特征向量按升序排序    [~,ind]=sort(diag(d));    Ds=d(ind,ind);    Vs=v(:,ind);
    V=fliplr(Vs);  % 排序后的特征向量V    D=rot90(Ds,2); % 排序后的特征值D
else                     % 考虑时空转换    S=X'*X;              % 计算协方差    [v,d]=eig(S);        % 使用eig函数求特征向量v和特征值d
    % 使用sort函数将特征值、特征向量按升序排序    [~,ind]=sort(diag(d));    Ds=d(ind,ind);    Vs=v(:,ind);
    VR=fliplr(Vs);  % 排序后的特征向量V    D=rot90(Ds,2);  % 排序后的特征值D
    % 求X*X'的特征值（《现代气候统计诊断与预测技术》109页）    VR=X*VR;    for i=1:length(D)        V(:,i)=VR(:,i)/sqrt(D(i,i));    endend
% 时间系数T=V'*X*sqrt(month_TP*year_TP);
% 特征向量的方差贡献diagonal=diag(D);total=0;for i=1:length(diagonal)    R(i)=diagonal(i)/sum(diagonal);  % 方差贡献率    total=total+diagonal(i);    G(i)=total/sum(diagonal);        % 累积方差贡献率end
% 标准化处理for i=1:length(diagonal)    V(:,i)=V(:,i)*sqrt(D(i,i));    % 特征向量乘以特征值的平方根    T(i,:)=T(i,:)/sqrt(D(i,i));    % 时间系数除以特征值的平方根end
%% 把删去的NaN数据补回[a,b]=size(V);
V_NaN=NaN(lon_TP*lat_TP,b);j=1;for i=1:lon_TP*lat_TP    if sum(X_logical(i,:))==0        V_NaN(i,:)=V(j,:);        j=j+1;    endend
%% 升维V_NaN=reshape(V_NaN,lon_TP,lat_TP,b);
% 保存数据AO=T_1;save ao_index_cla.mat AO