如何在data.vrt中引用lat/lon.vrt并写入元数据geoarrays？

嘿，刚好之前折腾过GDAL这种自定义地理参考的场景，我给你一步步捋清楚怎么操作～

你需要用文本编辑器（比如Notepad++、VS Code或者终端里的vim）打开data.vrt文件，做两处关键修改：

1. 添加GeoLocation标签 #

找到文件里的根元素<VRTDataset>，在它内部（比如放在所有<VRTRasterBand>标签的前面或者后面都可以）插入这段XML代码，用来指定经度和纬度对应的VRT文件：

<GeoLocation>
  <XBand>1</XBand>
  <YBand>1</YBand>
  <XDataset>lon.vrt</XDataset>
  <YDataset>lat.vrt</YDataset>
</GeoLocation>

这里的XBand和YBand填1就行，因为你的lon.tif和lat.tif都是单波段的栅格；XDataset和YDataset要填对路径——如果lat.vrt、lon.vrt和data.vrt在同一个文件夹里，直接填文件名就行，不然就得写绝对路径。

2. 写入GeoArrays元数据 #

接下来在<VRTDataset>里找<Metadata>部分（如果没有的话，就自己新建一个），添加GDAL专用的地理数组元数据，确保软件能正确识别经纬度和数据栅格的对应关系：

<Metadata domain="GDAL_GEOARRAYS">
  <MDI key="X_DATASET">lon.vrt</MDI>
  <MDI key="X_BAND">1</MDI>
  <MDI key="Y_DATASET">lat.vrt</MDI>
  <MDI key="Y_BAND">1</MDI>
  <MDI key="PIXEL_OFFSET">0</MDI>
  <MDI key="LINE_OFFSET">0</MDI>
  <MDI key="PIXEL_STEP">1</MDI>
  <MDI key="LINE_STEP">1</MDI>
</Metadata>

给你解释下这些参数：

• PIXEL_OFFSET/LINE_OFFSET：坐标数组的起始像素和行号，一般都是0（从栅格左上角开始算）
• PIXEL_STEP/LINE_STEP：坐标数组的采样间隔——因为你的lat/lon栅格和data栅格是一一对应的，所以填1就行；如果坐标栅格分辨率更低，比如每2个data像素对应1个坐标点，那这里就填2。

保存修改后的data.vrt，然后在终端里跑这个命令验证：

gdalinfo data.vrt

如果一切正常，你会在输出里看到类似这样的内容：

GeoLocation:
X: lon.vrt band 1
Y: lat.vrt band 1
同时在Metadata区域也能看到GDAL_GEOARRAYS的相关条目，这就说明配置成功了。

如果你的lat/lon栅格不是WGS84（EPSG:4326）投影，记得在lat.vrt和lon.vrt里也加上投影信息——比如在它们的<Metadata>里插入这段：

<Metadata domain="GEOLOCATION">
  <MDI key="CRS">EPSG:XXXX</MDI>
</Metadata>

把XXXX换成你实际的EPSG代码就行。

内容的提问来源于stack exchange，提问作者Ashley A.