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:TarimBasinmetaERA5_original';
% 输出tif根路径
tif_Path='E:TarimBasinmetaERA5_originalnc2tif';
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
2.Stackoverflow:MATLAB: How to calculate total precipitation per day using hourly data ? (netcdf)
本文发布于:2024-02-04 06:31:32,感谢您对本站的认可!
本文链接:https://www.4u4v.net/it/170701110553145.html
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。
留言与评论(共有 0 条评论) |