本文主要是介绍Matlab : nc批处理输出指定时刻tif,Hourly转Daily并输出tif(以ERA5的nc文件为例),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
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); %读取数据年份编号,以便于保存时以此编号储存tifnc_vart2m = ncread(strcat(nc_Path,nc_name),'t2m'); %2m temperature,knc_vard2m = ncread(strcat(nc_Path,nc_name),'d2m'); %2m dewpoint temperature,knc_varu10 = ncread(strcat(nc_Path,nc_name),'u10'); %10m u-component of wind,m/snc_varv10 = ncread(strcat(nc_Path,nc_name),'v10'); %10m v-component of wind,m/s% 直接输出11点气象数据t2m,d2m,Uz,用于计算ET0_hrn=0;time=size(nc_vart2m);%层数量doy=time(3)/24;t=time(3)-12;Uz=zeros(time(1),time(2),doy);for i=12:24:tt11=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_dayfor j=0:(doy-1)t2m_day=nc_vart2m(:,:,(24*j+1):(24*(j+1))); %grab an entire dayTmax=max(t2m_day,[],3);%find Daily maximum temperatureTmin=min(t2m_day,[],3);%find Daily minimum temperatured2m_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
这篇关于Matlab : nc批处理输出指定时刻tif,Hourly转Daily并输出tif(以ERA5的nc文件为例)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!