note

GDAL 入门指南

GDAL（Geospatial Data Abstraction Library，地理空间数据抽象库）是处理栅格地理空间数据的开源"瑞士军刀"。它可以读写数百种格式（GeoTIFF、JPEG、PNG、NetCDF、HDF5、ECW、COG 等），支持重投影、波段运算、裁剪、格式转换、金字塔构建等操作，是遥感与 GIS 领域的基石工具。

1. 安装

Python 环境

bash
pip install gdal          # 某些平台可能失败
pip install gdal==$(python -c "from osgeo import gdal; print(gdal.__version__)")

更推荐的方式：使用 conda，会自动装好 C++ 运行时依赖。

bash
conda install -c conda-forge gdal

验证安装：

python
from osgeo import gdal, osr
print(gdal.__version__)   # 例如 3.9.2

命令行工具

安装 GDAL 后还附带一组常用 CLI 工具：

工具用途
gdalinfo查看数据集元信息（尺寸、CRS、波段数等）
gdal_translate格式转换、简单参数调整
gdalwarp重投影 / 重采样
gdalbuildvrt生成虚拟数据集（VRT）
gdaladdo构建概览（overviews/金字塔）
gdaldemDEM 高程图派生（坡度、山体阴影等）
gdal_contour从栅格提取等值线
nearblack去除黑边
gdal_grid散点转栅格

2. 核心概念

  • Dataset：一个打开的栅格文件对象，所有操作的入口。
  • Band：单个波段（如 RGB 中的 R 波段）。通过 ds.GetRasterBand(i) 获取。
  • Block：磁盘上的存储块，ReadBlock/WriteBlock 按块 I/O，适合大文件。
  • Overview：降采样金字塔，加速大范围浏览显示。
  • NoData：无效值标记，通过 GetNoDataValue() 读取。
  • GeoTransform：仿射变换矩阵，把像素坐标映射到地理坐标：
    x_geo = gt[0] + x_pix * gt[1] + y_pix * gt[2]
y_geo = gt[3] + x_pix * gt[4] + y_pix * gt[5]
  • SRS / CRS：空间参考系（坐标系），由 osr.SpatialReference 表示，常见 WKT / EPSG 编码 / PROJ 字符串互转。

3. 最小示例：读写 GeoTIFF

3.1 读取

python
from osgeo import gdal

ds = gdal.Open("dem.tif")
if ds is None:
    raise FileNotFoundError("打不开 dem.tif")

print("尺寸:", ds.RasterXSize, "x", ds.RasterYSize)
print("波段数:", ds.RasterCount)
print("数据类型:", gdal.GetDataTypeName(ds.GetRasterBand(1).DataType))
print("GeoTransform:", ds.GetGeoTransform())
print("投影(WKT):", ds.GetProjection())

band = ds.GetRasterBand(1)
data = band.ReadAsArray()            # NumPy 数组 (y, x)
nodata = band.GetNoDataValue()
stats = band.ComputeStatistics(False)  # (min, max, mean, std)
print("统计:", stats)

3.2 写入

python
import numpy as np
from osgeo import gdal, osr

arr = np.random.randint(0, 255, size=(256, 256), dtype=np.uint8)

driver = gdal.GetDriverByName("GTiff")
out = driver.Create("out.tif", 256, 256, 1, gdal.GDT_Byte)
out.SetGeoTransform((440720, 30, 0, 3751320, 0, -30))  # UTM 31N 示例
srs = osr.SpatialReference()
srs.ImportFromEPSG(32631)
out.SetProjection(srs.ExportToWkt())
out.GetRasterBand(1).WriteArray(arr)
out.FlushCache()
del out

4. 常用操作

4.1 重投影（gdalwarp）

bash
gdalwarp -t_srs EPSG:4326 -r bilinear in_utm.tif out_wgs84.tif

Python：

python
from osgeo import gdal

src = gdal.Open("in.tif")
dst = gdal.GetDriverByName("GTiff").CreateCopy("out.tif", src, copyOnly=True)

# 用 WarpOptions 做真正的重采样
warp = gdal.AutoCreateWarpedVRT(src, src.GetGeoTransform(), src.GetProjection(),
                                None, "EPSG:4326", gdal.GRPC_BILINEAR)
gdal.GetDriverByName("GTiff").CreateCopy("out_wgs84.tif", warp, copyOnly=True)

4.2 裁剪子区域

