第九章：栅格分析步骤详解

栅格分析步骤（raster.*）处理栅格影像数据，共 5 个步骤，覆盖投影转换、裁剪、波段计算、统计和等值线提取。所有栅格步骤依赖 rasterio 库，基础安装已包含。

9.1 步骤总览

步骤 ID	名称	功能
`raster.reproject`	栅格投影转换	栅格数据重投影（rasterio warp）
`raster.clip`	栅格裁剪	用矢量掩膜裁剪栅格
`raster.calc`	栅格计算	波段数学运算（如 NDVI）
`raster.stats`	栅格统计	波段统计信息（min/max/mean/std）
`raster.contour`	等值线提取	从栅格生成等值线矢量数据

9.2 栅格数据结构说明

栅格步骤使用字典格式传递数据（raster_info），包含以下键：

键	类型	说明
`data`	ndarray	像素数组，shape `[bands, height, width]`
`transform`	Affine	仿射变换矩阵（geotransform）
`crs`	CRS	坐标参考系对象
`profile`	dict	rasterio 完整元信息

从 io.read_raster 读取后，这个字典通过 $step.output（或简写 $step）在步骤间传递；io.read_raster 的 stats 中另有 width、height、crs、band_count、dtype、bounds 等摘要信息（nodata 等完整元信息保存在 profile 内）。

9.3 `raster.reproject`：栅格投影转换

参数

参数	类型	必填	默认值	说明
`input`	raster_info	✅	—	输入栅格数据
`target_crs`	string	✅	—	目标坐标系（EPSG 代码）
`resampling`	string	❌	`nearest`	重采样方法：`nearest`/`bilinear`/`cubic`/`lanczos`

重采样方法选择

数据类型	推荐方法
分类数据（地类、土地利用）	`nearest`（最近邻，不改变像素值）
连续数值（DEM、温度）	`bilinear`（双线性，平滑过渡）
高质量连续数据	`cubic`（三次样条）
极高质量图像	`lanczos`

示例

- id: load-dem
  use: io.read_raster
  params: { path: "data/dem_wgs84.tif" }

- id: reproject-dem
  use: raster.reproject
  params:
    input: "$load-dem"
    target_crs: "EPSG:3857"
    resampling: "bilinear"

- id: save-reprojected
  use: io.write_raster
  params:
    input: "$reproject-dem"
    path: "output/dem_3857.tif"

9.4 `raster.clip`：栅格裁剪

参数

参数	类型	必填	默认值	说明
`input`	raster_info	✅	—	输入栅格数据
`mask`	geodataframe	✅	—	裁剪掩膜（矢量数据，需与栅格坐标系一致）
`crop`	boolean	❌	`true`	是否裁剪到掩膜范围
`nodata`	number	❌	`0`	NoData 填充值

示例

- id: load-satellite
  use: io.read_raster
  params: { path: "data/landsat_scene.tif" }

- id: load-boundary
  use: io.read_vector
  params: { path: "data/study_area.shp" }

# 确保矢量与栅格坐标系一致
- id: reproject-boundary
  use: vector.reproject
  params:
    input: "$load-boundary"
    target_crs: "EPSG:32650"    # 与 landsat_scene.tif 坐标系一致

- id: clip-satellite
  use: raster.clip
  params:
    input: "$load-satellite"
    mask: "$reproject-boundary"
    crop: true
    nodata: -9999

- id: save-clipped
  use: io.write_raster
  params:
    input: "$clip-satellite"
    path: "output/clipped_scene.tif"

9.5 `raster.calc`：栅格计算

raster.calc 是一个功能强大的波段计算步骤，使用 B1、B2… 引用各波段，支持 NumPy 数学运算。

参数

参数	类型	必填	默认值	说明
`input`	raster_info	✅	—	输入栅格数据
`expression`	string	✅	—	计算表达式，使用 B1, B2… 引用波段

安全性

表达式通过 utils/safe_eval.py 中的 AST 白名单求值：仅允许波段引用、数值运算（+、-、*、/、//、%、**）、比较与布尔运算，以及 np. 前缀的白名单数学函数（abs、sqrt、log、exp、where、clip、minimum、maximum、nan_to_num 等三角/对数/极值函数），其余函数调用一律拒绝，防止代码注入。

常用表达式示例

# NDVI (植被指数) = (近红外 - 红) / (近红外 + 红)
# Landsat 8: B5=近红外, B4=红
expression: "(B5 - B4) / (B5 + B4)"

# EVI (增强植被指数)
expression: "2.5 * (B5 - B4) / (B5 + 6 * B4 - 7.5 * B2 + 1)"

# 波段相加
expression: "B1 + B2 + B3"

# 归一化（np.* 白名单函数可用）
expression: "np.clip(B1 / 255.0, 0, 1)"

注意：表达式中只能使用 B1、B2… 形式引用波段，不支持 band1、Band1 等写法。B1.min() 之类的方法调用不在白名单内；只允许基本数学运算符、括号与 np.* 白名单函数。

NDVI 完整示例

