读取GLDAS-1版本数据集Noah水文模型

clc;
clear;

% 批读取NC文件的准备工作
datadir = 'D:\MATLAB\toolbox\fengweiigg-GRACE_Matlab_Toolbox-6e91b05\GRACE_data\gldas_noah\'; %指定批量数据所在的文件夹
filelist = dir([datadir,'*.nc4']);       %列出所有满足指定类型的文件
%  a = filelist(1).name;                 %查看要读取的文件的编号
%  b = filelist(2).name;
k=length(filelist);

for i = 1:k  %依次读取并处理
    
    % 批量读取NC文件
    ncFilePath = ['D:\MATLAB\toolbox\fengweiigg-GRACE_Matlab_Toolbox-6e91b05\GRACE_data\gldas_noah\',filelist(i).name]; %设定NC路径
    num = filelist(i).name(1:26); %读取数据编号，以便于保存时以此编号储存tif
      
    % 读取变量值
    lon=ncread(ncFilePath,'lon');          %读取经度信息（范围、精度）360
    lat=ncread(ncFilePath,'lat');          %读取维度信息150
    time=ncread(ncFilePath,'time');        %读取时间序列
    su=ncread(ncFilePath,'SoilMoi100_200cm_inst');      %获取土壤湿度
    can=ncread(ncFilePath,'CanopInt_inst');     %获取树冠蓄水量
    sw=ncread(ncFilePath,'SWE_inst');           %获取积雪量
    
    sum_pre=sum(su+can+sw,3)/180;
   
%     %  展示数据内部结构等信息
%      pcolor(lat,lon,sum_pre);
%      shading flat;                        %移除网格线，否则图上一片黑什么都没有
%      [x,y]=meshgrid(lon,lat);             %根据经纬度信息产生格网，3600列（经度），1800列（纬度）
%      phandle=pcolor(x,y,sum_pre');     %显示一个矩阵，其中x,y,sum_pre的行列数必须一致
%      shading flat;
%      colorbar

    % 存为tif格式
      data=flipud(sum_pre');       %很重要，这是镜像反转，否则最后的图像的南北朝向是错的
%     R = georasterref('RasterSize', size(data),'Latlim', [double(min(lat)) double(max(lat))], 'Lonlim', [double(min(lon)) double(max(lon))]);
%     geotiffwrite(['D:\MATLAB\toolbox\fengweiigg-GRACE_Matlab_Toolbox-6e91b05\GRACE_results\gldas-ys\',num,'.tif'],data,R);
%     disp([num,'done'])
    % 存为mat格式
    save('D:\MATLAB\toolbox\fengweiigg-GRACE_Matlab_Toolbox-6e91b05\GRACE_results\result-ys\gldas-ys','data');
    
    
  
end
disp('finish!')