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 通道进行显示的时候需要人为的选择相关的波段进行显示