【最新】冯伟老师GRACE工具箱改进版（无GUI界面）-可读取RL06数据

看到网上很多的同学和粉丝都提到了冯伟老师GRACE程序包的数据读取问题。趁着疫情隔离期间，花了一天时间完成了修改。这里我将详细介绍具体的代码调试过程。第一个专栏是如何读取自带的RL05数据，目前数据已经发布到RL06了，因此我将在这里介绍本人改进的可以读取RL06数据的一套代码。其实大多数函数都没有改变，只是修改了读取GSM数据、一阶项、二阶项的函数。下载地址：

链接：https://pan.baidu.com/s/1nWtM3tSA-pgdP4FUOCSTqg
提取码：UCAS

我们下载冯伟老师的工具包，我在这里把我调好的可以读取RL05数据的程序放在以下的链接，感兴趣的可以下载调试，后续会上传可以读取最新的RL06数据的程序包。 

1. 前期准备：

GRACE数据下载：http://icgem.gfz-potsdam.de/series/01_GRACE/CSR/CSR%20Release%2005
一阶项下载地址：https://podaac­w10n.jpl.nasa.gov/
二阶项下载地址：http://download.csr.utexas.edu/pub/slr/degree_2/

接着放入对应的目录下，C:\GRACE_Matlab_Toolbox\GRACE_data\RL06

接着放入对应的目录下，C:\GRACE_Matlab_Toolbox\GRACE_data\RL06

 2. 控制文件的修改

相比于RL05数据的控制文件，读取ICGEM网站的CSR RL06数据需要对应修改如下：

注意：控制文件中需要读取目录下所有GSM的文件名，手动修改是不可能的，可以采用以下的方法：在文件对应得目录下，新建txt文件，里面写入

DIR *.*/B>LIST.TXT

并将后缀名改为.bat，运行即可得到所有文件得txt文件，复制粘贴即可。

需要注意的是，我没有采用具有GUI界面的那个程序。下面切入正题：直接运行即可。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         A demo for GRACE toolbox                        %
%                         2022-11-29 Chistrong Wen                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
controlfile_path = 'GRAMAT_Control_File_csr_swenson.txt';
% GRACE_Matlab_Toolbox_preprocessing_core(controlfile_path)
%
load GRACE_resultsgrace_csr_2002_2017.mat
lon = 0.5:1:359.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);
A.lon = lon;
A.lat = lat;
A.rg  = flipud(grid_data_grace(:,:,47)).*100; % convert m to cm
wzq_plot(A)
xlabel('longitude', 'Fontname', 'Times New Roman', 'Fontsize',13);
ylabel('latitude', 'Fontname', 'Times New Roman', 'Fontsize',13);
caxis([-40,40])
h = colorbar
h.Title.String = 'cm';colormap jet
cb2.Title.FontName = 'times new roman';
cb2.Title.FontSize = 12;

最主要的函数修改如下：

（1）读取GRACE GSM数据

function    [cs,cs_sigma,int_year,int_month,meanday,time] = gmt_readgfc_ucas(pathname)
% Read the ICGEM gravity filed files
%
% INPUT:
%   pathname    filename with full path
%
% OUTPUT:
%   cs          spherical harmonic coefficients in CS-format
%   cs_sigma    spherical harmonic coefficients in CS-format (formal error)
%
% FENG Wei 24/03/2015
% State Key Laboratory of Geodesy and Earth's Dynamics
% Institute of Geodesy and Geophysics, Chinese Academy of Sciences
% fengwei@whigg.ac.cn

[dir_in,file_name,file_type]=fileparts(pathname);
if  ~strcmp(file_type,'.gfc') && ~strcmp(file_type,'.GFC')
    error('Format of gmt_readgfc function is wrong!');
end
% read the header
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);
    % In JPL GSM files, the maximum degree changes, 60 or 90
    if size(tline,2)>10 && strcmp(tline(1:10),'max_degree')
        degree_max=str2double(tline(11:end));
    end
