3. Python遥感数据处理基础：GDAL库安装与使用、栅格数据读写、波段操作、数据可视化

各位同学，今天我们来聊聊遥感数据处理中最基础、也是最绕不开的一个工具——GDAL。说实话，我刚开始接触遥感那会儿，面对一堆.tif文件，整个人是懵的。后来发现，只要搞定了GDAL，栅格数据就像自家后院一样熟悉。嗯，咱们今天就把这块硬骨头啃下来。

3.1 GDAL库：遥感数据的瑞士军刀

GDAL全称是Geospatial Data Abstraction Library，翻译过来就是地理空间数据抽象库。说白了，它就是一套能读写各种栅格和矢量格式的底层工具。我个人习惯把它比作瑞士军刀——功能多、体积小、哪儿都能用。

为什么遥感离不开GDAL？你想想看，遥感数据格式五花八门：GeoTIFF、ENVI、HDF、NetCDF……每个都有自己的存储结构。GDAL把这些格式统一成一个接口，你只需要学会一套API，就能处理几乎所有格式。我在项目中遇到过客户给了一个老旧的ERS格式数据，其他软件都打不开，最后还是GDAL搞定的。

3.2 GDAL安装：别踩这些坑

安装GDAL，说简单也简单，说麻烦也麻烦。我建议直接用conda安装，省心：

conda install -c conda-forge gdal

或者用pip：

pip install gdal

⚠️ 注意：pip安装GDAL经常遇到编译问题。我曾经在Windows上折腾了整整一下午，最后发现是Visual C++ Redistributable版本不对。如果你也遇到类似问题，直接换conda安装，或者去官网下载预编译的whl文件。

验证安装是否成功，打开Python试试：

from osgeo import gdal
print(gdal.__version__)

如果没报错，恭喜你，环境搭好了。

3.3 栅格数据读写：打开遥感世界的大门

读写栅格数据是基本功。我一般用gdal.Open()来打开文件，用gdal.GetDriverByName()来创建新文件。来看个例子：

from osgeo import gdal
import numpy as np

# 打开一个GeoTIFF文件
dataset = gdal.Open('landsat8.tif')
if dataset is None:
    print("文件打开失败，检查路径")
else:
    print(f"文件打开成功，大小：{dataset.RasterXSize} x {dataset.RasterYSize}")

读取数据时，我习惯用ReadAsArray()方法，它直接把栅格数据转成NumPy数组，方便后续处理：

# 读取第一个波段
band1 = dataset.GetRasterBand(1)
data = band1.ReadAsArray()
print(f"数据形状：{data.shape}，数据类型：{data.dtype}")

💡 小技巧：读取大文件时，别一次性全读进内存。我建议用ReadAsArray(xoff, yoff, xsize, ysize)分块读取，或者用ReadAsArray(buf_obj=...)指定输出缓冲区，能省不少内存。

写入数据也不复杂：

# 创建一个新的GeoTIFF文件
driver = gdal.GetDriverByName('GTiff')
out_dataset = driver.Create('output.tif', 
                            xsize=dataset.RasterXSize, 
                            ysize=dataset.RasterYSize, 
                            bands=1, 
                            eType=gdal.GDT_Float32)

# 设置地理参考信息
out_dataset.SetGeoTransform(dataset.GetGeoTransform())
out_dataset.SetProjection(dataset.GetProjection())

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

# 关闭文件
out_dataset.FlushCache()
out_dataset = None

这里有个坑：写完数据后一定要调用FlushCache()，否则数据可能没真正写入磁盘。我曾经因为这个丢过一整天的工作成果，从那以后每次写完都记得flush。

3.4 波段操作：多光谱数据的灵魂

遥感影像通常有多个波段，比如Landsat 8有11个波段。波段操作说白了就是对这些波段进行组合、提取、运算。我经常用波段组合来做假彩色合成，比如：

# 读取近红外、红、绿三个波段
nir = dataset.GetRasterBand(5).ReadAsArray()  # 近红外
red = dataset.GetRasterBand(4).ReadAsArray()  # 红
green = dataset.GetRasterBand(3).ReadAsArray() # 绿

# 计算NDVI（归一化植被指数）
ndvi = (nir.astype(float) - red.astype(float)) / (nir.astype(float) + red.astype(float) + 1e-10)

🔑 核心要点：波段操作时一定要注意数据类型。整数相除会得到整数结果，所以要先转成float。另外，分母加一个小数防止除零错误，这是我在实际项目中踩过的坑。

波段提取也很常用，比如我只想要某个特定波段：

# 提取第4波段（红光波段）
band4 = dataset.GetRasterBand(4)
band4_data = band4.ReadAsArray()

# 获取波段元数据
print(f"波段4最小值：{band4.GetMinimum()}")
print(f"波段4最大值：{band4.GetMaximum()}")
print(f"波段4描述：{band4.GetDescription()}")

3.5 数据可视化：让数据说话

数据可视化是遥感分析中不可或缺的一环。我习惯用matplotlib来展示结果，简单直观：

import matplotlib.pyplot as plt

# 显示单波段灰度图
plt.figure(figsize=(10, 6))
plt.imshow(band4_data, cmap='gray')
plt.colorbar(label='DN值')
plt.title('Landsat 8 第4波段')
plt.show()

多波段彩色合成更常用：

# 真彩色合成（红、绿、蓝波段）
red = dataset.GetRasterBand(4).ReadAsArray()
green = dataset.GetRasterBand(3).ReadAsArray()
blue = dataset.GetRasterBand(2).ReadAsArray()

# 归一化到0-1范围
rgb = np.stack([red, green, blue], axis=-1)
rgb = (rgb - rgb.min()) / (rgb.max() - rgb.min())

plt.figure(figsize=(12, 8))
plt.imshow(rgb)
plt.title('真彩色合成')
plt.axis('off')
plt.show()

💡 可视化小技巧：遥感影像的DN值范围通常很大，直接显示会偏暗。我建议做线性拉伸或直方图均衡化，效果会好很多。比如用np.percentile()去掉两端的异常值：

# 2%线性拉伸
p2, p98 = np.percentile(band4_data, (2, 98))
stretched = np.clip((band4_data - p2) / (p98 - p2), 0, 1)

下面这张图展示了本章的知识体系，帮你理清思路：

Python遥感数据处理基础：GDAL库知识体系 GDAL库 安装与配置 栅格数据读写 波段操作 数据可视化 conda install pip install 版本验证 gdal.Open() ReadAsArray() WriteArray() 波段提取 波段组合 NDVI计算 matplotlib 灰度/彩色显示

好了，以上就是GDAL处理遥感数据的基础内容。从安装到读写，从波段操作到可视化，每一步都是后续课程的基础。我个人建议你动手敲一遍代码，遇到问题很正常，多查查GDAL的官方文档，或者看看Stack Overflow上的讨论。记住，遥感数据处理没有捷径，但GDAL能帮你走得更稳。

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

← 上一章 📖 返回目录 下一章 →