【Python & RS】Rasterio 库安装+函数使用教程

Rasterio 是一个专门用于处理栅格数据（如遥感影像、地理信息系统中的栅格数据等）的 Python 库。它可以帮助我们读取、写入和操作地理空间数据，尤其是在遥感影像处理、地理数据分析等领域具有广泛的应用。

本文将通过详细的步骤，帮助你安装和使用 Rasterio 库，包括代码示例和图解，助你更好地学习栅格数据处理。

一、安装 Rasterio 库

在使用 Rasterio 库之前，首先需要安装它。你可以使用 pip 命令进行安装：

pip install rasterio

由于 Rasterio 依赖于 GDAL（Geospatial Data Abstraction Library），在某些平台上可能需要手动安装 GDAL。如果遇到问题，可以参考 Rasterio 安装文档 获取更多安装指导。

二、Rasterio 基本功能概述

Rasterio 提供了一系列工具，可以让我们读取、写入、处理栅格数据，并支持地理信息系统（GIS）中的常见数据格式，如 GeoTIFF。

核心功能：

• 读取栅格数据：支持多种栅格数据格式，如 GeoTIFF、HDF5、NetCDF 等。
• 栅格元数据处理：获取栅格图像的基本信息，如坐标参考系统（CRS）、地理坐标等。
• 图像切片与操作：可以对栅格数据进行子区域提取和数据变换。
• 写入栅格数据：将处理后的数据保存为不同的栅格文件格式。

三、Rasterio 常用函数和示例

1. 读取栅格数据

读取栅格数据的最常见方式是使用 rasterio.open() 打开文件，然后通过 .read() 方法读取图像数据。

示例代码

import rasterio

# 打开 GeoTIFF 文件
with rasterio.open('example.tif') as src:
    # 读取所有波段的栅格数据
    data = src.read()

    # 获取栅格的元数据
    print("CRS:", src.crs)
    print("Width, Height:", src.width, src.height)
    print("Bounds:", src.bounds)

说明：

• rasterio.open()：打开栅格文件。
• .read()：读取栅格数据，返回一个 numpy 数组，其中每个波段的数据都在不同的维度中。
• .crs：返回栅格的坐标参考系统（Coordinate Reference System）。
• .bounds：获取栅格数据的地理边界（即左下角和右上角的坐标）。

2. 读取指定波段的栅格数据

如果你的栅格数据包含多个波段（例如 RGB 图像），可以使用 .read(band_number) 读取特定波段的数据。

示例代码

with rasterio.open('example.tif') as src:
    # 读取第一波段的数据
    band1 = src.read(1)
    
    # 输出波段的最小值和最大值
    print(f"波段1 - Min: {band1.min()}, Max: {band1.max()}")

说明：

• src.read(1)：读取第一个波段的数据，返回一个二维的 numpy 数组。

3. 栅格数据的坐标变换

Rasterio 支持坐标系统的转换。例如，如果你需要将栅格数据从一个坐标参考系统（CRS）转换到另一个，可以使用 rasterio.warp 模块。

示例代码

from rasterio.warp import calculate_default_transform, reproject, Resampling

with rasterio.open('example.tif') as src:
    # 获取目标 CRS，假设目标是 EPSG:4326 (WGS 84)
    dst_crs = 'EPSG:4326'

    # 计算转换矩阵
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)

    # 创建目标栅格数据
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    # 执行栅格重投影
    with rasterio.open('reprojected.tif', 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

说明：

• calculate_default_transform()：计算从源 CRS 到目标 CRS 的变换。
• reproject()：执行栅格数据的重投影，改变栅格数据的坐标系统。

4. 写入栅格数据

使用 Rasterio 可以将处理后的栅格数据保存到新的文件中。

示例代码

import numpy as np

# 创建一个简单的数组作为栅格数据
data = np.random.random((100, 100))

# 设置栅格的元数据
kwargs = {
    'driver': 'GTiff',
    'count': 1,  # 波段数量
    'dtype': 'float32',
    'crs': 'EPSG:4326',
    'transform': rasterio.transform.from_origin(-180, 90, 1, 1),  # 假设栅格的左上角坐标为 (-180, 90)
    'width': 100,
    'height': 100
}

# 写入文件
with rasterio.open('output.tif', 'w', **kwargs) as dst:
    dst.write(data, 1)

说明：

• driver='GTiff'：指定输出文件格式为 GeoTIFF。
• .write(data, 1)：将数据写入第一个波段。

5. 栅格数据的掩膜（Mask）

有时我们只关心栅格中的某一部分数据，可以通过掩膜来获取特定区域的数据。

示例代码

from rasterio.mask import mask
import geojson

# 读取一个 GeoJSON 文件作为掩膜
with open('polygon.geojson') as f:
    geojson_data = geojson.load(f)

with rasterio.open('example.tif') as src:
    # 使用 GeoJSON 文件的几何来创建掩膜
    out_image, out_transform = mask(src, geojson_data['features'], crop=True)

    # 输出掩膜区域数据
    print(out_image)

说明：

• mask()：根据给定的几何掩膜提取栅格数据。
• geojson_data['features']：GeoJSON 数据中的多边形区域，作为掩膜区域。

四、常见问题与优化

1. 如何处理大文件？

对于大文件，可以使用 rasterio 的 逐块读取 或 内存映射 功能，避免内存溢出。使用 .read() 时，指定块读取（如 window）可以有效减少内存消耗。

2. 写入时的坐标系统不同怎么办？

确保写入栅格时的 CRS 和原始数据的 CRS 一致。如果需要转换，可以先进行 CRS 转换，然后再进行保存。

五、总结

Rasterio 是一个功能强大且易于使用的栅格数据处理库，适合处理遥感影像、地理数据分析和栅格图像的读写工作。通过本文的学习，你已经掌握了 Rasterio 的基本使用方法，包括栅格数据的读取、处理、写入和坐标变换等。

掌握 Rasterio 后，你可以轻松处理各种地理空间数据，支持进一步的遥感分析和 GIS 应用。