一个脚本文件读取CSR、JPL和GSFC的mascon产品，并绘制流域时序图

1、列出获取mascon数据的网站

CSR：https://www2.csr.utexas.edu/grace/

JPL：https://podaac.jpl.nasa.gov/dataset/TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.1_V3 (可能需要科学上网)

https://podaac.jpl.nasa.gov/grace/?tab=mission-objectives&sections=about%2Bdata

NASA-GSFC：https://earth.gsfc.nasa.gov/geo/data/grace-mascons

或者参考博文：整理分享GRACE领域可能常用的网站-CSDN博客

2、数据名称预处理

由于下载的JPL mascon数据为“.nc4”文件，于是将其修改为“.nc”文件；并将文件名开头修改为“JPL_”，以与其它机构的mascon文件名对应起来

3、读取数据

clear
clc
% info_CSR = ncinfo('E:\Data\mascon\data\CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.nc');
% info_JPL = ncinfo('E:\Data\mascon\data\GRCTellus.JPL.200204_202207.GLO.RL06M.MSCNv02CRI.nc');
% info_GSFC = ncinfo('E:\Data\mascon\data\gsfc.glb_.200204_202112_rl06v2.0_obp-ice6gd_halfdegree.nc');

Path = 'E:\Data\mascon\data\';			% 设置路径，记得加上最后的反斜杠
File = dir(fullfile(Path,'*.nc'));	        % dir 函数读取.nc格式的文件名 'name' 
Len = length(File);				% 读取文件数量

for i=1:Len
    full_path = strcat(Path,File(i).name);    % 拼接路径和文件名，并显示
    Var_name='time';
    eval(['time_' File(i,1).name(1,1:4) '=ncread(full_path,Var_name);']);%%依次读取变量
    Var_name='time_bounds';
    eval(['time_bounds_' File(i,1).name(1,1:4) '=ncread(full_path,Var_name);']);
    Var_name='lwe_thickness';
    eval(['lwe_thickness_' File(i,1).name(1,1:4) '=ncread(full_path,Var_name);']);
    num(i,1)=size(eval(['time_' File(i,1).name(1,1:4)]),1);
end

for i=1:Len
    for ii=1:num(i,1)
        temp1(:,:)=eval(['lwe_thickness_' File(i,1).name(1,1:4) '(:,:,ii);']);
        temp=temp1';
        eval(['mas_' File(i,1).name(1,1:4) '(:,:,ii)=flipud(temp);']);
    end
    clear temp temp1;
end

grid_CSR(:,:)=mas_CSR_(:,:,1);
gmt_grid2map(grid_CSR,-30,30,0,0,'EWH (cm)',['CSR 2002年4月'],20)
grid_JPL(:,:)=mas_JPL_(:,:,1);
gmt_grid2map(grid_JPL,-30,30,0,0,'EWH (cm)',['JPL 2002年4月'],20)
grid_GSFC(:,:)=mas_gsfc(:,:,1);
gmt_grid2map(grid_GSFC,-30,30,0,0,'EWH (cm)',['GSFC 2002年4月'],20)

以2002年4月的全球EWH格网点为例，

4、分别求取亚马逊、刚果河、长江和密西西比河流域的水储量变化，并绘制2002年4月至2010年12月的时序图

for i=1:Len
mas=eval(['mas_' File(i,1).name(1,1:4);]);
iindex=1;
for jj=1:4
    switch jj
        case 1
            mask = 'E:\Data\Basin\Amazon.vec';rows=211;
        case 2
            mask = 'E:\Data\Basin\Congo.vec';rows=181    ;
        case 3
            mask = 'E:\Data\Basin\Yangtze.vec';rows=12191;
        case 4
            mask = 'E:\Data\Basin\Mississippi.vec';rows=203;
    end
    eval(['TWSA_' File(i,1).name(1,1:4) '(:,jj)=gmt_grid2serie(mas,mask,[],rows);']);

end
end

FileNameTime02={'2002-04','2002-05','2002-08','2002-09','2002-10','2002-11','2002-12'};
FileNameTime03={'2003-01','2003-02','2003-03','2003-04','2003-05','2003-07','2003-08','2003-09','2003-10','2003-11','2003-12'};
FileNameTime04={'2004-01','2004-02','2004-03','2004-04','2004-05','2004-06','2004-07','2004-08','2004-09','2004-10','2004-11','2004-12'};
FileNameTime05={'2005-01','2005-02','2005-03','2005-04','2005-05','2005-06','2005-07','2005-08','2005-09','2005-10','2005-11','2005-12'};
FileNameTime06={'2006-01','2006-02','2006-03','2006-04','2006-05','2006-06','2006-07','2006-08','2006-09','2006-10','2006-11','2006-12'};
FileNameTime07={'2007-01','2007-02','2007-03','2007-04','2007-05','2007-06','2007-07','2007-08','2007-09','2007-10','2007-11','2007-12'};
FileNameTime08={'2008-01','2008-02','2008-03','2008-04','2008-05','2008-06','2008-07','2008-08','2008-09','2008-10','2008-11','2008-12'};
FileNameTime09={'2009-01','2009-02','2009-03','2009-04','2009-05','2009-06','2009-07','2009-08','2009-09','2009-10','2009-11','2009-12'};
FileNameTime10={'2010-01','2010-02','2010-03','2010-04','2010-05','2010-06','2010-07','2010-08','2010-09','2010-10','2010-11','2010-12'};
FileNameTime=[FileNameTime02,FileNameTime03,FileNameTime04,FileNameTime05,FileNameTime06,FileNameTime07,FileNameTime08,FileNameTime09,FileNameTime10];

num_file=size(FileNameTime,2);
FileNameTimeChar=char(FileNameTime);
int_year=zeros(num_file,1);int_month=zeros(num_file,1);%用于存储年和月
for i=1:num_file
    int_year(i)=str2double(FileNameTimeChar(i,1:4));
    int_month(i)=str2double(FileNameTimeChar(i,6:7));
end
time=int_year+(int_month-0.5)/12;

close all;
name={'Amazon','Congo','Yangtze','Mississippi'};
for jj=1:4
% jj=1;
figure(1)
subplot(2,2,jj)
plot(time,TWSA_CSR_(1:size(time,1),jj),'color',[0 1 0],'Linewidth',1.5);hold on;
plot(time,TWSA_JPL_(1:size(time,1),jj),'color',[0 0 1],'Linewidth',1.5);
plot(time,TWSA_gsfc(1:size(time,1),jj),'color',[1 0 0],'Linewidth',1.5);

% yline(0,'LineStyle','-','Color','red','LabelVerticalAlignment','middle','LabelHorizontalAlignment','center','LineWidth',1.5);
xlim([2002 2011]);
xticks([2002 2003 2004 2005 2006 2007 2008 2009 2010 2011]);
xlabel('时间 (年份)','Fontname','Time New Roman',"FontSize",15);   
l1=legend('CSR-mas','JPL-mas','GSFC-mas','location','northeast');
set(l1,'box','off',"FontSize",10,'Fontname','Time New Roman');
ylabel('TWSA (cm)','Fontname','Time New Roman',"FontSize",15);
box on;grid on;
title(name{1,jj},'Fontname','Time New Roman',"FontSize",15);
end

 注：mascon产品与GRACE Level-2月度球谐的缺值情况一致，如下图（自制）：

若有任何问题，请批评指正，欢迎点赞收藏+评论