Matlab 气象局网格；CO2；各点相关；气压海拔；EDW；SVD

小莹子

1.批量处理中国气象数据网格点数据
clc;
clear;

datadir='/data/Mete/';
file=dir([datadir,'*.txt']);
k=length(file);
% a=textread('SURF_CLI_CHN_TEM_MON_GRID_0.5-MEAN-201505.txt','%s','delimiter','\n');
% a(cellfun(@isempty,a))=[]  %删除空行
for i=1:k
    a=textread(file(i).name,'%s','delimiter','\n');
    A=a(7:78);            
    B=cell2mat(deblank(A));  
    C=(str2num(B))';   
    T=double(C);
    Lat=linspace(53.75,18.25,72);
    lat=double(Lat);
    Lon=linspace(72.25,135.75,128);
    lon=double(Lon);
    [X,Y] = meshgrid(lat,lon);
    lat1 = xlsread('/data/站点/1979-2018.xlsx','Sheet2','e2:e87');           %站点纬度
    lon1 = xlsread('/data/站点/1979-2018.xlsx','Sheet2','d2:d87');             %站点经度
    Tt= interp2(X,Y,T,lat1,lon1,'linear');  %linear 双线性插值法
    temp(:,:,i)=Tt(:,:);
end
tt=reshape(temp,[86 480]);
Ttb=reshape(tt(:,13:468),[86 12 38]);
T0=reshape(mean(Ttb(:,1:12,:),2),[86 38]);%年平均


2.txt存在空行
思考：matlab删除空行，重新输出为txt,替换原文件
a=textread('SURF_CLI_CHN_TEM_MON_GRID_0.5-MEAN-201505.txt','%s','delimiter','\n');
a(cellfun(@isempty,a))=[]  %删除空行


3.nc数据过大无法读取
在读取nc数据时显示：
错误使用 netcdflib
请求的288*192*32*600（4.0GB）数组超过预设的最大数组大小。

解决方法：裁取数据
a1='F:\MATLAB\Bishe\CO2\co2_AERmon_CESM2_historical_r10i1p1f1_gn_195001-199912.nc';

lat=ncread(a1,'lat');
Lat=find(lat>=25 & lat<=40);%纬度范围
lon=ncread(a1,'lon');
Lon=find(lon>=73 & lon<=106);%经度范围
lat_num=length(Lat);
lon_num=length(Lon);

c2=ncread(a1,'co2',[Lon(1) Lat(1) 1 1],[lon_num lat_num 32 600]);
OK！！！

小莹子

%%各个格点多年相关&显著性检验

clc;
clear all;
t=xlsread('/home/syrus/桌面/ERA5/T-原.xlsx','W');%存储为一维数据格式，即格点总数*时间
b=xlsread('/home/syrus/桌面/ERA5/alb-原.xlsx','W');
%循环读取
t1=t;%(:,1:18);
b1=b;%(:,1:18);
cor=corrcoef(t1(1,:),b1(1,:));
k=size(t1,1);
for i=1:k
[r p]=corrcoef(t1(i,:),b1(i,:));
R(:,:,:,:,i)=r(:,:,:,:);
P(:,:,:,:,i)=p(:,:,:,:);
end
R1=R(1,2,1,1,:);
Rr=reshape(R1,[5525 1]);
P1=P(1,2,1,1,:);
Pp=reshape(P1,[5525 1]);
RP=[Rr Pp];

小莹子

%根据气压求海拔高度

%this programe is to figure out altitude from surface pressure
clc;
clear all;
D1='/data/ERA5-t-alb-snow-pre/2018.nc';
D2='/data/ERA5-t-alb-snow-pre/2019.nc';
d1=ncread(D1,'sp');
d2=ncread(D2,'sp');
d=cat(3,d1,d2);%合并
pressure=mean(d(13:181,1:81,1:24),3);%平均
lat=ncread(D1,'latitude');
lat=lat(1:81,:);
lon=ncread(D1,'longitude');
lon=lon(13:181,:);
% load dataTP\pressure.mat; %pressure, unit:Pa

p0=101325;  %sea level pressure
g=9.8;
r=0.0065; %degree/1000m
R=287;
T0=288; %temperature is 15 degree
m=R*r/g;
% grid=22321; %totao grid points are 161*81=13041
% mon=12;  %there are only one year 2018
k1=169;  %how many longitude grid. Here is from 65 to 105E
k2=81;  %how many latitude grid. Here is from 45 to 25N
k3=492; %how many month. here is (1997-2019) 41 years
%figure out altitude from surface pressure
alt=zeros(k1,k2);
for i=1:k1   %longitude from 70 to 105E
    for j=1:k2  %latitide from 45 to 25N
      % pz=mean(pressure(i,j,:)); %mean surface pressure
       pz=pressure(i,j);
       pp=pz/p0;
       alt(i,j)=(1-pp^m)*T0/r;
    end
end
s=size(alt,3);%T
for i=1:169
    for j=1:81
         xx(((i-1)*81+j),1)=lon(i);
         xx(((i-1)*81+j),2)=lat(j);
           for s=1:s
              xx(((i-1)*81+j),(2+s))=alt(i,j,s);%T
           end
    end
end
%%筛选
latlon=xx(:,1:2);
latlonT=importdata('/data/ERA5-t-alb-snow-pre/latlon/latlon-5525.txt');
out=xx(:,3);
lia=ismember(latlon,latlonT,'rows');
index=find(lia>0);
out=out(index,:);

小莹子

%%海拔依赖性

