2. GDAL 栅格工具 (GDAL Tools)

来源说明 本文翻译自 Spatial Thoughts，原作者为 Ujaval Gandhi。

**原文链接：** https://courses.spatialthoughts.com/gdal-tools.html

**本文仅供学习交流，版权归原作者所有。**

观看模块视频 ↗

2.1 基础栅格处理 (Basic Raster Processing)

我们将通过处理 SRTM 高程栅格数据来开始学习基本的 GDAL 命令。在 命令提示符 窗口中，使用 cd 命令切换到 srtm 目录，该目录包含珠穆朗玛峰区域周围的 4 个独立 SRTM 图块。

cd srtm

使用 gdalinfo 命令查看单个图像的信息。

gdalinfo N28E086.hgt

一个有用的参数是 -stats，它会计算并显示图像统计信息。在栅格上运行它以获取图像中像素值的一些统计数据。

gdalinfo -stats N28E086.hgt

2.1.1 合并图块 (Merging Tiles)

现在我们将 4 个相邻的 SRTM 图块合并为一个栅格，以便我们可以一起处理它们。GDAL 提供了一种非常有用的格式，称为虚拟栅格 (Virtual Raster)。它允许我们创建一个扩展名为 .vrt 的虚拟文件，该文件是指向多个源文件的指针。.vrt 文件只是一个文本文件，因此它不占用磁盘空间，但允许我们像操作普通栅格文件一样运行任何 GDAL 命令。

首先，我们需要创建一个包含我们要合并的所有文件的文本文件。我们可以使用命令提示符上的 dir 命令列出匹配 *.hgt 模式的文件，并将输出重定向到一个文件。这里 /b 选项以 精简 (Bare) 模式运行命令，排除除文件名以外的所有信息。

Windows

dir /b *.hgt > filelist.txt

对于 Mac/Linux 系统，可以使用 ls 命令实现相同的效果。

Mac/Linux

ls *.hgt > filelist.txt

命令完成后，验证 filelist.txt 是否包含源图块的名称。

现在我们可以使用 gdalbuildvrt 程序从 filelist.txt 中的源文件创建一个虚拟栅格。

gdalbuildvrt -input_file_list filelist.txt merged.vrt

注意： 我们本可以使用命令 gdalbuildvrt merged.vrt *.hgt 一步完成此操作。但是，Windows 上的某些 GDAL 版本无法正确扩展 * 通配符，导致命令出错。建议在 Windows 上使用文件列表而不是通配符，以避免意外结果。[参考]

练习 1

你能找出合并后的栅格中的最高高程值是多少吗？由于这些栅格位于珠穆朗玛峰区域周围，MAXIMUM 值即为峰顶的高程。

2.1.2 格式转换 (Converting Formats)

让我们将 虚拟栅格 转换为 GeoTIFF 文件。gdal_translate 程序允许我们在 GDAL 支持的数百种数据格式之间进行转换。格式通常通过文件扩展名识别。或者，你也可以使用 -of 选项和格式的短名称（如 GTiff）来指定它。

gdal_translate -of GTiff merged.vrt merged.tif

2.1.3 压缩输出 (Compressing Output)

默认输出的 GeoTIFF 文件是未压缩的——这意味着每个像素的值都按原样存储在磁盘上。对于大型栅格，这会消耗大量磁盘空间。明智的做法是使用无损压缩算法来减小栅格大小，同时保持原始数据的完整保真度。GDAL 开箱即用地支持许多压缩算法，可以使用 -co 选项指定。最流行的无损压缩算法是 DEFLATE、LZW 和 PACKBITS。我们可以在数据集上尝试 DEFLATE 算法。

gdal_translate -of GTiff merged.vrt merged.tif -co COMPRESS=DEFLATE

未压缩的文件大小为 100+ MB。应用 DEFLATE 压缩后，文件大小减少到 75MB。我们可以通过指定其他选项进一步减小文件大小。PREDICTOR 选项有助于在相邻值相关时更好地压缩数据。对于高程数据，情况确实如此。TILED 选项将按块而不是按行压缩数据。

注意： 你可以使用续行符将长命令拆分为多行。Windows shell 使用 ^ 作为续行符，而 Mac/Linux shell 使用 \。

Windows

gdal_translate -of GTiff merged.vrt merged.tif ^
  -co COMPRESS=DEFLATE -co TILED=YES -co PREDICTOR=2

Mac/Linux

