MSWEP数据nc格式转tif格式

这里将数据从0.1°重采样到0.25°

%%
clc;
clear;

%%
fullpath = mfilename('fullpath');
[path, name] = fileparts(fullpath);

%%
ncpath = strcat(path, '\3hourly_nc\');
tifpath = strcat(path, '\tiffile_test\');
files = dir(strcat(ncpath, '*.nc'));
len = length(files);

for i = 1 : len
    filename = files(i).name;
    filename = filename(1 : size(filename, 2) - 3);
    ncfile = strcat(ncpath, filename, '.nc');
    ncinf = ncinfo(ncfile);
    
        
    if size(filename, 2) > 19
        filename = filename(1 : 19);
    end
    
    %
    stryear = filename(1 : 4);
    strdays = filename(5 : 7);
    strhour = filename(9 : 10);
    
    filedate = datetime(strcat(stryear, '-01-01')) + caldays(str2double(strdays) - 1);
    filedate = char(filedate);
    filedate = filedate([1 : 4, 6 : 7, 9 : 10]);
    filedate = strcat(filedate, strhour);
       
    OutputPath = strcat(tifpath, filedate, '.tif');
    
    if exist(OutputPath, 'file') == 2
        [num2str(i), '	', num2str(len), '	', filename, '	', filedate, '	continue']
        continue;
    end
    
    [num2str(i), '	', num2str(len), '	', filename, '	', filedate]
    
    preciSets = ncread(ncfile, 'precipitation');
    preciSets = rot90(preciSets);
    % preciSets = preciSets';
    lon = ncread(ncfile, 'lon');
    lat = ncread(ncfile, 'lat');
    
    [X, Y] = meshgrid(lon, lat);
    lat1 = 89.875 : -0.25 : -89.875;
    lon1 = -179.875: 0.25 : 179.875;
    [XI, YI] = meshgrid(lon1, lat1);
    preciSets1 = interp2(X, Y, preciSets, XI, YI, 'linear');
    
    Refference = georasterref('RasterSize', size(preciSets1), 'Lonlim', [-180 180], 'Latlim', [-90 90]);
    geotiffwrite(OutputPath, preciSets1, Refference);
end