Grace Mascon 数据处理

---

created: 2024-02-10T11:19
updated: 2024-02-11T10:55

内容主要参考
zhuanlan.zhihu.com/p/599446064

官网所下载的 nc 文件是全部的数据内容，包含了目前所以的年份的相关数据（因为没有单独的进行选择日期）

matlab直接ncread就行啦，arcgis可以读取nc文件

补充介绍文献部分

张岚，孙文科. 2022. 重力卫星 GRACE Mascon 产品的应用研究进展与展望. 地球与行星物理论评，53（1）：35-52. doi：10.19975/j.dqyxx.2021- 033

Zhang L, Sun W K. 2022. Progress and prospect of GRACE Mascon product and its application. Reviews of Geophysics and Planetary Physics, 53(1): 35-52. doi:10.19975/j.dqyxx.2021-033

通过文献阅读可得 Mascon 产品无需进行预处理即可进行计算处理可得
来自文献 张岚，2022 提到“GRACE 观测数据主要是以球谐系数的形式给出，需要应用者进行一系列预处理才可以得到对应的物理量. 为了克服此困难，也为了提高 GRACE 恢复重力场地空间分辨率，相关机构在近些年推出了新一代 GRACE 观测数据产品，即 Mascon 产品. 该产品的初衷是便于非大地测量和地球物理专业的人使用，比如水文学家、海洋学家，它无需进行任何后处理过程，使用上更加方便.”

代码部分

matlab 代码主要参照，经过适当改变可以进行使用
https://zhuanlan.zhihu.com/p/650622267

数据计算公式如下
CORRECTED_GRACE_MASCON =
GSU - MASCON_C20 + SLR_C20 + DEG1 - GIA + GAD

grace 卫星的相关数据官网

https://www2.csr.utexas.edu/grace/RL06_mascons.html

The GRACE and GRACE-FO solution time series are provided as single NetCDF timeseries. The CSR RL06.2 mascons are available in two formats.

官网提供的可视化工具
https://colab.research.google.com/drive/146ZbgEAJEOM5UxC2fnsraGm6DxPSMB2B?usp=sharing

CORRECTED_GRACE_MASCON 是经过计算过的，但是也提供了没有进行校正的部分，如图所示

matlab代码运行结果

代码内容为：

clc
clear
%% 读取识别
outpath='.';
% file='CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.nc';
file='CSR_GRACE_GRACE-FO_RL0602_Mascons_all-corrections.nc';
ncdisp(file) % 这里是展示输出所有的值

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 这里先读取是为了获取一些基本信息，比如 time_epoch 等初始信息

%% 数据处理 逐个输出
% 读取NetCDF文件
% ncfile = 'CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.nc';
ncfile = 'CSR_GRACE_GRACE-FO_RL0602_Mascons_all-corrections.nc';

% 这里是单独读取其中的 variables 结果
lat = ncread(ncfile, 'lat');
lon = ncread(ncfile, 'lon');
time = ncread(ncfile, 'time'); % 获取了相关时间
data = ncread(ncfile, 'lwe_thickness');

% 创建输出目录
output_dir = 'nctif';
if ~exist(output_dir, 'dir') % 判断是否是路径  exist(name, 'type')
    mkdir(output_dir);
end

% 循环遍历每个月的数据，并将其输出为tif格式
for i = 1:length(time)
    % 从时间戳中获取年份和月份
    date_num = datenum(2002, 1, 1) + time(i);
    % date_num = datetime(2002, 1, 1) + time(i);
    date_num=double(date_num);
    year = datestr(date_num, 'yyyy');
    month = datestr(date_num, 'mm');
    % year = datetime(date_num, 'yyyy');
    % month = datetime(date_num, 'mm');

    % 创建输出文件名
    output_filename = fullfile(output_dir, sprintf('%s_%s', year, month));

    % 获取当前月份的数据
    month_data = data(:,:,i);
    month_data=flipud(rot90(month_data,1));   %镜像反转，不反转的话最后的图像的南北朝向是错的
    % flipud - 将数组从上向下翻转
    %逆时针旋转90° 

    % 创建geotiff文件
    % 无非就是将矩阵输出为带有地理信息的 tiff 数据
    R=georasterref('RasterSize',size(month_data),'Latlim',[double(min(lat)) double(max(lat))],'Lonlim',[double(min(lon)) double(max(lon))]);
    % 写入真正的 TIf
    geotiffwrite([outpath, '\',output_filename,'.tif'],month_data,R) % 注意 R 的调用位置
end

disp("批量处理完成") 

补充知识
TIFF（Tagged Image File Format）是一种灵活的图像文件格式，它支持在文件中存储各种元数据和参考信息，不仅限于地理信息。

arcgis 数据处理

可以利用 arcgis pro 中的 根据NetCDF栅格图层来导入数据，注意输入 NetCDF 文件的路径导入

（必须要根据 netCDF （应该是调用相关的函数）来进行了相关操作）

由于波段具有224，所以 RGB 通道进行显示的时候需要人为的选择相关的波段进行显示