第1章：GDAL库入门与实践

各位同学好，我是你们这门课的老朋友。今天咱们来聊聊遥感数据处理中最基础、也最绕不开的一个工具——GDAL库。

说实话，我刚开始接触遥感那会儿，面对一堆.tif文件，完全不知道从哪下手。后来一位老前辈甩给我一句话：「小子，先把GDAL玩明白，其他的都好说。」嗯，这话我记到现在。

GDAL全称是Geospatial Data Abstraction Library，翻译过来就是「地理空间数据抽象库」。说白了，它就是个万能钥匙，能读写各种栅格和矢量数据格式。你想想看，遥感影像格式五花八门——GeoTIFF、HDF、NetCDF、IMG……要是每种格式都自己写解析器，那得累死。GDAL帮我们统一了接口，你只需要学会一套API，就能处理几十种格式。

1.1 GDAL的核心概念

要理解GDAL，得先搞明白三个核心概念：数据集、波段、仿射变换。我一个个来讲。

1.1.1 数据集（Dataset）

数据集是GDAL中最顶层的对象。你可以把它理解成「一张遥感影像文件」的抽象代表。打开一个.tif文件，返回的就是一个Dataset对象。这个对象里包含了影像的宽度、高度、波段数、投影信息、地理变换参数等等。

我在项目中遇到过一个问题：同事用ENVI导出的.dat文件，用GDAL打开时报错。后来发现是文件路径里有中文，GDAL对中文路径支持不太好。嗯，这里要注意，文件路径尽量用英文，省得给自己挖坑。

1.1.2 波段（Band）

波段就是影像的各个光谱通道。比如Landsat 8影像有11个波段，每个波段对应一个特定的波长范围。GDAL中，波段用索引访问，索引从1开始（不是0哦，这个坑我踩过）。

每个波段都有自己的属性：

• 波段编号：从1开始
• 数据类型：比如Byte、UInt16、Float32等
• 行列数：跟数据集的行列数一致
• NoData值：无效像素的标记值

1.1.3 仿射变换（GeoTransform）

仿射变换是连接「像素坐标」和「地理坐标」的桥梁。它由6个参数组成：

参数	含义
GT(0)	左上角像素的X坐标（经度或东向）
GT(1)	X方向的分辨率（像素宽度）
GT(2)	X方向的旋转系数（通常为0）
GT(3)	左上角像素的Y坐标（纬度或北向）
GT(4)	Y方向的旋转系数（通常为0）
GT(5)	Y方向的分辨率（像素高度，通常为负值）

为什么GT(5)通常是负值？因为影像的Y坐标是从上往下增加的，而地理坐标的Y是从南往北增加的。方向相反，所以分辨率是负的。这个细节我当年琢磨了好一会儿才想明白。

1.2 读取遥感影像元数据

元数据就是「关于数据的数据」。比如影像的拍摄时间、卫星名称、云量、投影信息等等。GDAL读取元数据非常方便：

from osgeo import gdal

# 打开数据集
dataset = gdal.Open('landsat8_b4.tif')

# 读取元数据
metadata = dataset.GetMetadata()
print('元数据：', metadata)

# 读取投影信息
projection = dataset.GetProjection()
print('投影：', projection)

# 读取仿射变换参数
transform = dataset.GetGeoTransform()
print('仿射变换：', transform)

# 读取影像尺寸
width = dataset.RasterXSize
height = dataset.RasterYSize
bands = dataset.RasterCount
print(f'宽度：{width}，高度：{height}，波段数：{bands}')

我曾经在处理一批国产卫星影像时，发现元数据里缺少投影信息。后来手动从文件名中提取了经纬度范围，重新构建了仿射变换参数。所以啊，元数据不是100%可靠的，有时候需要自己验证一下。

1.3 读取与写入波段数据

波段数据是遥感影像的核心——就是那些像素值。GDAL提供了两种读取方式：

1.3.1 读取整个波段

# 获取第1个波段
band = dataset.GetRasterBand(1)

