数据预处理与读取：GDAL库安装、读取GeoTIFF文件、获取影像元数据、波段信息提取

好了，咱们正式开始动手干活。

前面聊了那么多遥感数据的理论，什么分辨率、坐标系、存储格式……但说一千道一万，最终还是要落到代码上。你想想看，数据拿不到手里，一切都是空谈。

这一章，我们就来解决最核心的问题：怎么用Python把遥感数据读进来。我会带你从GDAL库的安装开始，一步步走到能提取出影像的元数据和波段信息。嗯，这里要注意，很多新手在这一步就卡住了，不是装不上库，就是读出来的数据是乱的。别担心，我把我在项目中踩过的坑都给你指出来。

3.1 GDAL库：遥感界的“瑞士军刀”

说到遥感数据的读取，GDAL（Geospatial Data Abstraction Library）是绕不开的。我个人习惯叫它“盖兜儿”，圈里人都这么叫。它几乎支持所有主流的栅格和矢量格式，GeoTIFF、IMG、HDF、NetCDF……你叫得上名字的，它基本都能处理。

说白了，它就是一座桥梁。把你的遥感数据和Python世界连接起来。

3.4.1 安装GDAL：一个“劝退”很多人的步骤

安装GDAL，可能是整个课程里最让人头疼的一步。我曾经在给一个新同事配环境时，折腾了整整一个下午。为什么会这样？因为GDAL底层是C++写的，它依赖很多本地的地理空间库。

这里我推荐两种最稳妥的方式：

⚠️ 避坑指南： 千万不要直接用 pip install gdal！这通常会失败，或者装上一个残缺版本。我曾经这么干过，结果编译报错，浪费了两个小时。

方法一：使用conda（强烈推荐）

如果你装了Anaconda或者Miniconda，这是最省心的方式。conda会自动帮你处理所有底层依赖。

conda install -c conda-forge gdal

等它跑完，基本就稳了。我所有项目都是用conda来管理GDAL的。

方法二：使用预编译的wheel文件

如果你坚持用pip，可以去 Christoph Gohlke的页面 下载对应你Python版本和系统的whl文件。然后：

pip install GDAL-3.4.3-cp39-cp39-win_amd64.whl

注意文件名里的 cp39 代表Python 3.9，别下错了。

💡 小技巧： 装完后，在Python里输入 from osgeo import gdal，如果不报错，就说明装好了。我每次装完都会先跑这一句，图个心安。

3.2 读取GeoTIFF：打开遥感世界的大门

GeoTIFF是最常见的遥感数据格式。它就是在普通TIFF图像的基础上，嵌入了地理坐标信息。所以，我们不仅能看图像，还能知道每个像素对应地球上的哪个位置。

读取一个GeoTIFF文件，核心就三步：

  1. 打开文件：用GDAL的 Open() 函数。
  2. 获取数据集：得到一个 Dataset 对象。
  3. 读取栅格波段：用 GetRasterBand() 方法。

来看代码：

from osgeo import gdal
import numpy as np

# 1. 打开文件
file_path = "path/to/your/landsat_image.tif"
dataset = gdal.Open(file_path)

# 2. 检查是否打开成功
if dataset is None:
    print("文件打开失败，请检查路径")
    exit()

print("文件打开成功！")

# 3. 读取第一个波段（索引从1开始）
band1 = dataset.GetRasterBand(1)
data_array = band1.ReadAsArray()

print(f"波段1的数据形状: {data_array.shape}")
print(f"数据类型: {data_array.dtype}")

你看，就这么几行，一个遥感影像的像素数据就已经躺在NumPy数组里了。之后你想做啥分析都行。

🔑 关键点： GDAL的波段索引是从 1 开始的，不是0！这是C语言的老传统了。我第一次用的时候，习惯性地写了 GetRasterBand(0)，结果返回了None，还以为是文件坏了。

3.3 获取影像元数据：给数据“验明正身”

光有像素数据还不够。你得知道这张图是什么时候拍的？坐标系是什么？分辨率是多少？这些信息都藏在元数据里。

我习惯把元数据比作影像的“身份证”。拿到它，你才能正确地使用这张图。

3.3.1 核心元数据提取

