利用GMT绘制逐月的GRACE趋势堆叠图

1.matlab中批量导出数据为txt文件，这里我导出的是CSR mascon 产品，会生成每个月的数据以及缺失月份的label.txt文件。

add = '/Volumes/Expansion/UCAS/';
data = mascon_csr;lon = double(data.lon);  % 720×1440
lat = double(data.lat);
rg = double(data.rg);    % 720×1440×243
tt = data.tt;            % 243×3，列为 [decimal_year, year, month][LonGrid, LatGrid] = meshgrid(lon(1, :), lat(:, 1));
LonVec = LonGrid(:);
LatVec = LatGrid(:);
N = numel(LonVec);% 构造唯一时间键
time_key = tt(:, 2) * 100 + tt(:, 3);  % e.g., 200301% 创建记录缺失的 cell 数组
missing_labels = {};% 初始化时间计数器
counter = 1;for y = 2002:2024for m = 1:12key = y * 100 + m;idx = find(time_key == key);if ~isempty(idx)rg_i = rg(:, :, idx(1));  % 只取第一个月记录，避免多帧拼接失败rgVec = rg_i(:);elsergVec = zeros(N, 1);% 记录缺失标签及顺序编号missing_labels{end+1, 1} = sprintf('%04d-%02d', y, m);missing_labels{end, 2} = counter;endout = [LonVec, LatVec, rgVec];filename = sprintf('data_%04d_%02d.txt', y, m);filepath = fullfile(add, filename);writematrix(out, filepath, 'Delimiter', 'space');fprintf('Saved: %s (Index: %d)\n', filename, counter);counter = counter + 1;end
end% 写入缺失月份标签文件
if ~isempty(missing_labels)label_path = fullfile(add, 'label.txt');fid = fopen(label_path, 'w');fprintf(fid, 'Missing Month\tIndex\n');for i = 1:size(missing_labels, 1)fprintf(fid, '%s\t%d\n', missing_labels{i,1}, missing_labels{i,2});endfclose(fid);fprintf('缺失标签已保存至：%s\n', label_path);
elsefprintf('无缺失月份，未生成 label.txt。\n');
end

2.在GMT中绘图

% Chistrong Wen
% University of Stuttgart
% 2025-6-18
gmt begin PhD_figure_gnss png,pdf E600
# 绘制底图
gmt set FORMAT_GEO_MAP = ddd:mm:ssF
# 推荐使用以下设置（兼容 GMT6）：
gmt set MAP_FRAME_TYPE plain
gmt set MAP_FRAME_PEN 0.8p
gmt set FONT_ANNOT_PRIMARY 10p,Helvetica,black
gmt set FONT_LABEL 10p,Helvetica,black
gmt set MAP_TICK_PEN 0.75p
gmt set MAP_TICK_LENGTH_PRIMARY -0.15c
gmt set MAP_ANNOT_OFFSET_PRIMARY 0.1c
gmt set MAP_ANNOT_OFFSET_PRIMARY = -0.15c
gmt set MAP_TICK_LENGTH_PRIMARY = -0.15c
gmt set FORMAT_ANNOT_PRIMARY = 2p,Times-Roman,black# 色标
gmt makecpt -Cpolar -T-10/10/1# 布局参数
rows=23
cols=12# 读取 label.txt 缺失编号（转为 0 起始）
missing_indices=$(awk 'NR > 1 {print $2 - 1}' label.txt)# 定义判断函数
is_missing() {for miss in $missing_indices; doif [[ $1 -eq $miss ]]; thenreturn 0fidonereturn 1
}gmt subplot begin ${rows}x${cols} -Fs6c/6c -R-180/180/-90/90 -JN6c -Baf -M0.3cfor i in $(seq 0 275); dogmt subplot set $i# 如果该编号在缺失列表中 → 跳过绘图if is_missing $i; thencontinuefiyear=$((2002 + i / 12))month=$((1 + i % 12))fname=$(printf "data_%04d_%02d.txt" $year $month)if [[ -f "$fname" ]]; thengmt xyz2grd "$fname" -R0/360/-90/90 -I1.0 -Gtmp_${i}.grdgmt grdimage tmp_${i}.grd -R-180/180/-90/90 -JN6c -Cgmt coast -A500 -W0.1p,blackfi
donegmt subplot end
gmt colorbar -DjBC+w36c/0.5c+o0c/-1.5c+h -Bxa5f5 -By+l"EWH (cm)" 
gmt endrm -f tmp_*.grd

运行结果：