clc;
clear all;
A=xlsread('/home/syrus/桌面/ERA5-sd.xlsx','slope');
A=A(:,3:6);

S0=find(A(:,1)<2000);                  
C=A(S0,:);  
a0=mean(C);

S1=find(A(:,1)>=2000&A(:,1)<2200);                  
C=A(S1,:);  
a1=mean(C);

S2=find(A(:,1)>=2200&A(:,1)<2400);
C=A(S2,:);
a2=mean(C);  

S3=find(A(:,1)>=2400&A(:,1)<2600);
C=A(S3,:);
a3=mean(C);  

S4=find(A(:,1)>=2600&A(:,1)<2800);
C=A(S4,:);
a4=mean(C);  

S5=find(A(:,1)>=2800&A(:,1)<3000);
C=A(S5,:);
a5=mean(C);  

S6=find(A(:,1)>=3000&A(:,1)<3200);  
C=A(S6,:);
a6=mean(C);  

S7=find(A(:,1)>=3200&A(:,1)<3400);
C=A(S7,:);
a7=mean(C);  

S8=find(A(:,1)>=3400&A(:,1)<3600);
C=A(S8,:);
a8=mean(C);  

S9=find(A(:,1)>=3600&A(:,1)<3800);
C=A(S9,:);
a9=mean(C);  

S10=find(A(:,1)>=3800&A(:,1)<4000);
C=A(S10,:);
a10=mean(C);  

S11=find(A(:,1)>=4000&A(:,1)<4200);
C=A(S11,:);
a11=mean(C);  

S12=find(A(:,1)>=4200&A(:,1)<4400);
C=A(S12,:);
a12=mean(C);  

S13=find(A(:,1)>=4400&A(:,1)<4600);
C=A(S13,:);
a13=mean(C);  

S14=find(A(:,1)>=4600&A(:,1)<4800);         
C=A(S14,:);                     
a14=mean(C);               

S15=find(A(:,1)>=4800&A(:,1)<5000);
C=A(S15,:);                     
a15=mean(C);               

S16=find(A(:,1)>=5000&A(:,1)<5200);         
C=A(S16,:);                     
a16=mean(C);               

S17=find(A(:,1)>=5200&A(:,1)<5400);         
C=A(S17,:);                     
a17=mean(C);               

% S18=find(A(:,1)>=5400&A(:,1)<5600);         
% C=A(S18,:);                     
% a18=mean(C);               

S18=find(A(:,1)>=5400);        
C=A(S18,:);                     
a18=mean(C);               



W=[a0
    a1
   a2
   a3
   a4
   a5
   a6
   a7
   a8
   a9
   a10
   a11
   a12
   a13
   a14
   a15
   a16
   a17
   a18
   ];

小莹子

%%SVD
附录：SVD代码：
clear
clc

nxx=1099;
nyy=1099;
X=xlsread('F:\MATLAB\TP-T2m\SVD\lcc.xlsx','s');
Y=xlsread('F:\MATLAB\TP-T2m\SVD\dl.xlsx','s');
X=X';
Y=Y';

%*标准化数据
for i=1:nxx
    Stdrain(i)=std(X(:,i));
    X(:,i)=(X(:,i)-mean(X(:,i)))./Stdrain(i);
end

for i=1:nyy
    Stdptaux(i)=std(Y(:,i));
    Y(:,i)=(Y(:,i)-mean(Y(:,i)))./Stdptaux(i);
end

X=X';%将标准化后的二维sst场资料进行转置后赋给X
Y=Y';%将标准化后的二维同期uwind场资料进行转置后赋给Y
A=X*Y';%得到一个协方差矩阵A
[U,S,V]=svd(A);%对A进行SVD分解
U=-U;
V=-V;
Train=U'*X;%SVD分解后左向量和其中的一原始场X乘积后得到一个时间序列Train  第一行就是第一模态的时间序列
Tptaux=V'*Y;%SVD分解后右向量和其中的一原始场Y乘积后得到一个时间序列Tptaux

%****************** temporal correlation coefficient ***********\时间相关系数
%第一个是第一对空间分布型时间系数之间的相关系数
s=1099;
for i=1:s
    Sx=sum(Train(i,:).^2);
    Sy=sum(Tptaux(i,:).^2);
    r(i)=S(i,i)/sqrt(Sx*Sy);
end
%***************** varience contribution *******************

for i=1:s
    sigma(i)=S(i,i).^2;
end
%总的方差贡献率
for i=1:s
   var_con(i)=sigma(i)/sum(sigma);

end
sum(sigma)

%*******************
for i=1:s
    Srain(i)=sum(Train(i,:).^2);
end
%左场对右场贡献率
for i=1:s
    varrain_con(i)=Srain(i)/sum(Srain);%计算时间序列Train的各个元素的方差贡献大小
end
%*******************
for i=1:s
    Sptaux(i)=sum(Tptaux(i,:).^2);
end
%右场对左场贡献率
for i=1:s
    varptaux_con(i)=Sptaux(i)/sum(Sptaux);%计算时间序列Tptaux的各个元素的方差贡献大小
end
%************** Heterogeneous correlation ******************
%格点气温的相关系数
for i=1:nxx
     R=corrcoef(Tptaux(1,:),X(i,:));
     Hetero_rrain1(i)=R(1,2);
end
%格点低云的相关系数
for i=1:nyy
    R=corrcoef(Train(1,:),Y(i,:));
    Hetero_ptu1(i)=R(1,2);
end

si18

楼主要是有示例数据就更好了

lanyf

感谢分享，十分有用

zihan6688

这样只给出程序，让人很难看懂的