ERA5 处理代码

小莹子

一、均值

clc;         
clear all;  

data='/data/ERA5-t-alb-snow-pre/';   
file=dir([data,'*.nc']);                              
k=length(file);                              
%循环读取文件
for s=1:k-1 %1979-2019,2020年去掉
  f=[data,file(s).name];
  ncid=netcdf.open(f,'NC_NOWRITE');     
  time=ncread(f,'time');
  t=ncread(f,'t2m')-273.15;% t2m; scale_factoe=0.0010483 add_offset=278.1858
  t=t(13:181,1:81,:);%65 to 105E  45 to 25N
  lat=ncread(f,'latitude');
  lat=lat(1:81,:);
  lon=ncread(f,'longitude');
  lon=lon(13:181,:);
  netcdf.close(ncid);        
  tem(:,:,:,s)=t(:,:,:);
end
Tt=reshape(tem,[169 81 492]);%41*12=492
a=ncread('/data/ERA5-t-alb-snow-pre/2020.nc','t2m')-273.15;
a=a(13:181,1:81,:);
T=cat(3,Tt,a);  %合并两个三维数据
Ttb=reshape(T(:,:,1:492),[169 81 12 41]);
T0=reshape(mean(Ttb(:,:,1:12,:),3),[169 81 41]);%年平均
T2=reshape(mean(Ttb(:,:,6:8,:),3),[169 81 41]);%夏季平均
%冬季 1-12----3月-2月
Ttt=reshape(T(:,:,3:494),[169 81 12 41]);
T4=reshape(mean(Ttt(:,:,10:12,:),3),[169 81 41]);%冬季平均 %3月-2月
%%气候态平均
s0=mean(T0,2);
s2=mean(T2,2);
s4=mean(T4,2);
smean=[s0 s2 s4];
%%%二维转一维
w=size(T0,3);%mean
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 w=1:w
            xx(((i-1)*81+j),(2+w))=T4(i,j,w);%mean T
           end
    end
end
%%筛选
latlon=xx(:,1:2);
latlonT=importdata('/data/ERA5-t-alb-snow-pre/latlon/latlon-5525.txt');
out=xx(:,3:43);%% xx(:,1:2)=[];%去掉前两列
lia=ismember(latlon,latlonT,'rows');
index=find(lia>0);
out=out(index,:);
%xlswrite('/home/syrus/桌面/ERA5-TP.xlsx',out,'sheet1');

mm=(mean(out,1))'; %多年均值
%slope结果t检验
for s=1:k
     y=B(s,:);
    [H,P,CI, STATS]=ttest2(x,y,0.05);
   S(:,s)=P;
   Ss=S';
end
% %matlab画图
[Lon Lat]=meshgrid(lon,lat);
% [lon lat]=meshgrid([70:0.5:105],[20:0.5:45]);
  TP=shaperead('/data/tibetan/TPBoundary_2500/TPBoundary_2500.shp');
TPX=[TP.X];
TPY=[TP.Y];
isin=inpolygon(Lon,Lat,TP.X,TP.Y);
isin2=isin';
%  for i=1:201%范围之外为空
%      for j=1:121
%          if(isin2(i,j)==0)
%              smean(i,j)=NaN;
%          end
%      end
%  end

