准备工作与环境搭建
安装 Python 与虚拟环境
对于 GIS 从业者而言，建立一个干净且可重复的开发环境是第一步。使用虚拟环境可以避免库版本冲突，确保后续的 Rasterio 与依赖可控稳定。
推荐采用 conda 或 venv 来创建一个专门的工作环境，并将 Python 版本锁定在 3.9～3.11 之间，以与主流遥感工具栈保持兼容。命名为 gis_env 以便于辨识，再统一安装核心包。
安装 Rasterio 与依赖
Rasterio 是处理栅格数据的核心库，它对 GDAL 的绑定较深，因此在安装时需要确保底层依赖就绪。优先使用 Conda 安装可避免 C 依赖问题，同时可通过 Pip 下载 Python 包。
基础安装示例见下方，确保在激活的环境中执行。Rasterio、Numpy 与 Matplotlib 是常用组合，便于数据处理与绘图。
# 使用 conda 创建并激活环境
conda create -n gis_env python=3.11
conda activate gis_env# 安装 rasterio、numpy、matplotlib 等
conda install -c conda-forge rasterio numpy matplotlib
Rasterio 的核心概念与卫星影像特点
Rasterio 的核心概念
数据集（dataset）是 Rasterio 读取栅格影像的基础对象，封装了影像的元数据、波段数据及投影信息。波段（band）是栅格数据的最小读取单元，常用于组合、计算和分析。
仿射变换（transform）描述像素网格与地理坐标系的映射关系，结合 CRS（坐标参考系），我们可以实现从像素坐标到地理坐标的转换。
卫星影像数据的波段、分辨率与坐标系
卫星影像通常包含多个波段，每个波段对应不同的光谱通道。选择合适的波段组合（如红光与近红外）可实现植被指数等分析，同时要关注分辨率与投影的一致性。
在读取数据时，影像的坐标系（CRS）与仿射变换是分析的关键要素，它们决定了后续裁切、重采样和对齐的正确性。
import rasterio# 打开影像并查看基本元数据
path = 'data/LC08_L1TP_023030_20200808_20200820_B4.TIF'
with rasterio.open(path) as src:print(src.name)print('波段数:', src.count)print('分辨率:', src.res)print('CRS:', src.crs)print('变换:', src.transform)
读取卫星影像元数据与波段提取
打开影像获取元数据
通过 rasterio.open 可以快速加载影像并读取元数据，例如波段数、图像尺寸、像元大小以及 nodata 值等信息。元数据是后续波段操作与数据对齐的基础。
掌握这些信息有助于评估影像质量、选择波段以及确定后续的处理策略。与 GIS 工作流高度相关的要点包括投影、变换和 nodata 设置。
读取单波段与多波段数据
影像中不同波段包含不同信息，单波段读取用于简单统计与可视化，多波段读取适合进行指数计算等高级分析。
在实际应用中，先读出需要的波段再进行数组级运算，确保不要将 nodata 值混入计算，以提升结果的准确性。
import numpy as np
import rasteriopath = 'data/LC08_L1TP_023030_20200808_20200820_B4.TIF' # 红波段示例
with rasterio.open(path) as src:red = src.read(1) # 读取第一波段（示例）nodata = src.nodataprint('红波段形状:', red.shape, 'nodata:', nodata)
数据处理：波段运算、掩膜与裁切
波段组合与简单指数计算
将不同波段组合进行运算是遥感分析的核心，例如常见的 NDVI、NDWI 等。使用 NumPy 数组进行逐元素计算，并注意处理 NA 值和除零情况。
计算结果通常为单通道栅格，保持原始数据类型或转换为浮点型以获得精度，再决定是否添加 nodata 掩膜。
import numpy as np# 假设 red 与 nir 已经从不同波段读取为 2D 数组
# nir = 近红外波段，red = 红波段
ndvi = (nir.astype(float) - red.astype(float)) / (nir + red)
ndvi = np.where((nir + red) == 0, -9999, ndvi) # 使用 -9999 作为 nodata
掩膜、裁切与重采样
在处理卫星影像时，常需要对感兴趣区域进行裁切，或对不同分辨率影像进行重采样以实现对齐。掩膜操作可将无效像元排除在分析之外，提升结果的可用性。
裁切与重采样的实现要点在于正确传递 变换（transform）与 CRS，以及选择合适的重采样算法（最近邻、双线性等）。
from rasterio.windows import Window# 假设 ndvi 已经准备好，进行简单裁切示例
with rasterio.open('data/NDVI.tif') as src:# 以左上角坐标和尺寸裁切一个子区域win = Window(100, 100, 500, 500)ndvi_crop = src.read(1, window=win)transform = src.window_transform(win)
数据输出与可视化
写出新栅格数据（GeoTIFF）
分析结果需要以栅格格式输出，常用 GeoTIFF 保存，保持 CRS、变换信息以及 nodata 设置的一致性，以便在 GIS 软件中进一步分析。
在写出新栅格时，要更新数据类型、波段数与 nodata，确保结果能正确被后续工具读取。
import rasterio# 使用与输入影像相同的 profile，更新为单波段浮点输出
out_path = 'output_ndvi.tif'
with rasterio.open('data/LC08_ndvi_input.tif') as src:profile = src.profile
profile.update(dtype=rasterio.float32, count=1, nodata=-9999)with rasterio.open(out_path, 'w', **profile) as dst:dst.write(ndvi.astype(rasterio.float32), 1)
快速可视化结果
将分析结果进行简单可视化有助于快速判断分析是否合理。Matplotlib 提供直观的色带映射，适用于展示 NDVI、NDWI 等指数，也便于在博客或报告中呈现。
在可视化时，选择合适的颜色映射和色标范围，避免误导读者对地表信息的解读。
import matplotlib.pyplot as pltplt.figure(figsize=(6, 6))
plt.imshow(ndvi, cmap='RdYlGn', vmin=-1, vmax=1)
plt.colorbar(label='NDVI')
plt.title('NDVI 可视化')
plt.axis('off')
plt.show()
