将xyz格式的GRACE数据转成geotiff格式

我们需要将xyz格式的文件转成geotiff便于成图，或者geotiff转成xyz用于数据运算，下面介绍如何实现这一操作，采用GMT和matlab两种方法。

1.GMT转换

我们先准备一个xyz文件，这里是一个降水文件。在gmt中采用以下的语句实现xyz转grd网格文件

xyz2grd DJF.txt -R-180/180/-90/90 -JN12c -I5 -Gm1.grd gmt grdsample m1.grd -Gm11.grd -I0.5

语句实现了将xyz转成grd文件，然后在Global mapper中打开，导出为geotiff格式即可，然后可以用于绘图和地理信息处理显示。

2.matlab程序实现

我们现在matlab中读取一个geotiff文件以观察其组成要素。

[A1,R] = geotiffread('Tibet_1000.tif');A1(A1<=0)=0;

可以看到，其由两部分组成，一个是数值矩阵（A1），另一个是投影信息（R）。因此我们如果想要将xyz数据转成geotiff格式的文件，同样需要准备两个信息：一个是数据矩阵，另一个是投影信息（当然需要自己设置）。

[A1,R] = geotiffread('Tibet_1000.tif');A1(A1<=0)=0; lon = linspace(65.0042,109.9958,5400); lat = linspace(20.0042,49.9958,3600); [lon,lat] = meshgrid(lon,lat);O1.lon = lon;O1.lat = lat;O1.rg = flipud(double(A1)); rg_plot(O1),title('10 km DEM'),colorbar

需要注意的是，一般我们所说的xyz数据是三列的，分别是

如果要转成geotiff文件，需要提前转成数值矩阵，即通常需要reshape一下（这针对于转换的影响为矩形）。

接下来，我们实现一个xyz转成geotiff的例子。

（1）先准备一个xyz数据，这里以MSSA插值的GRACE level03数据集为例，数据参考以下专栏：

分享一个月份连续的MSSA插值的GRACE level03数据集_grace数据集-CSDN博客

matlab先读取xyz数据，然后绘图，可见reshape的正确性：

% % load data A = load('GRACE_MSSA_2022_01.xyz'); %% reshape O.lon = reshape(A(:,1),181,361); O.lat = reshape(A(:,2),181,361); O.rg = reshape(A(:,3),181,361); wzq_plot(O)

wzq_plot函数如下，其中缺失的报错文件参加B站的置顶评论：

绘图函数的使用wzq_plot - 哔哩哔哩 (bilibili.com)

function wzq_plot(wzq) pcolor(wzq.lon,wzq.lat,wzq.rg) shading interp hold on; if max(wzq.lon(:)) < 200 coast=load('coastline-from-GMT-WNI.dat'); else coast=load('coastline-from-GMT-WNI-0-360.dat'); end plot(coast(:,1),coast(:,2),'k') hold off; colorbar colormap jet end

（2）配置投影信息

这里我们借用其他读取的geotiff文件的投影信息，然后按照实际情况进行修改配置，这里我们采用一个DEM的投影信息，我们需要修改的地方包括：

X Y的范围  采样分辨率  经纬度范围  栅格数量

然后我们得到配置好的新的投影信息,但是实际上有更简单的配置方法：

[A1,R] = geotiffread('Tibet_1000.tif');A1(A1<=0)=0; % % load data A = load('GRACE_MSSA_2022_01.xyz'); %% reshape O.lon = reshape(A(:,1),181,361); O.lat = reshape(A(:,2),181,361); O.rg = reshape(A(:,3),181,361); wzq_plot(O) %% setting new projection info R1.RasterInterpretation = 'Postings'; R1.XIntrinsicLimits = [1,361]; R1.YIntrinsicLimits = [1,181]; R1.SampleSpacingInLatitude = 1; R1.SampleSpacingInLongitude = 1; R1.LatitudeLimits = [-90,90]; R1.LongitudeLimits = [0,360]; R1.RasterSize = [181,361]; R1.AngleUnit = 'degree'; R1.ColumnsStartFrom = 'north'; R1.RowsStartFrom = 'east'; R1.CoordinateSystemType = 'geographic'; R1.AngleUnits = 'degrees'; latlim = [-90,90]; lonlim = [0,360]; R0 = georefcells(latlim,lonlim,size(O.rg)); % 设置一个地理坐标 geotiffwrite('GRACE_xyz2tiff.tif', O.rg, R0); %

然后我们可以得到GRACE_xyz2tiff.tif文件，需要注意的是，再次运行前需要删除之前生成的文件。

在global mapper中可以打开