GRACE RL06版本的数据预处理

function GRACE_RL06_preprocessing(controlfile_path)
% Read the Control File
fid=fopen(controlfile_path,'r');
num_file        = fscanf(fid,'%d',1);
radius_filter   = fscanf(fid,'%d',1);
destrip_method  = fscanf(fid,'%s',1);
% Check destriping method option
if ~strcmp(destrip_method,'NONE') ...
        && ~strcmp(destrip_method,'SWENSON') ...
        && ~strcmp(destrip_method,'CHAMBERS2007') ...
        && ~strcmp(destrip_method,'CHAMBERS2012') ...
        && ~strcmp(destrip_method,'CHENP3M6') ...
        && ~strcmp(destrip_method,'CHENP4M6') ...
        && ~strcmp(destrip_method,'DUAN')
    warndlg('Destrping method in the control file is wrong! (Options: NONE,SWENSON,CHAMBERS2007,CHAMBERS2012,CHENP3M4,CHENP4M6 and DUAN)','Warning');
    return;
end

option_gia      = fscanf(fid,'%s',1);
% Check GIA option
if ~strcmp(option_gia,'GIA_notRemoved') && ~strcmp(option_gia,'GIA_Removed_Geru')
    warndlg('GIA option in the control file is wrong! (Options: GIA_notRemoved and GIA_Removed_Geru)','Warning');
    return;
end

InputFileFormat = fscanf(fid,'%s',1);
% Check Input file format option
if ~strcmp(InputFileFormat,'GRACE') && ~strcmp(InputFileFormat,'ICGEM')
    warndlg('Format of input files in the control file is wrong! (Options: GRACE and ICGEM)','Warning');
    return;
end

OutputFileFormat_1= fscanf(fid,'%s',1);
% Check Output file format option
if strcmp(OutputFileFormat_1,'SH_MAT')
    OutputFileFormat_2= fscanf(fid,'%f',1);% max degree
    OutputFileFormat_3= fscanf(fid,'%s',1); % output name
elseif strcmp(OutputFileFormat_1,'GRID_MAT')
    OutputFileFormat_2= fscanf(fid,'%f',1); % max degree
    OutputFileFormat_3= fscanf(fid,'%s',1); % output name
    OutputFileFormat_4= fscanf(fid,'%f',1); % spatial resolution
    if OutputFileFormat_4~=1 && OutputFileFormat_4~=0.25
        warndlg('Resolution of output grids in the control file is wrong! (Options: 1 and 0.25)','Warning');
        return;
    end
else
    warndlg('Output format in the control file is wrong! (Options: SH_MAT and GRID_MAT)','Warning');
    return;
end

dir_c20         = fscanf(fid,'%s',1);
dir_c21_s21     = fscanf(fid,'%s',1);
dir_c22_s22     = fscanf(fid,'%s',1);
dir_degree_1    = fscanf(fid,'%s',1);
dir_in          = fscanf(fid,'%s',1);
dir_out         = fscanf(fid,'%s',1);
% Check files and directories
if ~exist(dir_c20,'file') && ~strcmp(dir_c20,'NAN')
    warndlg('Degree C20 file does not exist!','Warning');
    return;
end
if ~exist(dir_c21_s21,'file') && ~strcmp(dir_c21_s21,'NAN')
    warndlg('Degree C21 & S21 file does not exist!','Warning');
    return;
end
if ~exist(dir_c22_s22,'file') && ~strcmp(dir_c22_s22,'NAN')
    warndlg('Degree C22 & S22 file does not exist!','Warning');
    return;
end
if ~exist(dir_degree_1,'file') && ~strcmp(dir_degree_1,'NAN')
    warndlg('Degree 1 file does not exist!','Warning');
    return;
end
if ~exist(dir_in,'dir')
    warndlg('Directory of Input GRACE files does not exist!','Warning');
    return;
end
if ~exist(dir_out,'dir')
    warndlg('Directory of Output files does not exist!','Warning');
    return;
end

file_name=cell(num_file,1);% Name list of input GRACE files
for i = 1:num_file % read the name list of files in control file
    file_name{i} = fscanf(fid,'%s',1);
    if ~exist([dir_in,file_name{i}],'file')
        warndlg(strcat('Input GRACE files',[dir_in,file_name{i}],'does not exist!'),'Warning');
        return;
    end
end
fclose(fid);

%
% GSM-2_2003060-2003090_0030_UTCSR_0096_0005.gfc
% GSM-2_2004001-2004013_0013_UTCSR_0096_0005
% GSM-2_2002102-2002120_0018_UTCSR_0060_0004
%
% GSM-2_2002122-2002137_0013_EIGEN_G---_005a.gfc
% GSM-2_2003032-2003059_0028_EIGEN_G---_0005
% GSM-2_2003001-2003031_0029_JPLEM_0000_0005
% GSM-2_2002305-2002334_0028_JPLEM_0000_0005.gfc
%
%  ------------------------------------------------------
%   Initialize
%  ------------------------------------------------------

