Python读取遥感影像?Rasterio怎么用?

wowwwai GIS研习社 · 工具流程与项目排障

Python读取遥感影像?Rasterio怎么用? 这是很多 GIS 初学者从桌面软件转向 Python 自动化处理时遇到的第一个问题：影像明明可以在 QGIS 或 ArcGIS Pro 里打开，为什么用 Python 读取时却不知道从哪里拿到波段、坐标系、像元大小和地理范围？本文用 Rasterio 演示一套可复现的遥感影像读取流程，重点解决 GeoTIFF 影像读取、元数据查看、波段读取、窗口裁剪和结果验证这些常见任务。

引言：用 Python 读取遥感影像时先解决什么问题

在 GIS 工作中，遥感影像通常不是简单的图片，而是带有空间参考信息的栅格数据。它可能包含多个波段、投影坐标系、仿射变换参数、无效值、压缩方式和像元分辨率。

如果只是想“看一眼图像”，PIL 或 OpenCV 也能读取部分影像文件。但如果你的目标是做空间分析、计算 NDVI、批量裁剪、重投影或和矢量边界叠加，建议直接使用 Rasterio。Rasterio 是 Python GIS 生态中常用的栅格数据读写库，底层依赖 GDAL，语法比直接调用 GDAL 更适合 Python 用户。

本文默认你要读取的是常见 GeoTIFF 遥感影像，例如 Sentinel、Landsat、无人机正射影像或从 GIS 软件导出的 tif 文件。

背景：为什么遥感影像不能只当普通图片读取

普通图片通常关注宽度、高度和颜色值，而遥感影像还必须关注空间信息。对 GIS 用户来说，读取影像时至少要确认以下内容：

• 坐标系 CRS：影像使用的是 WGS84 经纬度坐标，还是某个投影坐标系。
• 仿射变换 transform：像素行列号如何转换为真实地理坐标。
• 像元大小 resolution：一个像元代表地面上的多大范围。
• 波段数量 count：单波段 DEM、多光谱影像、RGB 影像的读取方式不同。
• NoData 值：哪些像元代表无效区域，不能直接参与统计。
• 数据类型 dtype：uint8、uint16、float32 会影响显示、计算和存储。

这也是很多人第一次用 Python 读取遥感影像时出错的原因：代码虽然读出了数组，但没有理解数组与空间位置之间的关系。Rasterio 的优势就在于，它既能读取像元数组，也能保留 GIS 分析所需的空间元数据。

原理：Rasterio 读取遥感影像的核心对象

Rasterio 读取影像时，通常会先打开一个数据集对象。这个对象类似于 GIS 软件中的一个栅格图层，里面包含影像的波段、范围、坐标系和像元值。

最基础的读取逻辑如下：

import rasterio

image_path = "example.tif"

with rasterio.open(image_path) as src:
    print(src.crs)
    print(src.width, src.height)
    print(src.count)
    band1 = src.read(1)

这里需要注意两点：

• src.read(1) 表示读取第 1 个波段，不是 Python 数组中的第 0 个索引。
• 读取结果 band1 是一个 NumPy 数组，形状通常是 行数 × 列数，也就是 height × width。

Rasterio 中最常用的几个属性如下：

属性	含义	常见用途
src.crs	坐标参考系统	判断是否需要重投影
src.transform	仿射变换参数	像素坐标与地理坐标转换
src.bounds	影像地理范围	检查是否覆盖研究区
src.count	波段数量	判断单波段还是多波段
src.width、src.height	列数和行数	判断影像尺寸
src.nodata	无效值	做统计前剔除无效像元
src.profile	完整栅格配置	写出新影像时复用参数

步骤：Python读取遥感影像 Rasterio怎么用

步骤一：安装 Rasterio

如果你使用的是 Conda 环境，推荐通过 conda-forge 安装 Rasterio，这样 GDAL 依赖更容易匹配：

conda create -n gis-raster python=3.11
conda activate gis-raster
conda install -c conda-forge rasterio numpy matplotlib