# 读取整个波段数据（返回numpy数组）
data = band.ReadAsArray()
print('数据类型：', data.dtype)
print('数据形状：', data.shape)
print('最小值：', data.min(), '最大值：', data.max())

1.3.2 读取部分区域

如果影像很大（比如几百MB甚至几个GB），一次性读取全部数据会撑爆内存。这时候可以按需读取：

# 读取从(100, 200)开始，宽300像素、高400像素的区域
subset = band.ReadAsArray(100, 200, 300, 400)
print('子区域形状：', subset.shape)

1.3.3 写入波段数据

写入数据需要先创建一个新的数据集，或者使用已有的数据集：

# 创建新的GeoTIFF文件
driver = gdal.GetDriverByName('GTiff')
out_dataset = driver.Create('output.tif', width, height, 1, gdal.GDT_Float32)

# 设置仿射变换和投影
out_dataset.SetGeoTransform(transform)
out_dataset.SetProjection(projection)

# 写入数据
out_band = out_dataset.GetRasterBand(1)
out_band.WriteArray(new_data)

# 设置NoData值
out_band.SetNoDataValue(-9999)

# 关闭数据集（释放资源）
out_dataset = None

我个人习惯在写入完成后，再读回来验证一下数据是否正确。别嫌麻烦，我曾经因为忘记设置NoData值，导致后续的植被指数计算出现了大量异常值，排查了半天才发现问题。

1.4 地理坐标与像素坐标转换

这是遥感数据处理中最常见的操作之一。比如你有一个GPS点（经纬度），想知道它在影像上的像素位置；或者反过来，你知道影像上的某个像素，想知道它的地理坐标。

转换公式其实很简单，就是仿射变换的逆运算：

def pixel_to_geo(col, row, transform):
    """像素坐标转地理坐标"""
    x = transform[0] + col * transform[1] + row * transform[2]
    y = transform[3] + col * transform[4] + row * transform[5]
    return x, y

def geo_to_pixel(x, y, transform):
    """地理坐标转像素坐标"""
    # 解线性方程组
    det = transform[1] * transform[5] - transform[2] * transform[4]
    col = (transform[5] * (x - transform[0]) - transform[2] * (y - transform[3])) / det
    row = (-transform[4] * (x - transform[0]) + transform[1] * (y - transform[3])) / det
    return int(col), int(row)

# 示例：将像素(500, 300)转为地理坐标
lon, lat = pixel_to_geo(500, 300, transform)
print(f'像素(500, 300)对应的地理坐标：({lon:.4f}, {lat:.4f})')

# 示例：将地理坐标(116.4, 39.9)转为像素坐标
col, row = geo_to_pixel(116.4, 39.9, transform)
print(f'地理坐标(116.4, 39.9)对应的像素位置：({col}, {row})')

GDAL也提供了现成的方法，不用自己手写公式：

# 使用GDAL内置方法
from osgeo import osr

# 创建坐标转换对象
source_srs = osr.SpatialReference()
source_srs.ImportFromEPSG(4326)  # WGS84经纬度

target_srs = osr.SpatialReference()
target_srs.ImportFromWkt(dataset.GetProjection())

transform_coord = osr.CoordinateTransformation(source_srs, target_srs)

# 转换坐标
lon, lat = 116.4, 39.9
x, y, z = transform_coord.TransformPoint(lon, lat)
print(f'投影坐标：({x:.2f}, {y:.2f})')

# 再转像素坐标
col = int((x - transform[0]) / transform[1])
row = int((y - transform[3]) / transform[5])
print(f'像素位置：({col}, {row})')

知识体系总览

下面这张图是我自己画的，把本章的核心知识点串起来了。你保存下来，以后忘了就翻出来看看。

好了，这一章的内容就到这里。GDAL的核心概念其实不复杂，但需要多动手练习。我建议你找一张Landsat或Sentinel-2的影像，按照上面的代码一步步操作，把元数据、波段数据、坐标转换都跑一遍。遇到问题很正常，多试几次就熟练了。

记住：遥感数据处理，GDAL是第一关。这一关过了，后面的植被指数计算、生态评估才能顺利展开。

---

公众号:蓝海资料掘金营，微信deep3321