ITSG、COST-G、Tongji和WHU Level-2数据产品读取绘图（Matlab）

数据介绍：

ICGEM International Center for Global Gravity Field Models (gfz-potsdam.de)

ITSG 2018：Institute of Geodesy at Graz University of Technolog（格拉茨理工大学大地测量研究所） 2018版本，最高60阶球谐系数

COST-G：时变重力场解决方案组合服务解RL01版本，最高96阶球谐系数

Time-Variable Gravity Field Modeling and Simulation from Present and Future Gravity Satellite Missions | ISSI Team led by Wei Feng (CN) (issibern.ch)

Tongji：同济大学2022版本，最高96阶球谐系数

WHU：武汉大学基于卫星间位势差的约束GRACE月重力场模型

代码中将所有球谐系数均人为截断至60阶，只对球谐系数做最基本的替换低阶项和扣除平均重力场处理，其余包括GIA改正、滤波处理和泄漏改正等均未涉及；

clearvars -except
addpath('H:\CSDN\Control\');
fid=fopen('Control_Tongji.txt','r');name='Tongji';
% fid=fopen('Control_WHU.txt','r');name='WHU';
% fid=fopen('Control_ITSG.txt','r');name='ITSG';
% fid=fopen('Control_COST_G.txt','r');name='COST-G';
num_file= fscanf(fid,'%d',1);
dir_c20 = fscanf(fid,'%s',1);dir_degree_1= fscanf(fid,'%s',1);
dir_in = fscanf(fid,'%s',1);dir_out= fscanf(fid,'%s',1);

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=[FileNameTime04,FileNameTime05,FileNameTime06,FileNameTime07,FileNameTime08,FileNameTime09,FileNameTime10];
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;

file_name=cell(num_file,1);
for i = 1:num_file
    file_name{i} = fscanf(fid,'%s',1);
end
fclose(fid);

lmax=60;
cs= zeros(num_file,lmax+1,lmax+1);
cs_sgi= zeros(num_file,lmax+1,lmax+1);
cs_res= zeros(num_file,lmax+1,lmax+1);
cs_mss= zeros(num_file,lmax+1,lmax+1);

hwait=waitbar(0,'Waiting>>>>>>>>');  
for ii=1:num_file
str=['Processing...',num2str(ii),'/',num2str(num_file),'    '];
hwait=waitbar(ii/num_file,hwait,str,'Name','Reading');
pathname=strcat(dir_in,file_name{ii,1});
head_index=0;
fid = fopen(pathname,'r');
tline = fgetl(fid);
while size(tline,2)<11 || ~strcmp(tline(1:11),'end_of_head')
      head_index = head_index+1;    tline = fgetl(fid);
end
fclose(fid);

[~, l, m, Clm, Slm, Clm_sigma, Slm_sigma] = textread(pathname,'%s %u %u %f %f %f %f','headerlines',head_index+1);

for i = 1:1891  %%截断至60阶
% for i=1:length(l)    
    sc_tmp(l(i)+1,lmax+1-m(i)) = Slm(i);
    sc_tmp(l(i)+1,lmax+1+m(i)) = Clm(i);
end
cs(ii,:,:)=gmt_sc2cs(sc_tmp);
end

cs_replace=cs;
[cs_replace] = gmt_replace_degree_1(dir_degree_1,cs_replace,int_year,int_month,num_file);
[cs_replace] = gmt_replace_C20(dir_c20,cs_replace,int_year,int_month,num_file);

cs_mean = mean(cs_replace(1:72,:,:)); %扣除2004年1月-2009年12月的均值

for i=1:num_file
    cs_res(i,:,:)  = cs_replace(i,:,:)-cs_mean(1,:,:);
end

for i=1:num_file
cs_tmp(:,:) = cs_res(i,:,:);
cs_mss(i,:,:)=gmt_gc2mc(cs_tmp);
end

grid=gmt_cs2grid(cs_mss,0,1);

for ii=1:num_file
    if int_year(ii)==2004
        if int_month(ii)==2
            grid_1(:,:)=grid(:,:,ii);
            figure(1);
            gmt_grid2map(grid_1*100,-30,30,0,0,'EWH (cm)',[name '  ' num2str(int_year(ii),'%02d') '-' num2str(int_month(ii),'%02d') ],20);
            set(gcf,'unit','centimeters','position',[1,2,30,15]);
        end
    end
end