[c,h]=contourf(lon,lat,slope');
set(h,'color','none');
shading flat
colormap(jet)%jet  hot hsv
colorbar
hold on

%  contour(lon,lat,smean',[0.0 0.0],'w','LineWidth',3); % k黑色 b蓝色 w白色
%  hold on

plot(TPX,TPY,'-k','linewidth',2);
hold on
set(gca,'LineWidth',2,'FontSize',15','Ylim',[25 45],'Xlim',[60 105]...
    ,'XTick',[60:5:105],'XTicklabel',{'60E','65E','70E','75E','80E','85E','90E','95E','100E','105E'}...
    ,'YTick',[25:5:45],'YTicklabel',{'25N','30N','35N','40N','45N'});
set(gca,'XGrid','off');
set(gca,'YGrid','off');
title(['1980-2019 ERA5 Y-slope '],'FontSize',16);
saveas(gcf,['/home/syrus/桌面/出图/w ERA5 .png']);

二、斜率
clc;         
clear all;  

data='/data/ERA5-t-alb-snow-pre/';   
file=dir([data,'*.nc']);                              
k=length(file);                              
%循环读取文件
for s=1:k-1 %1979-2019,2020年去掉
  f=[data,file(s).name];
  ncid=netcdf.open(f,'NC_NOWRITE');     
  time=ncread(f,'time');
  t=ncread(f,'t2m')-273.15;% t2m; scale_factoe=0.0010483 add_offset=278.1858
  t=t(13:181,1:81,:);%65 to 105E  45 to 25N
  lat=ncread(f,'latitude');
  lat=lat(1:81,:);
  lon=ncread(f,'longitude');
  lon=lon(13:181,:);
  netcdf.close(ncid);        
  tem(:,:,:,s)=t(:,:,:);
end
Tt=reshape(tem,[169 81 492]);%41*12=492
a=ncread('/data/ERA5-t-alb-snow-pre/2020.nc','t2m')-273.15;
a=a(13:181,1:81,:);
T=cat(3,Tt,a);  %合并两个三维数据
Ttb=reshape(T(:,:,1:492),[169 81 12 41]);
T0=reshape(mean(Ttb(:,:,1:12,:),3),[169 81 41]);%年平均
T2=reshape(mean(Ttb(:,:,6:8,:),3),[169 81 41]);%夏季平均
%冬季 1-12----3月-2月
Ttt=reshape(T(:,:,3:494),[169 81 12 41]);
T4=reshape(mean(Ttt(:,:,10:12,:),3),[169 81 41]);%冬季平均 %3月-2月
%%figure out the linear trend最小二乘法直线
x=1979:1:2019;
x=x';
for i=1:169
    for j=1:81
        y1=reshape(T4(i,j,:),41,1);%%%%
        [a b]=linearfit(x,y1);
        slope(i,j)=b*10;
    end
end
%斜率出图 转后文
%%%二维转一维
w=size(slope,3);
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 w=1:w
               xx(((i-1)*81+j),(2+w))=slope(i,j,w);
           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,:);
%xlswrite('/home/syrus/桌面/ERA5-TP.xlsx',out,'sheet1');

% %matlab画图
[Lon Lat]=meshgrid(lon,lat);
% [lon lat]=meshgrid([70:0.5:105],[20:0.5:45]);
  TP=shaperead('/data/tibetan/TPBoundary_2500/TPBoundary_2500.shp');
TPX=[TP.X];
TPY=[TP.Y];
isin=inpolygon(Lon,Lat,TP.X,TP.Y);
isin2=isin';
%  for i=1:201
%      for j=1:121
%          if(isin2(i,j)==0)
%              smean(i,j)=NaN;
%          end
%      end
%  end

[c,h]=contourf(lon,lat,slope');
set(h,'color','none');
shading flat
colormap(jet)%jet  hot hsv
colorbar
hold on

%  contour(lon,lat,smean',[0.0 0.0],'w','LineWidth',3); % k黑色 b蓝色 w白色
%  hold on

plot(TPX,TPY,'-k','linewidth',2);
hold on
set(gca,'LineWidth',2,'FontSize',15','Ylim',[25 45],'Xlim',[60 105]...
    ,'XTick',[60:5:105],'XTicklabel',{'60E','65E','70E','75E','80E','85E','90E','95E','100E','105E'}...
    ,'YTick',[25:5:45],'YTicklabel',{'25N','30N','35N','40N','45N'});
set(gca,'XGrid','off');
set(gca,'YGrid','off');
title(['1980-2019 ERA5 W-slope '],'FontSize',16);
saveas(gcf,['/home/syrus/桌面/出图/w ERA5 .png']);
% %斜率
% A=xx(:,3:43);
% x=1979:1:2019;
% k1=size(A,1);
% for j =1:k1
%     y=A(j,:);
%     [a m b]=regression(x,y);
%     m=m*10;
%     M(:,j)=m;
% end
%     N=M';

二.2 linearfit
%%
function [a,b]=linearfit(x,y)
xy=x.*y; x2=x.^2;
x_mean=mean(x);
y_mean=mean(y);
xy_mean=mean(xy);
x2_mean=mean(x2);
b=(xy_mean-x_mean*y_mean)/(x2_mean-x_mean^2);
a=y_mean-b*x_mean;
return

vrcsds

这个应该就是地面或海平面吧，猜测对应EC计算时使用的地形数据

小莹子

vrcsds 发表于 2020-6-8 16:19
这个应该就是地面或海平面吧，猜测对应EC计算时使用的地形数据

说是用地表气压计算个海拔高度

小莹子

  t=ncread(f,'t2m')-273.15;%%处理结果比真实值高和
  t=ncread(f,'t2m')*0.0010483+278.1858-273.15;% t2m; scale_factoe=0.0010483 add_offset=278.1858  %%%处理结果不合实际，说用ncread提取之后就不用scale