lmax=OutputFileFormat_2; % maximum degree of output files

cs              = zeros(num_file,lmax+1,lmax+1);
cs_sgi          = zeros(num_file,lmax+1,lmax+1);
int_year        = zeros(num_file,1);
int_month       = zeros(num_file,1);
time            = zeros(num_file,1);


cs_res          = zeros(num_file,lmax+1,lmax+1);
cs_destrip      = zeros(num_file,lmax+1,lmax+1);
cs_mss          = zeros(num_file,lmax+1,lmax+1);
cs_fltr         = zeros(num_file,lmax+1,lmax+1);

if strcmp(OutputFileFormat_1,'GRID_MAT')
    if strcmp(OutputFileFormat_4,'1')
        grid_data_grace = zeros(180,360,num_file);
    elseif strcmp(OutputFileFormat_4,'0.25')
        grid_data_grace = zeros(721,1440,num_file);
    end
end
%
% cs_tmp          = zeros(lmax+1,lmax+1);
% sc_tmp          = zeros(lmax+1,2*lmax+1);
% grid_tmp        = zeros(180,360);

% region_plot = zeros(num_file,1);
% region_plot_csr_error = zeros(num_file,1);


hwait=waitbar(0,'Waiting>>>>>>>>');
% ----------------------------------------
% Read the GRACE Files
% ----------------------------------------
    for ii=1:num_file
        waitbar(0.1*ii/num_file,hwait,'Read the GRACE Files');
        [cs(ii,:,:),cs_sgi(ii,:,:),int_year(ii),int_month(ii),time(ii)] = gmt_readgsm_RL06(dir_in,file_name{ii},lmax,'ICGEM');
    end

% ----------------------------------------
% Read and Replace C20 and Degree_1
% ----------------------------------------
waitbar(0.2,hwait,'Replace Degree 2 and Degree 1');
pause(0.5);

cs_replace=cs;
if ~strcmp(dir_degree_1,'NAN')
    [cs_replace,tag] = gmt_replace_degree_1(dir_degree_1,cs_replace,int_year,int_month,num_file);
    if tag==0
        warndlg('Format of degree 1 file is wrong!','Warning');
        close(hwait);
        return;
    end
end
if ~strcmp(dir_c20,'NAN')
    [cs_replace,tag]  = gmt_replace_C20(dir_c20,cs_replace,int_year,int_month,num_file);
    if tag==0
        warndlg('Format of C20 file is wrong!','Warning');
        close(hwait);
        return;
    end
end
if ~strcmp(dir_c21_s21,'NAN')
    [cs_replace,tag] = gmt_replace_C21_S21_C22_S22(dir_c21_s21,cs_replace,int_year,int_month,num_file);
    if tag==0
        warndlg('Format of C21 & S21 file is wrong!','Warning');
        close(hwait);
        return;
    end
end
if ~strcmp(dir_c22_s22,'NAN')
    [cs_replace,tag] = gmt_replace_C21_S21_C22_S22(dir_c22_s22,cs_replace,int_year,int_month,num_file);
    if tag==0
        warndlg('Format of C22 & S22 file is wrong!','Warning');
        close(hwait);
        return;
    end
end


% ----------------------------------------
% Remove average
% ----------------------------------------
waitbar(0.3,hwait,'Remove the average gravity field');
pause(0.5);

if num_file>1
    cs_mean = mean(cs_replace);
    for i=1:num_file
        cs_res(i,:,:)  = cs_replace(i,:,:)-cs_mean(1,:,:);
    end
else % only one input files
    cs_res=cs_replace;
end

% ----------------------------------------
% Destriping
% ----------------------------------------
for i=1:num_file
    %     waitbar(0.3+0.4*i/num_file,hwait,'Do destriping');
    waitbar(0.3+0.5*i/num_file,hwait,'Convert geoid to mass & Do destriping');
    cs_tmp(:,:) = cs_res(i,:,:);
    
    % ----------------------------------------
    % Geoid to Mass
    % ----------------------------------------
    cs_tmp(:,:)=gmt_gc2mc(cs_tmp);

%     %     % TEST
%     if i==17
%         ttt(:,:)=gmt_gc2mc(cs_tmp);
%         cs_destrip(i,:,:) = gmt_destriping_test(ttt,destrip_method);
%     end
    
% Do destriping and get the mass coefficients
    cs_mss(i,:,:) = gmt_destriping(cs_tmp,destrip_method);

end
pause(0.5);

