matlab读取GLDAS数据，并绘制全球格网图

全球陆地水同化模型（Global Land Data Assimilation System, GLDAS）由戈达德航天飞行中心和国家环境预测中心联合开发，包括四个模型：NOAH、MOSAIC、VIC和CLM。GLDAS模型输出的全球土壤水分、雪水当量、冠层水分和其他地表变量构成地表含水量的主要部分。本文以GLDAS_Noah格网分辨率为1°×1°的数据为例。

第一部分：数据下载和了解数据所包含的信息

1、登录earth data网站（GES DISC Search: Showing 1 - 25 of 42 datasets associated with gldas_noah (nasa.gov)）（不要科学上网），直接搜索gldas_noah，找到时间分辨率为1个月，空间分辨率为1°×1°的数据产品（GLDAS_NOAH10_M v2.1）。

2、可以发现数据在全球格网上，纬度的范围为90°--60°，说明产品不包含南极圈附近30°带宽范围内的数据。点击Subset/Get data下载数据。

3、选择2002年4月至2009年12月的数据为例，get data生成数据下载链接，可以先选中“download instructions”中的“download the list of links .txt”，该文本文件中包含所要下载所有数据的网址链接，后期可以利用IDM软件快速下载。（IDM抓取网页下载数据的教程可参照：利用IDM直接抓取网页上的数据，人为选择进行下载-CSDN博客）

4、下载完数据后，在matlab中了解数据所包含哪些信息，用到matlab自带“ncinfo”函数，查看变量info中的系列变量，对于各种变量的具体含义可以通过name简单了解，深入了解可以点开相应的“Attributes”变量或下载“README_GLDAS2.pdf”文档。

info = ncinfo('E:\Code\GLDAS\Noah\GLDAS_NOAH10_M.A200204.021.nc4');

第二部分：读取数据，并绘图

1、首先构造Contr.txt（控制文档），其结构如下：

（这件简单解释一下为什么下载2002年4月至2009年12月数据，控制文档中的数据有缺失，因为这是为了和GRACE卫星的数据缺失月份对应上，而人为做的删减操作。）

2、控制文件构造完后，即可按照如下代码进行数据读取和绘图：

clc;clearvars -except;
addpath('E:\Code\GLDAS\Noah');
info = ncinfo('E:\Code\GLDAS\Noah\GLDAS_NOAH10_M.A200204.021.nc4');
fid=fopen('Contr1.txt','r');
num_file= fscanf(fid,'%d',1);
lat= fscanf(fid,'%d',1);
lon= fscanf(fid,'%d',1);
dir_in=fscanf(fid,'%s',1);
file_name=cell(num_file,1);
for i = 1:num_file 
    file_name{i} = fscanf(fid,'%s',1);
end
gldas=zeros(num_file,lat,lon);
SWE=zeros(num_file,lat,lon);
SoilMoi0_10cm_inst=zeros(num_file,lat,lon);
SoilMoi100_200cm_inst=zeros(num_file,lat,lon);
SoilMoi10_40cm_inst=zeros(num_file,lat,lon);
SoilMoi40_100cm_inst=zeros(num_file,lat,lon);
CanopInt_inst=zeros(num_file,lat,lon);
time=zeros(num_file,1);

for ii=1:num_file
File = dir(fullfile(dir_in,'*.nc4'));	        % dir 函数读取.nc格式的文件名 'name' ，路径‘folder’ 等信息
    file = strcat(dir_in,File(ii).name);    % 拼接路径和文件名，并显示
    SWE(ii,:,:)=rot90(ncread(file,'SWE_inst'));
    SoilMoi0_10cm_inst(ii,:,:)=rot90(ncread(file,'SoilMoi0_10cm_inst'));
    SoilMoi10_40cm_inst(ii,:,:)=rot90(ncread(file,'SoilMoi10_40cm_inst'));
    SoilMoi40_100cm_inst(ii,:,:)=rot90(ncread(file,'SoilMoi40_100cm_inst'));
    SoilMoi100_200cm_inst(ii,:,:)=rot90(ncread(file,'SoilMoi100_200cm_inst'));
    CanopInt_inst(ii,:,:)=rot90(ncread(file,'CanopInt_inst'));
    time(ii,:,:)=ncread(file,'time');
end

gldas=zeros(num_file,180,360);
for j=1:num_file
    gldas(j,1:150,:)=(SWE(j,:,:)+SoilMoi0_10cm_inst(j,:,:)+SoilMoi100_200cm_inst(j,:,:)+...
        SoilMoi10_40cm_inst(j,:,:)+SoilMoi40_100cm_inst(j,:,:)+CanopInt_inst(j,:,:));%单位为m
end%积雪量、树冠总蓄水和四层土壤含水量累加

for j=1:num_file
    gldas_1(:,:)=gldas(j,:,:);
    gldas_1Q=gldas_1(:,1:180);gldas_1H=gldas_1(:,181:360);
    gldas_2=[gldas_1H,gldas_1Q]; gldas_2(find(isnan(gldas_2)==1))=0;%找到缺失值NaN赋值为0
    gldas(j,:,:)=gldas_2(:,:);
end
gldasm=mean(gldas(19:90,:,:));
for j=1:num_file
    gldas(j,:,:)=gldas(j,:,:)-gldasm(1,:,:);
end
    for i=1:num_file
        temp(:,:)=gldas(i,:,:);
        gldas_1(:,:,i)=temp(:,:);
    end

ii=1;
grid_1(:,:)=gldas_1(:,:,ii);
figure(2);
gmt_grid2map(grid_1*0.1,-30,30,0,0,'cm',['GLDAS Noah']);

PS：1、最后绘图用到冯伟老师编写的用于GRACE数据绘图的gmt_grid2map函数

2、如何批量复制文件名称并保存至同一个txt中，可以新建一个txt，并输入如下内容

@Echo off
dir /b>test.txt

再将txt后缀名改为.bat，然后将.bat文件置于需要复制文件名的路径下，运行.bat即可