下面这段代码，几乎是我每个遥感项目的“开场白”：

# 获取影像的列数、行数、波段数
cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount

print(f"影像大小: {cols} x {rows} 像素")
print(f"波段数量: {bands}")

# 获取地理变换参数（仿射变换矩阵）
geotransform = dataset.GetGeoTransform()
print(f"左上角坐标 (X, Y): ({geotransform[0]:.2f}, {geotransform[3]:.2f})")
print(f"像素宽度 (X方向分辨率): {geotransform[1]:.2f}")
print(f"像素高度 (Y方向分辨率): {geotransform[5]:.2f}")

# 获取投影信息（WKT格式）
projection = dataset.GetProjection()
print(f"投影坐标系: {projection[:50]}...")  # 只打印前50个字符

这里有个细节，GetGeoTransform() 返回的是一个包含6个元素的元组。它描述的是图像坐标和地理坐标之间的转换关系。你不需要背下来，但要知道它能告诉你：这张图的左上角在哪里，每个像素代表地面多大面积。

💡 我的经验： 拿到数据后，我第一件事就是打印 RasterXSize 和 RasterYSize。如果一张Landsat影像显示只有几百像素，那大概率是只读了缩略图，或者文件没读全。

3.4 波段信息提取：深入数据的“灵魂”

多光谱影像有多个波段。每个波段记录了不同波长的电磁波信息。比如，近红外波段对植被特别敏感，热红外波段能反映地表温度。

提取波段信息，就是要把每个波段的“个性”摸清楚。

3.4.1 遍历所有波段

我通常会写一个循环，把每个波段的基本信息都打印出来：

for i in range(1, bands + 1):
    band = dataset.GetRasterBand(i)
    
    # 波段元数据
    band_description = band.GetDescription()
    band_no_data = band.GetNoDataValue()
    band_min = band.GetMinimum()
    band_max = band.GetMaximum()
    
    # 如果没计算过统计值，可以手动计算
    if band_min is None or band_max is None:
        band_min, band_max = band.ComputeRasterMinMax()
    
    print(f"波段 {i}:")
    print(f"  描述: {band_description if band_description else '无'}")
    print(f"  数据类型: {gdal.GetDataTypeName(band.DataType)}")
    print(f"  无效值: {band_no_data}")
    print(f"  像素值范围: {band_min:.2f} - {band_max:.2f}")
    print("-" * 30)

这段代码会输出类似这样的信息：

波段 描述 数据类型 像素值范围
1 Coastal aerosol UInt16 0 - 65535
2 Blue UInt16 0 - 60000
3 Green UInt16 0 - 58000
4 Red UInt16 0 - 62000
5 NIR UInt16 0 - 65000
⚠️ 注意： 很多Landsat数据的无效值（NoData）是 -9999 或者 0。如果你不处理这些值，直接做统计分析，结果会完全跑偏。我见过有人算NDVI，结果出来一堆负数，就是因为没屏蔽无效值。

3.5 本章知识体系：一张图看懂

说了这么多，我们来梳理一下。下面这张SVG图，展示了数据预处理与读取的核心逻辑：

数据预处理与读取：核心流程 1. 安装GDAL库 conda / pip / whl 2. 打开GeoTIFF gdal.Open() 3. Dataset对象 影像的“句柄” A. 提取元数据 • 行列数：RasterXSize/YSize • 地理变换：GetGeoTransform() • 投影信息：GetProjection() B. 提取波段信息 • 遍历波段：GetRasterBand() • 数据类型、无效值 • 像素值范围 C. 读取像素数据 ReadAsArray() → NumPy数组

这张图把整个流程串起来了。从安装GDAL开始，到打开文件，再到三个分支：提取元数据、提取波段信息、读取像素数据。你跟着这个流程走，基本不会迷路。

3.6 写在最后

数据读取是遥感分析的基石。这一步走得稳，后面的处理和分析才能顺利。我个人觉得，花点时间把GDAL的安装和基本API搞透，是非常值得的。

嗯，这一章的内容就到这里。代码不多，但都是干货。你可以在自己的电脑上找一张GeoTIFF试试看，把元数据和波段信息都打印出来。相信我，当你看到那些数字在屏幕上滚动时，会很有成就感。

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

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