pipeline:
  name: "NDVI 计算"
  description: "计算研究区 NDVI，输出为 GeoTIFF"

  variables:
    input_path: "data/landsat8_multispectral.tif"

  steps:
    - id: load-multispectral
      use: io.read_raster
      params: { path: "${input_path}" }

    - id: load-boundary
      use: io.read_vector
      params: { path: "data/study_area.shp" }

    - id: clip-scene
      use: raster.clip
      params:
        input: "$load-multispectral"
        mask: "$load-boundary"
        crop: true

    # B5 = 近红外波段（Landsat 8 第5波段）
    # B4 = 红波段（Landsat 8 第4波段）
    - id: calc-ndvi
      use: raster.calc
      params:
        input: "$clip-scene"
        expression: "(B5 - B4) / (B5 + B4)"

    - id: save-ndvi
      use: io.write_raster
      params:
        input: "$calc-ndvi"
        path: "output/ndvi.tif"

  outputs:
    ndvi_min: "$calc-ndvi.min"
    ndvi_max: "$calc-ndvi.max"
    ndvi_mean: "$calc-ndvi.mean"

9.6 `raster.stats`：栅格统计

参数

参数	类型	必填	默认值	说明
`input`	raster_info	✅	—	输入栅格数据
`band`	number	❌	`1`	统计的波段编号（1-indexed）

输出

属性	类型	说明
`output`	dict	统计结果字典（与 stats 内容一致）
`stats.min`	float	最小值（忽略 NaN）
`stats.max`	float	最大值
`stats.mean`	float	平均值
`stats.std`	float	标准差
`stats.band`	int	实际统计的波段编号
`stats.pixel_count`	int	参与统计的像素数量

示例

- id: load-dem
  use: io.read_raster
  params: { path: "data/dem.tif" }

- id: dem-stats
  use: raster.stats
  params:
    input: "$load-dem"

outputs:
  elevation_min: "$dem-stats.min"
  elevation_max: "$dem-stats.max"
  elevation_mean: "$dem-stats.mean"

9.7 `raster.contour`：等值线提取

raster.contour 从栅格数据中提取等值线，输出为矢量 GeoDataFrame（LineString 类型）。

参数

参数	类型	必填	默认值	说明
`input`	raster_info	✅	—	输入栅格数据
`interval`	number	✅	—	等值线间隔
`band`	number	❌	`1`	使用的波段编号（1-indexed）
`base`	number	❌	`0`	基准值（等值线从该值开始按 interval 递增）

输出

属性	类型	说明
`output`	GeoDataFrame	等值线矢量数据（LineString），含 `elevation` 字段表示等值
`stats.contour_count`	int	生成的等值线数量
`stats.interval`	number	实际使用的等值线间隔
`stats.min_elevation`	float	最小等值
`stats.max_elevation`	float	最大等值

示例

- id: load-dem
  use: io.read_raster
  params: { path: "data/dem.tif" }

# 每 100 米一条等高线
- id: contour-100m
  use: raster.contour
  params:
    input: "$load-dem"
    interval: 100
    band: 1
    base: 0

- id: save-contours
  use: io.write_vector
  params:
    input: "$contour-100m"
    path: "output/contours_100m.geojson"
    format: "GeoJSON"

9.8 栅格分析综合流水线示例

pipeline:
  name: "遥感影像综合分析"
  description: "对多光谱影像进行裁剪、NDVI 计算、统计、等值线提取"

  variables:
    scene_path: "data/scene_multiband.tif"
    study_area: "data/study_area.shp"

  steps:
    - id: load-scene
      use: io.read_raster
      params: { path: "${scene_path}" }

    - id: load-boundary
      use: io.read_vector
      params: { path: "${study_area}" }

    - id: clip-to-study-area
      use: raster.clip
      params:
        input: "$load-scene"
        mask: "$load-boundary"
        crop: true

    - id: calc-ndvi
      use: raster.calc
      params:
        input: "$clip-to-study-area"
        expression: "(B5 - B4) / (B5 + B4)"

    - id: ndvi-statistics
      use: raster.stats
      params:
        input: "$calc-ndvi"

    - id: save-ndvi
      use: io.write_raster
      params:
        input: "$calc-ndvi"
        path: "output/ndvi.tif"

    # 每 0.1 提取一条 NDVI 等值线
    - id: ndvi-contours
      use: raster.contour
      params:
        input: "$calc-ndvi"
        interval: 0.1
        base: -1.0
      on_error: skip   # 等值线提取失败不中断主流程

    - id: save-contours
      use: io.write_vector
      params:
        input: "$ndvi-contours"
        path: "output/ndvi_contours.geojson"
        format: "GeoJSON"
      when: "$ndvi-contours.contour_count > 0"

  outputs:
    ndvi_min: "$ndvi-statistics.min"
    ndvi_max: "$ndvi-statistics.max"
    ndvi_mean: "$ndvi-statistics.mean"

9.9 本章小结

本章介绍了 5 个栅格分析步骤：

raster.reproject：栅格投影转换，支持 bilinear/cubic/nearest 等重采样方法
raster.clip：用矢量掩膜裁剪栅格，需注意矢量与栅格坐标系一致性
raster.calc：波段数学运算，用 B1/B2… 引用波段，安全 AST 求值
raster.stats：统计波段的 min/max/mean/std 等信息
raster.contour：从栅格提取等值线为矢量数据，含 elevation 属性字段

导航：← 第八章：矢量分析步骤｜第十章：空间分析步骤 →