如果你使用 pip，也可以尝试：

pip install rasterio numpy matplotlib

在 Windows 环境下，如果 pip 安装失败，通常是 GDAL 依赖或编译环境问题。对 GIS 初学者来说，优先使用 Conda 会更省时间。

步骤二：打开 GeoTIFF 并查看基础信息

下面这段代码适合用于第一次检查遥感影像。它不会修改原始文件，只读取基本元数据。

import rasterio

image_path = "example.tif"

with rasterio.open(image_path) as src:
    print("文件驱动:", src.driver)
    print("宽度:", src.width)
    print("高度:", src.height)
    print("波段数:", src.count)
    print("坐标系:", src.crs)
    print("仿射变换:", src.transform)
    print("地理范围:", src.bounds)
    print("NoData:", src.nodata)
    print("数据类型:", src.dtypes)
    print("像元分辨率:", src.res)

运行后，重点看三个信息：

• 坐标系是否为空：如果 src.crs 是 None，说明影像可能没有正确写入空间参考。
• 范围是否合理：经纬度影像范围一般在 -180 到 180、-90 到 90 附近；投影坐标则可能是米级坐标。
• 波段数是否符合预期：RGB 影像通常是 3 个波段，多光谱影像可能更多。

步骤三：读取单个波段为 NumPy 数组

读取第 1 个波段：

import rasterio

with rasterio.open("example.tif") as src:
    band1 = src.read(1)

print(type(band1))
print(band1.shape)
print(band1.min(), band1.max())

如果影像很大，直接计算 min() 和 max() 可能会受到 NoData 值影响。例如 NoData 为 -9999 时，最小值往往没有实际意义。建议先做掩膜处理。

import numpy as np
import rasterio

with rasterio.open("example.tif") as src:
    band1 = src.read(1)
    nodata = src.nodata

if nodata is not None:
    valid = band1 != nodata
    print("有效像元最小值:", band1[valid].min())
    print("有效像元最大值:", band1[valid].max())
else:
    print("最小值:", band1.min())
    print("最大值:", band1.max())

步骤四：读取多波段影像

如果是 RGB 或多光谱影像，可以一次性读取所有波段：

import rasterio

with rasterio.open("example.tif") as src:
    data = src.read()

print(data.shape)

这里的 data.shape 通常是：

(波段数, 行数, 列数)

这和很多图像库的 (行数, 列数, 通道数) 不一样。用于 Matplotlib 显示时，经常需要转换维度。

import rasterio
import matplotlib.pyplot as plt
import numpy as np

with rasterio.open("rgb.tif") as src:
    rgb = src.read([1, 2, 3])

rgb_show = np.transpose(rgb, (1, 2, 0))

plt.imshow(rgb_show)
plt.title("RGB image")
plt.axis("off")
plt.show()

如果显示结果全黑或全白，可能不是代码错了，而是像元值范围不适合直接显示。例如 uint16 影像的值可能在 0 到 10000 之间，需要先做拉伸。

步骤五：按窗口读取大影像，避免内存爆掉

遥感影像可能非常大，直接 src.read() 会一次性把所有数据读入内存。对于无人机正射影像、全国尺度栅格或高分辨率影像，建议使用窗口读取。

import rasterio
from rasterio.windows import Window

with rasterio.open("large_image.tif") as src:
    window = Window(col_off=0, row_off=0, width=1024, height=1024)
    tile = src.read(1, window=window)

print(tile.shape)

窗口参数含义如下：

• col_off：从第几列开始读取。
• row_off：从第几行开始读取。
• width：读取多少列。
• height：读取多少行。

窗口读取适合切片、分块统计、深度学习样本制作和批量处理大影像。

步骤六：根据地理范围读取影像局部区域

很多 GIS 任务不是按像素位置读取，而是按地理坐标范围读取。例如你知道研究区外接矩形，需要从影像中提取对应区域。可以使用 from_bounds。

import rasterio
from rasterio.windows import from_bounds