% figure(2)
% ttt1(:,:)=cs_res(17,:,:);
% ttt2=gmt_gc2mc(ttt1);
% plot(ttt2(:,2),'-s')
% 
% hold on
% ttt(:,:)=cs_destrip(17,:,:);
% plot(ttt(:,2),'-r')
% ----------------------------------------
% Geoid to Mass
% ----------------------------------------
% for i=1:num_file
%     waitbar(0.7+0.1*i/num_file,hwait,'Convert Gravity Geoid to Equivalent Water Height (unit:meter)');
%     cs_tmp(:,:) = cs_destrip(i,:,:);
%     cs_mss(i,:,:) = gmt_gc2mc(cs_tmp);
% end
% pause(0.5);

% figure(2)
% ttt1(:,:)=cs_res(17,:,:);
% ttt2=gmt_gc2mc(ttt1);
% plot(ttt2(:,2),'-s')
% 
% hold on
% ttt(:,:)=cs_mss(17,:,:);
% plot(ttt(:,2),'-r')


% disp('. remove Uplift PGR');
% grid_pgr  = readpgr_new( 'GIA_n100_uplift_0km.txt');

% ----------------------------------------
% Remove GIA effect
% ----------------------------------------
waitbar(0.8,hwait,'Remove the GIA effect');
if strcmp(option_gia,'GIA_Removed_Geru')
    % GIA file must be in Matlab Path
    grid_pgr  = gmt_readpgr('GIA_n100_mass_0km.txt'); 
    cs_pgr = gmt_grid2cs(grid_pgr,lmax);
    for ii=1:num_file
        cs_tmp(:,:)         = cs_mss(ii,:,:);
        % remove PGR from stokes coefficients
        cs_tmp              = cs_tmp - (time(ii)-mean(time))*cs_pgr;
        cs_mss(ii,:,:)      = cs_tmp(:,:);
    end
end
pause(0.5);

% ----------------------------------------
% Gaussian filtering
% ----------------------------------------        
waitbar(0.9,hwait,'Do Gaussian filtering');
cs_fltr=gmt_gaussian_filter(cs_mss,radius_filter);
pause(0.5);

% ----------------------------------------
% Save the output files
% ----------------------------------------
% from number to string
for ii=1:num_file
    str_year{ii}     = num2str(int_year(ii));
    str_month{ii}   = num2str(int_month(ii),'%02u');
end
% Save SH coefficients
if strcmp(OutputFileFormat_1,'SH_MAT')
    cs_grace=cs_fltr(:,1:OutputFileFormat_2+1,1:OutputFileFormat_2+1);
    save(strcat(dir_out, OutputFileFormat_3,'.mat'),'cs_grace','str_year','str_month','time');
    waitbar(1,hwait,'Save the Files');
    close(hwait);
    return;
end
% Create Grids
for ii=1:num_file
    cs_tmp(:,:)         = cs_fltr(ii,:,:);
    if OutputFileFormat_4==1
        grid_tmp=gmt_cs2grid(cs_tmp,0,1);
    elseif OutputFileFormat_4==0.25
        grid_tmp=gmt_cs2grid(cs_tmp,0,0.25);
    end
    grid_data_grace(:,:,ii) = grid_tmp(:,:);
end
% Save Grids
if strcmp(OutputFileFormat_1,'GRID_MAT')    
    waitbar(1,hwait,'Save the grid file in MAT format');
    save(strcat(dir_out, OutputFileFormat_3,'.mat'),'grid_data_grace','str_year','str_month','time');
    close(hwait);
    return;
end
% if strcmp(OutputFileFormat_1,'GRID_GRD')
%     for ii=1:num_file
%         waitbar(0.9+0.1*ii/num_file,hwait,'Save the monthly grid files');
%         grid_tmp(:,:)=grid_grace(:,:,ii);
%         if OutputFileFormat_4==1
%             grdwrite2(0.5:359.5,-89.5:89.5,double(flipud(grid_tmp)),strcat(dir_out,'GSM_',str_year{ii},str_month{ii},'.grd'));
%         elseif OutputFileFormat_4==0.25
%             grdwrite2(0:0.25:359.75,-90:0.25:90,double(flipud(grid_tmp)),strcat(dir_out,'GSM_',str_year{ii},str_month{ii},'.grd'));
%         end
%     end
    
    %     write_fid = fopen(strcat(dir_out,'CSR_GSM_',str_year{ii},str_month{ii},'.txt'),'w');
    %     for j=1:180
    %         for i=1:360
    % %             fprintf(write_fid,'%4.1f %4.1f %f',lam(i),90.-theta(j),grid(j,i,ii)-(time(ii)-mean(time))*grid_pgr(j,i) );
    %             fprintf(write_fid,'%4.1f %4.1f %f',lam(i),90.-theta(j),grid(j,i,ii) );
    %             fprintf(write_fid,'\n');
    %         end
    %     end
    %     fclose(write_fid);
    
% end

close(hwait);

基于中科院冯伟研究员的GRACE开源工具包修改
所需要的控制文件也一样。