Matlab : nc批处理输出指定时刻tif,Hourly转Daily并输出tif(以ERA5的nc文件为例)
该代码可以实现多个nc文件的循环批处理,上代码!
补充:
(1)1年的nc文件构成是:
(2)nc文件命名例如:ERA2020,所以截取年份时是4:7;
(3)使用时请根据自己的数据情况和需求更改代码适配。
%%
% @LastEditors : 喵小鱼~
% @Date : 2021-11-10 20:00
% @LastEditTime : 2021-11-20 16:50
% @Description : nc2tif,hourly2daily,year by year
% @Objective : ERA5小时数据转日数据,并输出11点数据
%%
clc;
clear all;
% 输入nc根路径
nc_Path='E:\TarimBasin\meta\ERA5_original\';
% 输出tif根路径
tif_Path='E:\TarimBasin\meta\ERA5_original\nc2tif\';
nc_list=dir(strcat(nc_Path,'*.nc'));%列出路径下所有nc文件
% ERA5小时数据转日数据,并输出11点数据
for yy=1:length(nc_list)
nc_name=nc_list(yy).name;
info=ncinfo(strcat(nc_Path,nc_name)); %读取nc文件信息
year =nc_list(yy).name(4:7); %读取数据年份编号,以便于保存时以此编号储存tif
nc_vart2m = ncread(strcat(nc_Path,nc_name),'t2m'); %2m temperature,k
nc_vard2m = ncread(strcat(nc_Path,nc_name),'d2m'); %2m dewpoint temperature,k
nc_varu10 = ncread(strcat(nc_Path,nc_name),'u10'); %10m u-component of wind,m/s
nc_varv10 = ncread(strcat(nc_Path,nc_name),'v10'); %10m v-component of wind,m/s
% 直接输出11点气象数据t2m,d2m,Uz,用于计算ET0_hr
n=0;
time=size(nc_vart2m);%层数量
doy=time(3)/24;
t=time(3)-12;
Uz=zeros(time(1),time(2),doy);
for i=12:24:t
t11=nc_vart2m(:,:,i);
dt11=nc_vard2m(:,:,i);
u10=nc_varu10(:,:,i);
v10=nc_varv10(:,:,i);
u11=sqrt(u10.^2+v10.^2);%先计算再采样与先采样再计算会改变精度吗?应该不会
n=n+1;
t11_file=fullfile(strcat(tif_Path,year,'\t11\'));
dt11_file=fullfile(strcat(tif_Path,year,'\dt11\'));
u11_file=fullfile(strcat(tif_Path,year,'\u11\'));
t11_Path=strcat(t11_file,'t11_doy',year,num2str(n,'%03d'),'.tif');
dt11_Path=strcat(dt11_file,'dt11_doy',year,num2str(n,'%03d'),'.tif');
u11_Path=strcat(u11_file,'u11_doy',year,num2str(n,'%03d'),'.tif');
t11 = rot90(t11);
dt11 = rot90(dt11);
u11 = rot90(u11);
Ref=georasterref('RasterSize',size(t11),'Latlim',[30 50],'Lonlim',[70 100]);
Ref.ColumnsStartFrom = 'south';
geotiffwrite(t11_Path,t11,Ref);
geotiffwrite(dt11_Path,dt11,Ref);
geotiffwrite(u11_Path,u11,Ref);
disp(i);
end
% 小时数据转日平均数据&日最值数据,Tmean,Tmax,Tmin,dTmean,meanUz用于计算ET0_hr & ET_day
for j=0:(doy-1)
t2m_day=nc_vart2m(:,:,(24*j+1):(24*(j+1))); %grab an entire day
Tmax=max(t2m_day,[],3);%find Daily maximum temperature
Tmin=min(t2m_day,[],3);%find Daily minimum temperature
d2m_day=nc_vard2m(:,:,(24*j+1):(24*(j+1)));
dTmean=sum(d2m_day,3)./24;
u10_day=nc_varu10(:,:,(24*j+1):(24*(j+1)));
v10_day=nc_varv10(:,:,(24*j+1):(24*(j+1)));
meanu10=sum(u10_day,3)./24;
meanv10=sum(v10_day,3)./24;
meanUz=sqrt(meanu10.^2+meanv10.^2);
Tmax_file=fullfile(strcat(tif_Path,year,'\maxTem\'));
Tmin_file=fullfile(strcat(tif_Path,year,'\minTem\'));
dTmean_file=fullfile(strcat(tif_Path,year,'\dTmean\'));
meanUz_file=fullfile(strcat(tif_Path,year,'\meanUz\'));
Tmax_Path=strcat(Tmax_file,'Tmax_doy',year,num2str(j+1,'%03d'),'.tif');
Tmin_Path=strcat(Tmin_file,'Tmin_doy',year,num2str(j+1,'%03d'),'.tif');
dTmean_Path=strcat(dTmean_file,'dTmean_doy',year,num2str(j+1,'%03d'),'.tif');
meanUz_Path=strcat(meanUz_file,'meanUz_doy',year,num2str(j+1,'%03d'),'.tif');
Tmax = rot90(Tmax);
Tmin = rot90(Tmin);
dTmean = rot90(dTmean);
meanUz = rot90(meanUz);
Ref=georasterref('RasterSize',size(Tmax),'Latlim',[30 50],'Lonlim',[70 100]);
Ref.ColumnsStartFrom = 'south';
geotiffwrite(Tmax_Path,Tmax,Ref);
geotiffwrite(Tmin_Path,Tmin,Ref);
geotiffwrite(dTmean_Path,dTmean,Ref);
geotiffwrite(meanUz_Path,meanUz,Ref);
disp(j+1);
end
end
参考:
1.CSDN:matlab将nc数据转换为tif
https://blog.csdn.net/qq_34149805/article/details/73024600
2.Stackoverflow:MATLAB: How to calculate total precipitation per day using hourly data ? (netcdf)
https://stackoverflow.com/questions/58658489/matlab-how-to-calculate-total-precipitation-per-day-using-hourly-data-netcdf