02-Matlab读取从ICGEM下载的GRACE(-FO)卫星重力数据
一、数据准备
首先将从ICGEM官网（http://icgem.gfz-potsdam.de/home）下载的GRACE(-FO)数据解压存至一个文件夹下：
二、数据格式介绍
从ICGEM下载的GRACE（-FO）文件有着相同的数据格式，其文件名如“GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600.gfc”, 其中“GSM”表示该文件对应的是大地水准面球谐位系数，“02095-2002120”该数据文件对应的数据观测时间是2002年第95天至2002年第120天，“UTCSR”表示的是该文件由CSR机构发布。
打开任意一期数据文件，可发现数据由头文件、球谐系数组成:
头文件介绍了GRACE数据的技术说明文档以及其他的基本信息，此处不再进行详细阐述。红色框中阐述了一些基本常数，如地球半径、地球密度常数、球谐系数文件的最大阶数等。最后是极为重要的Stocks球谐系数，主要包含7列数据：第一列无具体含义、第二列L对应阶数、第三列M对应次数、第四列C对应Cnm球谐系数、第五列S对应Snm球谐系数，第六列sigma C对应Cnm的误差，第7列sigma S对应Snm的误差。
指定读取球谐系数的最大阶数，如maxDegree=60,并指定读取的数据时间范围，如timegap=[200204 202012]，将下载的GRACE（-FO）数据文件夹的路径赋给dir_in，如“dir_in='H:\matlab_text_code\GRACE_and_GRACE_Fo\CSR_06_Grace_and_Grace_Fo'；”，运行[GCSR_cnm,GCSR_snm,GCSR_time,GCSR_time_ym]=hfl_readgsm(dir_in,'ICGEM',maxDegree,timegap);即可获取对应时间段的GRACE（-FO）球谐系数，以及对应的时间。
附上hfl_readgsm.m和hfl_readstocks.m函数：
function [cnm,snm,time,time_ym]=hfl_readgsm(dir_in,type,lmax,timegap)
% Read the gravity filed files
%
% INPUT:
% dir_in full path
%
% OUTPUT:
% cnm&snm spherical harmonic coefficients|c\| and |s\|
% time fraction in the year
% time_ym year+month
%
switch (type)
case 'ICGEM'
files=dir(fullfile(dir_in,'*.gfc'));
lenfiles=length(files);
cnm=zeros(lmax+1,lmax+1,lenfiles);snm=cnm;time_bounds=zeros(lenfiles,2);
h = waitbar(0,'Loading the GRCAE data...');
for i=1:lenfiles
filename=fullfile(dir_in,files(i).name);
time_bounds(i,:)=[str2double(files(i).name(7:13)) str2double(files(i).name(15:21))];
[cnm(:,:,i),snm(:,:,i)]=hfl_readstocks(filename,lmax);
waitbar(i/lenfiles,h,['loading' num2str(i) '/ ' num2str(lenfiles)]);
end
close(h)
[time,time_ym]=gmt_TimeBoundsToGRACETime(time_bounds);
pos=(time_ym>=timegap(1)&time_ym<=timegap(2));
cnm=cnm(:,:,pos);
snm=snm(:,:,pos);
time=time(pos);
time_ym=time_ym(pos);
otherwise
warndlg('Format of filetype is only ICGEM!','Warning');
end
end
function [cnm,snm]=hfl_readstocks(filename,lmax)
% Read stocks potential coeffcients from file.
%
% Input:
% filename: file of potential coeffcients in ICGEM format (*.gfc).
% Output:
% cnm,snm: matricies containing the potential coefficients.
% The dimensions are (n + 1) x (n + 1) with n = lmax.
% The lower triangular matrix elements cnm(n + 1, m + 1) and snm(n + 1, m + 1)
% contain the potential coefficients of degree n and order m.
% open the file
fid = fopen(filename);
hasErrors = true;
% read header
while 1
line = fgetl(fid);
if ~ischar(line), break, end
if isempty(line), continue, end
keyword = textscan(line,'%s',1);
if(strcmp(keyword{1}, 'max_degree'))
cells = textscan(line,'%s%d',1);
maxDegree = cells{2};
end
if(strcmp(keyword{1}, 'errors'))
cells = textscan(line,'%s%s',1);
if(strcmp(cells{2}, 'no'))
hasErrors = false;
end
end
if(strcmp(keyword{1}, 'end_of_head'))
break
end
end
if (lmax==maxDegree)
% init the output matricies.
cnm = zeros(lmax+1, lmax+1);
snm = zeros(lmax+1, lmax+1);
% read potential coefficients
while 1
line = fgetl(fid);
if ~ischar(line), break, end
if isempty(line), continue, end
if(hasErrors)
cells = textscan(line,'%s%d%d%f%f%f%f');
else
cells = textscan(line,'%s%d%d%f%f');
end
if(cells{2}+1<=lmax+1&&cells{3}+1<=lmax+1)
cnm(cells{2}+1, cells{3}+1) = cells{4};
snm(cells{2}+1, cells{3}+1) = cells{5};
else
continue;
end
end
else
warndlg('Format of file is wrong, please check the maxDegree!','Warning');
end
fclose(fid);
end
