GDAL 的命令行功能又多又强，用熟了很香，但参数一多就容易混：到底是 gdalwarp 还是 gdal_translate？-te 和 -projwin 哪个是裁剪？COG、overviews、内部分块、NoData、重采样、对齐……每次去翻官方文档都嫌太长。所以在此记录自己常用、并且验证可用的命令模板。

坐标系转换（重投影）

最常用命令就是 gdalwarp：

# 把 input.tif 转成 EPSG:4326
gdalwarp \
  -t_srs EPSG:4326 \
  -r bilinear \
  -dstnodata -9999 \
  input.tif output_4326.tif
  • -t_srs EPSG:4326：目标坐标系
  • -r bilinear：连续型数据用双线性（影像、高程等）
  • -dstnodata -9999：新栅格的 NoData 填什么值

如果源栅格本身没投影信息，需要顺带声明一下原始坐标系：

gdalwarp \
  -s_srs EPSG:3857 \
  -t_srs EPSG:4326 \
  input_raw.tif output_fixed.tif

生成 COG（Cloud Optimized GeoTIFF）

GDAL 新版有直接的 COG driver，用起来很省心：

# 直接从原始 tif 转 COG
gdal_translate \
  input.tif output_cog.tif \
  -of COG \
  -co "COMPRESS=DEFLATE" \
  -co "BLOCKSIZE=512" \
  -co "RESAMPLING=AVERAGE"
  • -of COG：直接生成符合 COG 规范的 GeoTIFF
  • BLOCKSIZE：内部 tile 边长（256 / 512 视情况而定）
  • COMPRESS=DEFLATE：无损压缩，通用保险

如果想走“手搓”路线，也可以是：

  1. 用 gdal_translate 生成 TILED GeoTIFF
  2. 用 gdaladdo 生成 overviews
  3. 再做 copy / 优化

但一般现在直接 -of COG 更省事。

生成 Overviews（金字塔）

前提：GeoTIFF 已经 TILED=YES。

# 给 big.tif 做金字塔（2、4、8、16 级）
gdaladdo \
  -r average \
  big.tif 2 4 8 16
  • -r average：连续型数据常用
  • 对分类 / 掩膜类数据可以用 -r mode 或 nearest

如果打算配合 Web 服务，常见的是搞到 2 4 8 16 32 之类。

内部切片（inner tile）

本质是 在 GeoTIFF 里启用内部分块，方便按块读取。

# 把大图转成 256 x 256 的内部分块
gdal_translate \
  -of GTiff \
  -co "TILED=YES" \
  -co "BLOCKXSIZE=256" \
  -co "BLOCKYSIZE=256" \
  input.tif tiled_256.tif
  • BLOCKXSIZE / BLOCKYSIZE：控制单块大小（常见 256 / 512）
  • 建议配合压缩一起用，例如 -co "COMPRESS=DEFLATE"

如果是给 Web 地图后端/切片服务用，这一步基本是标配。

修改 NoData 值

只改元数据里的 NoData（不动像素本身）

适合“原本没 NoData，或 NoData 写错”的场景：

# 修改为 -9999，原地编辑
gdal_edit.py -a_nodata -9999 input.tif

或者输出到新文件：

gdal_translate \
  -a_nodata -9999 \
  input.tif output_nodata.tif

把某个像素值“变成 NoData”（需要栅格计算）

比如：当前 0 代表空值，想改成 NoData：

gdal_calc.py \
  -A input.tif \
  --calc="A" \
  --NoDataValue=0 \
  --outfile=output_nodata_0.tif

再配合 gdal_edit.py 把 0 这个值声明为 NoData。

简单栅格计算

gdal_calc.py 的日常用法：一行命令搞定加减乘除。

例 1：单位转换（比如数值 × 0.001）

gdal_calc.py \
  -A input.tif \
  --A_band=1 \
  --calc="A * 0.001" \
  --NoDataValue=-9999 \
  --outfile=output_scaled.tif

例 2：把小于 0 的值都设为 0

gdal_calc.py \
  -A input.tif \
  --calc="(A*(A>=0))" \
  --NoDataValue=-9999 \
  --outfile=output_clipped.tif

例 3：两个栅格相减（B - A）

gdal_calc.py \
  -A ras1.tif -B ras2.tif \
  --calc="B - A" \
  --NoDataValue=-9999 \
  --outfile=delta.tif

基本模式就是：

把 A/B/C 当变量，calc 里写 Numpy 风格表达式。

栅格对齐（分辨率 + 网格）

对齐到固定分辨率 & 网格（-tr + -tap）

# 把 input.tif 重采样到 30m 分辨率，并对齐到 30m 网格
gdalwarp \
  -t_srs EPSG:xxxx \
  -tr 30 30 \
  -tap \
  -r bilinear \
  -srcnodata -9999 \
  -dstnodata -9999 \
  input.tif aligned_30m.tif
  • -tr 30 30：目标分辨率
  • -tap：使像元网格对齐到“分辨率的整数倍”网格（避免半个像元偏移）

对齐到“基准栅格”（典型工作流）

  1. 用 gdalinfo base.tif 看基准图的：

    • 分辨率（Pixel Size）
    • Extent（Corner Coordinates）

  2. 用这些参数喂给 gdalwarp：

gdalwarp \
  -t_srs EPSG:xxxx \
  -te xmin ymin xmax ymax \
  -tr xres yres \
  -r bilinear \
  -srcnodata -9999 \
  -dstnodata -9999 \
  input.tif aligned_to_base.tif

实践中可以写个小脚本，从 gdalinfo 自动获取 -te & -tr。

按给定范围裁剪栅格