with rasterio.open("example.tif") as src:
    minx, miny, maxx, maxy = src.bounds

    center_minx = minx + (maxx - minx) * 0.25
    center_maxx = minx + (maxx - minx) * 0.75
    center_miny = miny + (maxy - miny) * 0.25
    center_maxy = miny + (maxy - miny) * 0.75

    window = from_bounds(
        center_minx,
        center_miny,
        center_maxx,
        center_maxy,
        transform=src.transform
    )

    subset = src.read(1, window=window)

print(subset.shape)

这里必须确保输入的范围坐标和影像坐标系一致。如果影像是投影坐标系，不能直接传入经纬度范围；如果研究区边界是 WGS84，经常需要先重投影到影像 CRS。

步骤七：读取影像并写出新的 GeoTIFF

实际项目中，经常需要读取一个波段，做计算后写出新栅格。写出时不要只保存数组，还要保留坐标系、transform 和 NoData 等信息。

import rasterio
import numpy as np

input_path = "example.tif"
output_path = "band1_copy.tif"

with rasterio.open(input_path) as src:
    band1 = src.read(1)
    profile = src.profile.copy()

    profile.update(
        count=1,
        dtype=band1.dtype,
        compress="lzw"
    )

    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(band1, 1)

这段代码的关键是 profile。它包含了写出 GeoTIFF 所需的大部分空间参数。如果不复用或正确更新 profile，很容易生成“能打开但位置不对”的影像。

常见坑：Rasterio读取遥感影像容易出错的地方

坑一：波段编号从 1 开始，不是从 0 开始

NumPy 数组索引从 0 开始，但 Rasterio 的波段编号从 1 开始。读取第一个波段应写：

src.read(1)

不要写成：

src.read(0)

坑二：数组维度和图像显示习惯不同

Rasterio 多波段读取结果是 (bands, rows, cols)，而 Matplotlib 显示 RGB 通常需要 (rows, cols, bands)。如果颜色错乱或报维度错误，先检查维度顺序。

坑三：NoData 没有处理，统计结果不可信

很多影像边缘或空白区域会使用 NoData 值。直接对整幅影像求均值、最小值、最大值，可能会把无效像元算进去。更稳妥的写法是读取掩膜数组：

import rasterio

with rasterio.open("example.tif") as src:
    band1_masked = src.read(1, masked=True)

print(band1_masked.mean())
print(band1_masked.min())
print(band1_masked.max())

masked=True 会返回 NumPy 的 MaskedArray，有助于自动忽略无效像元。

坑四：坐标系为空或坐标系不一致

如果 src.crs 是 None，说明影像没有可用 CRS。此时不能可靠地做坐标转换、叠加分析或按地理范围裁剪。

如果矢量边界和影像 CRS 不一致，不能直接用边界坐标去裁剪影像。应先把矢量数据重投影到影像坐标系，或者把范围坐标转换到影像 CRS。

坑五：直接读取超大影像导致内存不足

如果影像尺寸达到几万行几万列，一次性读取所有波段会占用大量内存。建议优先使用窗口读取、分块遍历或生成影像金字塔后再用于浏览显示。

方法比较：Rasterio、GDAL、OpenCV、PIL 该怎么选

工具	适合场景	优点	局限
Rasterio	Python GIS 栅格读写、GeoTIFF 处理、窗口读取	语法清晰，保留 CRS、transform、bounds 等空间信息	复杂格式转换和底层高级功能仍可能需要 GDAL
GDAL	专业栅格格式转换、重投影、批处理	功能强，格式支持广	Python API 对初学者不够直观
OpenCV	图像识别、计算机视觉、深度学习预处理	图像处理能力强，速度快	通常不关注遥感影像的空间参考信息
PIL	普通图片读取和简单格式转换	轻量，容易上手	不适合严肃 GIS 栅格分析
QGIS	可视化检查、手动裁剪、交互式分析	直观，适合验证结果	批量自动化能力不如 Python 脚本灵活