gdal_translate -of GTiff merged.vrt merged.tif \
  -co COMPRESS=DEFLATE -co TILED=YES -co PREDICTOR=2

现在生成的文件小得多，只有 39MB。

2.1.4 设置 NoData 值 (Setting NoData Values)

gdalinfo 程序的输出显示原始数据的 NoData 值设置为 -32768。我们可以设置一个新的 NoData 值。-a_nodata 选项允许我们指定新值。

Windows

gdal_translate -of GTiff merged.vrt merged.tif ^
  -co COMPRESS=DEFLATE -co TILED=YES -co PREDICTOR=2 -a_nodata -9999

Mac/Linux

gdal_translate -of GTiff merged.vrt merged.tif \
  -co COMPRESS=DEFLATE -co TILED=YES -co PREDICTOR=2 -a_nodata -9999

2.1.5 编写云优化 GeoTIFF (COG)

一种名为云优化 GeoTIFF (COG) 的新格式使访问大量影像变得更加容易和便于分析。云优化 GeoTIFF 的行为就像普通的 GeoTIFF 影像，但你无需在本地下载整个图像，而是可以访问托管在云服务器上的影像的部分，并将其流式传输到 QGIS 等客户端。

GDAL 通过指定 -of COG 选项，使得创建 COG 文件变得非常容易。指定 -co NUM_THREADS=ALL_CPUS 有助于通过使用所有可用的 CPU 进行压缩和创建内部概览图来加速创建过程。

GDAL COG 驱动程序默认启用了以下创建选项 -co：

具有内部切片 (即 -co TILED=YES)
LZW 压缩 (即 -co COMPRESS=LZW)
自动选择预测器 (即 -co PREDICTOR=YES 为数据类型选择合适的预测器)

Windows

gdal_translate -of COG merged.vrt merged_cog.tif ^
  -co COMPRESS=DEFLATE -co PREDICTOR=YES -co NUM_THREADS=ALL_CPUS ^
  -a_nodata -9999

Mac/Linux

gdal_translate -of COG merged.vrt merged_cog.tif \
  -co COMPRESS=DEFLATE -co PREDICTOR=YES -co NUM_THREADS=ALL_CPUS \
  -a_nodata -9999

2.2 处理高程数据 (Processing Elevation Data)

GDAL 附带了 gdaldem 实用程序，它提供了一套用于可视化和分析数字高程模型 (DEM) 的工具。该工具支持以下模式：

Hillshade (山体阴影)
Slope (坡度)
Aspect (坡向)
Color-relief (彩色地形)
Terrain Ruggedness Index (TRI, 地形耐用性指数)
Topographic Position Index (TPI, 地形位置指数)
Roughness (粗糙度)

需要注意的一点是，DEM 的 x、y 和 z 单位应该相同。如果你的数据使用的是地理坐标系（如 EPSG:4326），而高度单位是米，则必须使用 -s 选项指定比例值。

2.2.1 创建山体阴影 (Creating Hillshade)

让我们从合并的 SRTM 数据集创建一个山体阴影图。hillshade 模式创建一个具有漂亮阴影浮雕效果的 8 位栅格。此数据集的 X 和 Y 单位是度，Z 单位是米。所以我们指定 111120 作为比例值。

gdaldem hillshade merged.tif hillshade.tif -s 111120

gdaldem 支持多种山体阴影算法。除了默认算法外，它目前还包括以下算法：

组合阴影 (-combined)：坡度和斜向阴影的组合。
多向阴影 (-multidirectional)：来自 225 度、270 度、315 度和 360 度方位的山体阴影的组合。

让我们使用 -multidirectional 选项创建一个山体阴影图。

gdaldem hillshade merged.tif hillshade_combined.tif -s 111120 -multidirectional

2.2.2 创建彩色地形 (Creating Color Relief)

彩色地形图是一种高程图，其中不同的高程范围被着以不同的颜色。color-relief 模式可以使用文本文件中提供的颜色和高程范围创建彩色地形图。

我们需要使用文本文件提供颜色映射。创建一个名为 colormap.txt 的新文件，内容如下，并将其保存在与 merged.tif 相同的目录中。文本文件的格式是 elevation, red, green, blue 值。

1000,101,146,82
1500,190,202,130
2000,241,225,145
2500,244,200,126
3000,197,147,117
4000,204,169,170
5000,251,238,253
6000,255,255,255