这里主要记两种思路：

  1. 按坐标范围（bbox）裁剪：用 gdalwarp 或 gdal_translate 的窗口参数
  2. 按矢量边界裁剪：用 -cutline

用 gdalwarp 按 bbox 裁剪

# 按给定坐标范围裁剪（xmin ymin xmax ymax）
gdalwarp \
  -te xmin ymin xmax ymax \
  -te_srs EPSG:4326 \
  -r bilinear \
  -srcnodata -9999 \
  -dstnodata -9999 \
  input.tif clipped_bbox.tif
  • -te xmin ymin xmax ymax：裁剪范围的外接矩形（边界盒）
  • -te_srs EPSG:4326：这些坐标是在哪个坐标系里（很重要）
  • 其他参数可以顺手带上重采样和 NoData 设置

提示：使用 gdalinfo input.tif 可以先看看源栅格的大致范围，避免手抄坐标填错。

用 gdal_translate 按窗口裁剪（projwin）

# 按左上 / 右下坐标裁剪
gdal_translate \
  -projwin ulx uly lrx lry \
  -projwin_srs EPSG:4326 \
  input.tif clipped_projwin.tif
  • ulx uly：左上角（Upper Left X / Upper Left Y）
  • lrx lry：右下角（Lower Right X / Lower Right Y）
  • 一样可以通过 gdalinfo 获取这些值

简单来说： gdalwarp 更偏“重投影 + 裁剪 + 对齐”， gdal_translate 更偏“在原坐标系里抠一块出来”。

使用矢量边界裁剪栅格

gdalwarp -cutline clip.shp -crop_to_cutline input.tif clipped.tif

从 DEM 生成坡度（slope）与坡向（aspect）

只要 DEM 有了，高程有了，计算坡度/坡向基本就是 gdaldem 的标准操作。

为了保证数值正确，需要特别注意的是：DEM 原始栅格必须是投影坐标

从 DEM 生成坡度栅格

默认输出单位是 度：

# 输出坡度（单位：度）
gdaldem slope \
  dem.tif slope_deg.tif \
  -s 1.0 \
  -of GTiff

  • -s：z-factor，高程单位和 XY 单位不一致时要调整

    • 如果 DEM 是“米”，投影也是“米”（UTM 这类），一般 -s 1.0 就行
    • 如果是地理坐标（经纬度、单位是度），就要给一个换算系数，这里略过，遇到再查

如果更习惯 百分比坡度：

# 输出坡度（单位：百分比）
gdaldem slope \
  dem.tif slope_pct.tif \
  -p \
  -s 1.0 \
  -of GTiff
  • -p：表示以百分比形式输出坡度

从 DEM 生成坡向栅格

gdaldem aspect \
  dem.tif aspect.tif \
  -of GTiff \
  -zero_for_flat \
  -trigonometric
  • 坡向默认单位是度（0–360°）
  • -zero_for_flat：对于完全平坦区域，输出 0（方便后续处理）

  • -trigonometric：方位角按数学中的逆时针方向，从正 x 轴（一般是东）开始

    • 不加的话是“地理方式”（0° 指北，顺时针）

如果只想要某种编码方式，可以看下 gdaldem aspect 的帮助，有一些额外选项（比如设定平坦区域的值）。

坡向的后处理，使用 -1 表示平地

为了不把“北向=0°”和“平地”搞混，习惯上，我会将表示平地的值设为 -1。需要两步实现：

  1. 先正常生成 aspect（平地 = NoData）

不要加 -zero_for_flat，让平地区域保持为 NoData：

gdaldem aspect \
  dem.tif aspect_raw.tif \
  -of GTiff \
  -trigonometric

这时平地通常是 NoData=-9999（可以用 gdalinfo aspect_raw.tif | grep NoData 确认一下）。

  1. 用 gdal_calc 把 “平地 NoData” 改成 -1

假设 NoData 是 -9999，用一个简单的表达式重编码：

gdal_calc.py \
  -A aspect_raw.tif \
  --calc="(A==-9999)*-1 + (A!=-9999)*A" \
  --outfile=aspect_flat_minus1.tif \
  --NoDataValue=-9999

解释一下表达式：

  • A==-9999：平地区域（原来的 NoData） → 变成 1
  • (A==-9999)*-1：平地 → -1
  • (A!=-9999)*A：非平地 → 保持原值

  • 两项相加，得到：

    • 平地像元：-1
    • 其余像元：原始 aspect 值（0–360）

如果不想再保留 NoData（全部都是真实值），最后可以顺手去掉 nodata 标记：

gdal_edit.py -unsetnodata aspect_flat_minus1.tif

这样最终结果就是：平地 = -1，其它像元 = 正常坡向角度。

虚拟栅格 & 拼接

从目录新建虚拟栅格（不真正合并，轻量拼接）

# 把 assets 目录下所有 tif 拼成一个 VRT（虚拟栅格）
gdalbuildvrt mosaic.vrt assets/*.tif

常用场景：先用 VRT 试试看效果，后面再决定要不要导出成实体大图。

虚拟栅格转实体 GeoTIFF（带压缩、内部分块）

gdal_translate \
  -of GTiff \
  -co "COMPRESS=JPEG" \
  -co "TILED=YES" \
  mosaic.vrt mosaic.tif
  • COMPRESS=JPEG：适合影像类（照片型），有损但体积小。
  • TILED=YES：为后续 overviews / COG 之类打基础。

8. 一些顺手的通用参数（记在这里省查文档）

  • -co "BIGTIFF=YES"：超大图（> 4GB）时记得加
  • -co "NUM_THREADS=ALL_CPUS"：能并行就并行一点
  • -r nearest / bilinear / cubic / lanczos：重采样常用选项
  • -multi：gdalwarp 多线程