end
fclose(fid);
% re-read the file, skip the comment lines
[~, l, m, Clm, Slm, Clm_sigma, Slm_sigma] = textread(pathname,'%s %u %u %f %f %f %f','headerlines',head_index+1);
for i = 1:length(l)
    sc_tmp(l(i)+1,degree_max+1-m(i)) = Slm(i);
    sc_tmp(l(i)+1,degree_max+1+m(i)) = Clm(i);
    sc_sigma_tmp(l(i)+1,degree_max+1-m(i)) = Slm_sigma(i);
    sc_sigma_tmp(l(i)+1,degree_max+1+m(i)) = Clm_sigma(i);
end

cs=gmt_sc2cs(sc_tmp);
cs_sigma=gmt_sc2cs(sc_sigma_tmp);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% modified 2022-11-29 GET TIME
% Get time tag GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600.gfc
year1 = str2num(file_name(7:10));
year2 = str2num(file_name(15:18));
day1  = str2num(file_name(11:13));
day2  = str2num(file_name(19:21));
if year1 == year2
    meanday = (day1+day2)/2;
else
    if (day1+(366-day1+day2)/2)>365 % in latter year
        year1   = year1 + 1;
        meanday = day2-(366-day1+day2)/2;
    else
        meanday = day1+(366-day1+day2)/2;  % in former year
    end
end
time        = year1 + meanday/365.;
meanday     = round(meanday);
int_year    = year1;
int_month   = gmt_get_mon_day(year1,meanday+1);
end

（2）读取一阶项

function [ cs_replace,tag] = gmt_replace_degree_1(dir_in,cs,int_year,int_month,num_file )

% Read the degree 1 files and replace the corresponding original Stokes
% coefficients with them
%
% INPUT:
%   dir_in      full path
%   cs          spherical harmonic coefficients in CS-format
%   int_year    year
%   int_month   mont
%   num_file    number of files
%
% OUTPUT:
%   cs_replace  spherical harmonic coefficients with replaced degree one
%   tag         check
%
% FENG Wei 22/03/2015
% State Key Laboratory of Geodesy and Earth's Dynamics
% Institute of Geodesy and Geophysics, Chinese Academy of Sciences
% fengwei@whigg.ac.cn
cs_replace=zeros(size(cs));
[~, FILE_NAME,~]=fileparts(dir_in);
ind = 0;
if (strcmp(FILE_NAME,'TN-13_GEOC_CSR_RL06'))
    tag=1;
    % read the header
    %---------------------------modified-----------------------------------
    head_index=0;
    fid2 = fopen(dir_in,'r');
    while ~feof(fid2)
        str = fgetl(fid2);
        if fid2 == -1
            ('Error opening the file');
        end
        if (isempty(str) | str==' ')
           continue
        end % to find the blankspace, and skip it!
        if ~ischar(str)
           break,
        end
        if strcmp(str(1:6),'GRCOF2')
            ind = ind+1;
            a=sscanf(str,'%s %d %d %f %f %f %f %f %f %f');

        if(mod(ind,2)==0)
            D_C(ind,2) = a(9);
            D_S(ind,2) = a(10);
            D_C(ind,4) = a(11);
            D_S(ind,3) = a(12);
            l(ind) = a(2);
        else
            D_C(ind,1) = a(9);
            D_S(ind,1) = a(10);
            D_C(ind,3) = a(11);
            D_S(ind,3) = a(12);
            l(ind) = a(2);
        end
        end
    end
    % Replace Degree 1
    for ii=1:num_file
        cs_replace(ii,:,:) = cs(ii,:,:);
        for jj=1:length(l)
            if (int_year(ii)==year(jj) && int_month(ii)==mon(jj) && l(jj)==1 && m(jj)==0)
                cs_replace(ii,2,1) = aa(jj);
            elseif (int_year(ii)==year(jj) && int_month(ii)==mon(jj) && l(jj)==1 && m(jj)==1)
                cs_replace(ii,2,2) = aa(jj);
                cs_replace(ii,1,2) = bb(jj);
            end
        end
    end
else
    tag=0;
    %     warndlg('Format of degree 1 file is wrong!','Warning');
end

（3）读取二阶项，C20没有改变，只是修改了C21S21C22S22的函数。