我们可以提供此文件来创建彩色高程图。

gdaldem color-relief merged.tif colormap.txt colorized.tif

Color Relief

练习 2

将彩色地形图保存为 PNG 格式。

2.3 处理航空影像 (Processing Aerial Imagery)

使用 cd 命令切换到 naip 目录，该目录包含 JPEG2000 格式的单个航空影像图块。

cd naip

2.3.1 从源图块创建预览图像

源影像经过高度压缩并覆盖了大片区域。与其在查看器中加载全分辨率图块，不如创建一个预览镶嵌图，帮助我们评估数据的覆盖范围和质量。

我们首先从所有 .jp2 图块创建一个 虚拟栅格 镶嵌图。

Windows

dir /b *.jp2 > filelist.txt

Mac/Linux

ls *.jp2 > filelist.txt

我们正在处理光度 RGB 图块，并希望对结果应用 JPEG 压缩。为了避免 JPEG 伪影并保持 nodata 值，最好添加一个 Alpha 波段，其中包含所有有效像素的掩码。我们使用 -addalpha 选项来创建 Alpha 波段。

gdalbuildvrt -addalpha -input_file_list filelist.txt naip.vrt

我们可以使用 -outsize 选项并指定百分比来生成较小尺寸的预览。下面的命令生成 JPEG 格式的 2% 预览图像。由于我们有一个 4 波段图像，我们使用 -b 选项指定用于生成 JPG 图像的波段。

gdal_translate -b 1 -b 2 -b 3 -of JPEG -outsize 2% 2% naip.vrt naip_preview.jpg

2.3.2 创建图块索引 (Create a Tile Index)

在处理大量影像图块时，生成图块索引非常有用。图块索引层是一个多边形层，包含每个图块的边界框。gdaltindex 程序从输入文件列表创建图块索引。

gdaltindex -write_absolute_path index.shp --optfile filelist.txt

2.3.3 镶嵌并裁剪到 AOI

假设我们要将所有源图块镶嵌成单个图像，并将其裁剪到给定的 AOI（感兴趣区域）。

我们可以使用 gdalwarp 实用程序并使用 -cutline 选项来裁剪栅格。

gdalwarp -cutline aoi.shp  -crop_to_cutline naip.vrt aoi_cut.tif -co COMPRESS=DEFLATE -co TILED=YES

我们还可以向输出文件添加 JPEG 压缩以减小文件大小。请注意，JPEG 是一种有损压缩方法，可能会导致图像镶嵌的边缘出现伪影。为了防止这种情况，我们指定 --config GDAL_TIFF_INTERNAL_MASK YES 选项，该选项使用 Alpha 波段中的掩码。

Windows

gdal_translate aoi_cut.tif aoi.tif ^
  -co COMPRESS=JPEG -co TILED=YES -co PHOTOMETRIC=YCBCR -co JPEG_QUALITY=75 ^
  -b 1 -b 2 -b 3 -mask 4 --config GDAL_TIFF_INTERNAL_MASK YES

Mac/Linux

gdal_translate aoi_cut.tif aoi.tif \
  -co COMPRESS=JPEG -co TILED=YES -co PHOTOMETRIC=YCBCR -co JPEG_QUALITY=75 \
  -b 1 -b 2 -b 3 -mask 4 --config GDAL_TIFF_INTERNAL_MASK YES

2.3.4 创建概览图 (Creating Overviews)

如果你尝试将生成的栅格加载到查看器中，你会发现渲染需要很长时间。这是因为查看器正在以原始分辨率渲染所有像素。常见的解决方案是创建 金字塔图块 或 概览图 (Overviews)。GDAL 提供了 gdaladdo 实用程序来创建概览图。

你可以使用默认选项运行 gdaladdo (GDAL-Add-Overview) 命令来创建内部概览图。

gdaladdo aoi.tif

默认概览图使用最近邻重采样。我们可以尝试使用 -r bilinear 选项进行双线性插值。

gdaladdo -r bilinear --config COMPRESS_OVERVIEW JPEG aoi.tif

2.4 处理卫星影像 (Processing Satellite Imagery)

本节展示如何使用 Landsat-8 卫星数据并创建各种衍生产品。使用 cd 命令切换到 landsat8 目录。

波段号	波段名称
B2	蓝色 (Blue)
B3	绿色 (Green)
B4	红色 (Red)
B5	近红外 (Near Infrared)
B8	全色 (Panchromatic)