python
sub = ds.GetRasterBand(1).ReadAsArray(xoff=100, yoff=100, xsize=500, ysize=500)

或整幅 ROI：

bash
gdal_translate -projwin 440720 3751320 441220 3750820 in.tif crop.tif

4.3 格式转换

bash
gdal_translate -of JPEG -co "QUALITY=90" in.tif out.jpg
gdal_translate -ot Float32 -scale 0 255 0 1 in.tif normalized.tif

4.4 拼接 / VRT

bash
gdalbuildvrt mosaic.vrt scene1.tif scene2.tif scene3.tif
gdalwarp -t_srs EPSG:32631 -r near mosaic.vrt mosaic.tif

4.5 构建金字塔（Overviews）

bash
gdaladdo -r average in.tif 2 4 8 16 32

Python：

python
ds = gdal.Open("in.tif", gdal.GA_Update)
for b in range(1, ds.RasterCount + 1):
    ds.BuildOverviews('AVERAGE', [2, 4, 8, 16])

4.6 多波段合成（RGB → 单文件）

python
r = gdal.Open("red.tif").GetRasterBand(1).ReadAsArray()
g = gdal.Open("green.tif").GetRasterBand(1).ReadAsArray()
b = gdal.Open("blue.tif").GetRasterBand(1).ReadAsArray()

rgb = np.stack([r, g, b], axis=0)   # (3, h, w)
driver = gdal.GetDriverByName("GTiff")
out = driver.Create("rgb.tif", rgb.shape[1], rgb.shape[2], 3, gdal.GDT_UInt16)
for i, band_arr in enumerate(rgb, start=1):
    out.GetRasterBand(i).WriteArray(band_arr)
out.SetGeoTransform(gdal.Open("red.tif").GetGeoTransform())
out.SetProjection(gdal.Open("red.tif").GetProjection())
out.FlushCache()

4.7 内存 Dataset（无需落盘）

python
mem_drv = gdal.GetDriverByName("MEM")
mem_ds = mem_drv.Create("", width, height, bands, gdal.GDT_Float32)
mem_ds.GetRasterBand(1).WriteArray(arr)
# 可直接作为其它函数的输入

5. 坐标系统（OSR）

python
from osgeo import osr

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)                # WGS84 经纬度
print(srs.ExportToWkt())
print(srs.ExportToProj4())              # +proj=longlat +datum=WGS84 ...

# 任意两种 SRS 互转
srs2 = osr.SpatialReference()
srs2.ImportFromEPSG(32631)              # UTM 31N

transformer = osr.CoordinateTransformation(srs, srs2)
lon, lat = 120.0, 30.0
x, y, _ = transformer.TransformPoint(lon, lat)
print(x, y)                             # UTM 坐标

GDAL 3.x 起内部使用 OGR 的 CRS 对象，ImportFromEPSG 已能正确处理多数情况；老代码里的 SetWellKnownGeogCS("WGS84") 仍可用但推荐 EPSG 方式。

6. 性能技巧

  1. 分块读取：大文件避免一次性 ReadAsArray()，改用 ReadBlock(blockX, blockY) 或 RasterIO 带窗口读取。
  2. RASTERIO_FLAGS：gdal.GRF_ReadOnly、gdal.GRF_RastNoData 等标志位控制行为。
  3. 驱动选项（Creation Options）：写 GeoTIFF 时常用 -co COMPRESS=LZW -co TILED=YES -co BLOCKSIZE=256。
  4. 缓存：gdal.SetCacheUnlimited() 可放开内存缓存限制（注意机器内存）。
  5. 多线程：gdalwarp 自带 -wo NUM_THREADS=8；Python 侧可配合 concurrent.futures 对分块并行。
  6. Virtual Raster (VRT)：不落地即可组合多个文件，I/O 按需读取。

7. 调试与排错

  • gdalinfo -stats file.tif：快速确认文件能否被识别、是否有 NoData。
  • GDAL_DEBUGON=1：环境变量打印详细日志。
  • CPL_CURL_VERBOSE=1：调试网络数据源（HTTP/S3/WMS）。
  • 报错 No such file or directory 多半是驱动未注册，先 gdal.AllRegister()。
  • Windows 上 DLL 找不到：确保用 conda 安装的 gdal，或将 bin 目录加入 PATH。

8. 进一步学习

本指南基于 GDAL 3.x API 编写。不同版本间个别函数签名可能有差异，请以你安装的版本对应文档为准。

更新于:

note

Copyright © w