function [cs_replace,tag,FILE_NAME] = gmt_replace_C21_S21_C22_S22(dir_in,cs,int_year,int_month,num_file)

% Read the C21, S21, C22, S22 files and replace the corresponding original
% Stokes coefficients with them
%
% Read C21_S21_RL05.txt or C22_S22_RL05.txt file provided from CSR's FTP site ftp://ftp.csr.utexas.edu/pub/slr/degree_2/
% Replace C21 & S21 or C22 & S22 in GRACE GSM files with those in
% C21_S21_RL05.txt or C22_S22_RL05.txt
% Restored AOD should be removed from values in C21_S21_RL05.txt and C22_S22_RL05.txt
%
% INPUT:
%   dir_in      full path
%   cs          spherical harmonic coefficients in CS-format
%   int_year    year
%   int_month   mont
%   num_file    number of files
%
% OUTPUT:
%   cs_replace  spherical harmonic coefficients with replaced C21, S21, C22, S22
%   tag         check
%

%
% FENG Wei 22/03/2015
% State Key Laboratory of Geodesy and Earth's Dynamics
% Institute of Geodesy and Geophysics, Chinese Academy of Sciences
% fengwei@whigg.ac.cn

cs_replace=cs;
time=zeros(num_file,1);
[~, FILE_NAME,~]=fileparts(dir_in)

if strcmp(FILE_NAME,'C21_S21_RL06') || strcmp(FILE_NAME,'C22_S22_RL06')
    tag=1;
    % read the header
    ind = 0;
    fid2 = fopen(dir_in,'r');
%     fid2 = fopen(file_name, 'r');
    while ~feof(fid2)
        str = fgetl(fid2);
        if fid2 == -1
        ('Error opening the file');
        end
        if (isempty(str) | str==' '|str == '#')
        continue
        end % to find the blankspace, and skip it!
        if ~ischar(str)
        break,
        end
        if strcmp(str(2),'2')
            ind = ind+1;
            a=sscanf(str,' %f %f %f %f %f %f %f %f %f');
            time_year(ind) = a(1); % time
            C21_C22(ind) = a(2); % C21 from SLR (normalized)
            S21_S22(ind) = a(3); % S21 from SLR (normalized)
%             c21(ind,4) = a(4)*1e-10; % Solution sigma for C21 (1E-10)
%             c21(ind,5) = a(5)*1e-10; % Solution sigma for S21 (1E-10)
            C21_C22_aod(ind) = a(6)*1e-10; % Solution sigma for C21 (1E-10)
            S21_S22_aod(ind) = a(7)*1e-10; % Solution sigma for S21 (1E-10)

        end
    end

    % re-read the file, skip the comment lines
%     [time_year, C21_C22, S21_S22, ~, ~, C21_C22_aod, S21_S22_aod ] = textread(dir_in,'%f %f %f %f %f %f %f ',head_index)
    C21_C22_size=max(size(C21_C22));
    for ii=1:num_file
        % Approximate mid-point of monthly solution, See C21_S21_RL05.txt
        time(ii) = int_year(ii) + gmt_get_days(int_year(ii),int_month(ii),15)/365;
        for jj=1:C21_C22_size
            if abs(time(ii)-time_year(jj))<=0.03
                if strcmp(FILE_NAME,'C21_S21_RL06')
                    cs_replace(ii,3,2) = C21_C22(jj)-C21_C22_aod(jj)*(1E-10);% restored AOD has to be removed;
                    cs_replace(ii,1,3) = S21_S22(jj)-S21_S22_aod(jj)*(1E-10);
                elseif strcmp(FILE_NAME,'C22_S22_RL06')
                    cs_replace(ii,3,3) = C21_C22(jj)-C21_C22_aod(jj)*(1E-10);% restored AOD has to be removed
                    cs_replace(ii,2,3) = S21_S22(jj)-S21_S22_aod(jj)*(1E-10);
                end
            end
        end
    end
else
%     warndlg('The name of C21&S21 or C22&S22 file is wrong!','ERROR');
    tag=0;
end
end

希望对大家有帮助！有问题，欢迎加入交流群