cd landsat8

2.4.1 将单个波段合并为 RGB 合成图

让我们通过组合三个不同的波段（红、绿和蓝）来创建一个 RGB 合成图像。这里我们必须使用 -separate 选项，它告诉命令将每个图像放在单独的波段中。

Windows

gdalbuildvrt -o rgb.vrt -separate ^
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif ^
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B3.tif ^
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B2.tif

Mac/Linux

gdalbuildvrt -o rgb.vrt -separate \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B3.tif \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B2.tif

然后使用 gdal_translate 合并它们。

gdal_translate rgb.vrt rgb.tif -co PHOTOMETRIC=RGB -co COMPRESS=DEFLATE

2.4.2 应用直方图拉伸和色彩校正

生成的合成图看起来很暗且对比度低。我们可以应用直方图拉伸来增加对比度。这是使用 -scale 选项完成的。由于大多数像素的值在 0 到 0.3 之间，我们可以选择这些作为最小值和最大值，并应用对比度拉伸使它们从 0 变为 255。

gdal_translate -scale 0 0.3 0 255 -ot Byte rgb.tif rgb_stretch.tif

RGB Composite with Linear Stretch

我们还可以应用非线性拉伸。gdal_translate 有一个 -exponent 选项。选择 0 到 1 之间的指数值将增强低强度值——从而产生更亮的图像。让我们尝试指数值 0.5。

gdal_translate -scale 0 0.3 0 255 -exponent 0.5 -ot Byte rgb.tif rgb_stretch.tif

RGB Composite with Exponential Stretch

2.4.3 栅格代数 (Raster Algebra)

对于栅格代数运算，GDAL 提供了栅格计算器程序 gdal_calc.py。输入栅格使用 A-Z 中的任何字母指定。表达式使用 --calc 选项指定，它支持 NumPy 语法和函数。

下面的示例计算 NDVI (归一化植被指数)。

Windows

python %CONDA_PREFIX%\Scripts\gdal_calc.py ^
  -A RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B5.tif ^
  -B RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif ^
  --outfile ndvi.tif --calc="(A-B)/(A+B)" --NoDataValue=-999

Mac/Linux

gdal_calc.py -A RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B5.tif \
  -B RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif \
  --outfile ndvi.tif --calc="(A-B)/(A+B)" --NoDataValue=-999

练习 3

创建一个包含近红外、红和绿波段的 NRG 合成图像。对结果应用对比度拉伸并将其保存为 PNG 图像。

NRG Composite

2.4.4 全色锐化 (Pan Sharpening)

Landsat-8 卫星在红、绿、蓝波段产生 30m 空间分辨率的图像，而在全色波段产生 15m 空间分辨率的图像。我们可以使用全色波段的更高空间分辨率来提高其他波段的分辨率，从而产生具有更多细节的更清晰图像。这个过程称为 全色锐化。

GDAL 附带了一个脚本 gdal_pansharpen.py，它实现了 Brovey 算法来计算全色锐化输出。

Windows

python %CONDA_PREFIX%\Scripts\gdal_pansharpen.py RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B8.tif ^
  rgb.tif pansharpened.tif -r bilinear -co COMPRESS=DEFLATE -co PHOTOMETRIC=RGB

Mac/Linux

gdal_pansharpen.py RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B8.tif \
  rgb.tif pansharpened.tif -r bilinear -co COMPRESS=DEFLATE -co PHOTOMETRIC=RGB

我们可以应用与之前相同的对比度拉伸并比较输出。你会注意到生成的合成图更加清晰，可以更好地解析场景中的细节。

Windows

gdal_translate -scale 0 0.3 0 255 -exponent 0.5 -ot Byte -a_nodata 0 ^
  pansharpened.tif pansharpened_stretch.tif

Mac/Linux

gdal_translate -scale 0 0.3 0 255 -exponent 0.5 -ot Byte -a_nodata 0 \
  pansharpened.tif pansharpened_stretch.tif

2.5 处理 WMS 图层 (Processing WMS Layers)

GDAL 支持从各种 Web 服务读取数据，包括 Web 地图服务 (WMS) 图层。

2.5.1 列出 WMS 图层

我们可以获取相应服务的 URL 并使用 gdalinfo 列出所有可用图层。

gdalinfo "WMS:https://sedac.ciesin.columbia.edu/geoserver/ows?version=1.3.0"

2.5.2 创建服务描述文件

