【GDAL-Python】6-在Python中使用GDAL合并栅格数据

1-介绍

（1）教程内容：介绍如何使用 Python 中的 GDAL 将同一目录中的一堆栅格切片合并到单个 geotiff 文件中；

（2）技术路线：使用GDAL合并栅格数据有两个选项：
1） 在Python脚本中使用子进程调用gdal_merge.py （终端也可以）
2） 创建一个虚拟栅格vrt，然后使用translate转换为geotiff
PS:在合并影像文件较多且数据量大时推荐2），1）方法会加载所有文件，对硬件的内存占用较高。

1.1 虚拟栅格vrt介绍

虚拟栅格是一种将一堆现有数据集组织到一个文件中的方法，也就是告诉影像文件所在位置，以merge.vrt文件为例，内存2.48k，里面内容如下所示

<VRTDataset rasterXSize="5024" rasterYSize="4851">
  <SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>
  <GeoTransform> -8.1622514831999993e+01,  3.3507924922972912e-06,  0.0000000000000000e+00,  3.0274203637999999e+01,  0.0000000000000000e+00, -3.3507923440937610e-06</GeoTransform>
  <VRTRasterBand dataType="UInt16" band="1">
    <NoDataValue>65535</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">dem00.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="2513" RasterYSize="2427" DataType="UInt16" BlockXSize="2513" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2513" ySize="2427" />
      <DstRect xOff="0" yOff="2413.02479225609" xSize="2512.7121477539" ySize="2426.72185112651" />
      <NODATA>65535</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">dem01.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="2513" RasterYSize="2423" DataType="UInt16" BlockXSize="2513" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2513" ySize="2423" />
      <DstRect xOff="111.667314776095" yOff="0" xSize="2512.12721150489" ySize="2422.15845285225" />
      <NODATA>65535</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">dem10.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="2512" RasterYSize="2428" DataType="UInt16" BlockXSize="2512" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2512" ySize="2427.61982265555" />
      <DstRect xOff="2400.35663756777" yOff="2422.53508019045" xSize="2512.87449740808" ySize="2428.46491980955" />
      <NODATA>65535</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">dem11.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="2512" RasterYSize="2424" DataType="UInt16" BlockXSize="2512" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2512" ySize="2424" />
      <DstRect xOff="2511.52168310533" yOff="7.64535589435378" xSize="2512.28568147848" ySize="2424.27586248922" />
      <NODATA>65535</NODATA>
    </ComplexSource>
  </VRTRasterBand>
</VRTDataset>

1.2 合并要点

合并文件需要满足：（1）合并文件不支持不含地理参考的文件；（2）相同的波段数

地理参考：遥感影像的地理参考是指为遥感影像设置的空间坐标系统和投影参数，这些信息使得影像能够在地球表面上正确地定位。地理参考包括了GeoTransform和Projection两个主要参数。
1）GeoTransform是一个六元素元组，包含了左上角的x坐标、水平分辨率、旋转角度等信息，用于描述影像在水平方向上的空间分布特性；
2）Projection参数则定义了影像的空间坐标系统，如WGS_1984等，这决定了影像如何在地球表面进行投影

2-代码实现

2.1 数据介绍

采用的数据是5-在Python中使用GDAL将栅格数据分割成相等的部分均匀分割的4块。

2.2 代码实现

from osgeo import gdal
import glob
import subprocess
import os

demList=glob.glob('dem[0-9][0-9].tif')
print(demList)

# =============================1 利用gdal.Translate合并========================
vrt=gdal.BuildVRT('merged.vrt',demList)
gdal.Translate('mergedDEM2.tif',vrt)

# =============================2.利用gdal_merge.py========================
#由于gdal_merge.py不是可执行文件，而是脚本，可以在终端运行
# cmd = "gdal_merge.py -ps 10 -10 -o mergedDEM.tif"
# subprocess.call(cmd.split()+demList)

vrt=None

2.3 效果显示

效果如下，其在接边处还有待处理。（TODO)

3.参考资料

（1）gdal应用-gdalbuildvrt
（2）Python GDAL工具使用及使用VRT格式数据处理
（3）gdal_merge使用说明