Matlab读取Swarm球谐系数，并绘制EWH全球格网图（存在疑问）

ICGEM官网下载 COST-G发布的40×40的球谐系数

close all;
clearvars -except;
% addpath('E:\Code\Tool\Function\GRACE_functions');
dir_degree_1 = 'E:\Code\GRACE_data\Degree_1\deg1_coef.txt';
dir_c20      = 'E:\Code\GRACE_data\Degree_2\C20_RL06.txt';
myDir_Swarm        = 'E:\Data\Swarm\EGF\'; % COST-G 
myFiles_gfc_Swarm  = dir(fullfile(myDir_Swarm,'*.gfc'));
num_file    = length(myFiles_gfc_Swarm);

lmax=40;name='Swarm';
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','SSM');
 pathname=strcat(myDir_Swarm,myFiles_gfc_Swarm(ii,1).name);
 [cs(ii,:,:),cs_sgi(ii,:,:),int_year(ii),int_month(ii),time(ii)] =gmt_readgfc(pathname);
 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(:,:,:));
for i=1:num_file
    cs_res(i,:,:)  = cs_replace(i,:,:)-cs_mean(1,:,:);%扣除平均值
end
CS_(:,:)=cs_res(1,:,:);

for i=1:num_file
cs_tmp(:,:) = cs_res(i,:,:);
cs_tmp(:,:)=gmt_gc2mc(cs_tmp);
cs_mss(i,:,:) = gmt_destriping(cs_tmp,'NONE');
end
eval(['grid_' name '=gmt_cs2grid(cs_mss,750,1);']);

ii=1;
eval(['grid_1(:,:)=grid_' name '(:,:,ii);']);
gmt_grid2map(grid_1.*100,-30,30,1,0,'cm',[name   ' ' num2str(int_year(ii),'%02d') '-' num2str(int_month(ii),'%02d')],20)

未滤波处理结果：

 750km高斯滤波结果：

 从结果可以看出40阶的Swarm球谐系数反演的噪声太大，即使经过强约束的高斯滤波后效果仍然不理想，不太清楚是Swarm数据本身所致还是代码实现有误？还请有懂的友友帮忙解答下，不胜感激！

点赞支持，一键分享。欢迎评论交流！