GDAL 还可以从 WMS 图层创建服务描述 XML 文件。许多 GIS 程序（包括 QGIS）将这些 XML 文件识别为有效的栅格图层。

gdal_translate -of WMS "WMS:https://sedac.ciesin.columbia.edu/geoserver/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=grand-v1%3Agrand-v1-reservoirs-rev01&CRS=CRS:84&BBOX=-153.037,-45.881,176.825,70.398" reservoirs.xml

2.5.3 下载 WMS 图层

如果需要离线访问 WMS 图层，我们可以从 WMS 图层创建地理配准栅格。我们可以使用 gdalwarp 和 -tr 选项来指定离线栅格所需的输出分辨率。

gdalwarp -tr 0.1 0.1 reservoirs.xml reservoirs.tif -co COMPRESS=DEFLATE -co TILED=YES

我们还可以通过使用 -te 选项指定范围来提取特定区域的更高分辨率提取。

gdalwarp -tr 0.005 0.005 -te 68.106 6.762 97.412 37.078 reservoirs.xml reservoirs_india.tif -co COMPRESS=DEFLATE -co TILED=YES

练习 4

创建一个命令，获取 reservoirs.xml 文件并为印度地区创建一个名为 reservoirs_india_reprojected.tif 的 GeoTiff 文件，分辨率为 500m，CRS 为 WGS 84 / India NSF LCC (EPSG:7755)。

2.6 地理配准 (Georeferencing)

GDAL 命令行实用程序在批量地理配准任务中非常有用。我们将学习 2 种不同的地理配准/变形图像技术。

2.6.1 使用边界框坐标进行地理配准

如果已知创建这些图像时使用的边界框坐标和 CRS，请使用 gdal_translate 命令分配地理参考信息。-a_ullr 选项允许你使用左上角 (UL) 和右下角 (LR) 坐标指定边界框坐标。

Windows

gdal_translate -a_ullr -180 90 180 -90 -a_srs EPSG:4326 ^
  earth_at_night.jpg earth_at_night.tif ^
  -co COMPRESS=JPEG -co TILED=YES -co PHOTOMETRIC=RGB

Mac/Linux

gdal_translate -a_ullr -180 90 180 -90 -a_srs EPSG:4326 \
  earth_at_night.jpg earth_at_night.tif \
  -co COMPRESS=JPEG -co TILED=YES -co PHOTOMETRIC=RGB

2.6.2 使用 GCP 进行地理配准

另一种地理配准图像的方法是使用地面控制点 (GCP)。给定一组 GCP，gdalwarp 可以使用各种变换类型对图像进行地理配准。

第一步是使用 gdal_translate 实用程序将这些 GCP 存储在图像元数据中。

Windows

gdal_translate -gcp 418 893 70 15 -gcp 380 2432 70 5 -gcp 3453 2434  90 5 ^
  -gcp 3407 895 90 15 -gcp 2662 911 85 15 ^
  1870_southern-india.jpg india-with-gcp.tif

Mac/Linux

gdal_translate -gcp 418 893 70 15 -gcp 380 2432 70 5 -gcp 3453 2434  90 5 \
  -gcp 3407 895 90 15 -gcp 2662 911 85 15 \
  1870_southern-india.jpg india-with-gcp.tif

现在 GCP 已存储在图像中，我们准备进行地理配准。假设地图的 CRS 是基于 Everest 1830 基准面的地理 CRS，我们选择 EPSG:4042 作为目标 CRS。

让我们尝试 多项式 1 (Polynomial 1) 变换并检查结果。

Windows

gdalwarp -t_srs EPSG:4042 -order 1 -tr 0.005 0.005  ^
  india-with-gcp.tif india-reprojected-polynomial.tif ^
  -co COMPRESS=JPEG -co JPEG_QUALITY=50 -co PHOTOMETRIC=YCBCR

Mac/Linux

gdalwarp -t_srs EPSG:4042 -order 1 -tr 0.005 0.005 \
  india-with-gcp.tif india-reprojected-polynomial.tif \
  -co COMPRESS=JPEG -co JPEG_QUALITY=50 -co PHOTOMETRIC=YCBCR

练习 5

在 india-with-gcp.tif 文件上尝试薄板样条 (Thin-plate-spline) 变换，并将结果保存为 india-reprojected-tps.tif 文件。注意，薄板样条选项在 gdalwarp 命令中通过 -tps 提供。