如果你的任务是 Python读取遥感影像 并继续做 GIS 分析，Rasterio 通常是首选。如果只是把影像当普通图片输入模型，可以结合 OpenCV，但仍建议用 Rasterio 先确认空间信息和数据范围。

检查清单：读取遥感影像后如何确认结果正确

每次用 Rasterio 读取遥感影像后，建议按下面清单检查一遍：

• 是否能成功打开文件，没有路径、权限或格式错误。
• src.crs 是否存在，坐标系是否符合项目要求。
• src.bounds 是否落在合理的地理范围内。
• src.width、src.height 是否和 GIS 软件中看到的尺寸一致。
• src.count 是否和预期波段数一致。
• src.nodata 是否已识别，统计前是否排除了 NoData。
• 多波段影像显示时，是否正确处理了维度顺序。
• 大影像是否采用窗口读取，避免一次性占用过多内存。
• 写出新 GeoTIFF 时，是否保留了 CRS、transform、dtype 和 nodata。
• 输出结果是否能在 QGIS 或 ArcGIS Pro 中正确叠加到原始数据位置。

FAQ：Python读取遥感影像常见问题

Q1：Rasterio 能读取哪些遥感影像格式？

Rasterio 基于 GDAL，常见 GeoTIFF、部分 IMG、VRT 等栅格格式通常都可以读取。但具体支持情况取决于本机 GDAL 编译时启用的驱动。实际项目中，GeoTIFF 是最推荐的交换格式。

Q2：为什么 Rasterio 读取的数组和 QGIS 显示效果不一样？

QGIS 显示影像时通常会自动做拉伸、设置色带或忽略 NoData，而 Rasterio 读取的是原始像元值。如果直接用 Matplotlib 显示 uint16 遥感影像，可能会偏黑或偏白。需要根据数据范围做归一化或百分比拉伸。

Q3：Python读取遥感影像时为什么 src.crs 是 None？

这说明影像文件中没有正确写入坐标参考系统，或者该格式的空间参考没有被当前环境识别。可以先在 QGIS 中查看图层属性，确认是否存在 CRS。如果确实缺失，需要从数据来源确认正确坐标系后再赋予，而不是随便猜一个。

Q4：Rasterio 怎么读取指定范围内的影像？

如果你有的是像素行列号，可以使用 Window。如果你有的是地理坐标范围，可以使用 from_bounds 生成窗口。前提是范围坐标必须和影像 CRS 一致。

Q5：读取多波段影像时，为什么 shape 是 bands、rows、cols？

Rasterio 遵循栅格数据集的波段组织方式，所以返回结果是 (波段数, 行数, 列数)。如果要用于普通 RGB 显示，通常需要用 np.transpose 转成 (行数, 列数, 波段数)。

Q6：Rasterio 适合做 NDVI 计算吗？

适合。NDVI 本质上是读取红光波段和近红外波段后进行数组计算。但要注意波段顺序、NoData、数据类型和除零问题。计算完成后写出结果时，需要保留原影像的空间参考信息。

import rasterio
import numpy as np

with rasterio.open("multiband.tif") as src:
    red = src.read(3).astype("float32")
    nir = src.read(4).astype("float32")
    profile = src.profile.copy()

ndvi = (nir - red) / (nir + red + 1e-6)

profile.update(count=1, dtype="float32")

with rasterio.open("ndvi.tif", "w", **profile) as dst:
    dst.write(ndvi, 1)

结论：Rasterio怎么用，关键是同时读懂数组和空间信息

学习 Python读取遥感影像 时，不要只停留在把 tif 读成 NumPy 数组。对 GIS 工作来说，更重要的是理解 Rasterio 返回的 CRS、transform、bounds、NoData 和 profile。它们决定了影像能否正确叠加、裁剪、统计和写出。

如果你刚开始使用 Rasterio，建议先掌握四个动作：打开影像查看元数据、读取单波段、读取多波段、按窗口读取大影像。再进一步学习掩膜裁剪、重投影、分块处理和指数计算。这样就能把遥感影像从“能打开”推进到“能自动化分析”。