Matlab : nc批处理输出指定时刻tif,Hourly转Daily并输出tif(以ERA5的nc文件为例)

1 篇文章 0 订阅

Matlab : nc批处理输出指定时刻tif,Hourly转Daily并输出tif(以ERA5的nc文件为例)
该代码可以实现多个nc文件的循环批处理,上代码!
补充:
(1)1年的nc文件构成是:
ERA5的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

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值