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

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

9.1 步骤总览

9.2 栅格数据结构说明

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

从 io.read_raster 读取后，这个字典通过步骤引用传递，metadata 字段也会额外存储 transform、crs、dtype、nodata、count、width、height、bounds。

9.3 raster.reproject：栅格投影转换

参数

重采样方法选择

示例

- 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"
    metadata: "$reproject-dem.metadata"

9.4 raster.clip：栅格裁剪

参数

示例

- 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"
    metadata: "$clip-satellite.metadata"

9.5 raster.calc：栅格计算

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

参数

安全性

表达式通过 safe_eval.py 中的 AST 白名单求值，只允许数值运算，禁止任意函数调用，防止代码注入。

常用表达式示例

# 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"

# 归一化
expression: "(B1 - B1.min()) / (B1.max() - B1.min())"  # 注意：B1.min()语法不支持，需要用数值

注意：表达式中只能使用 B1、B2... 形式引用波段，不支持 band1、Band1 等写法。表达式通过 AST 白名单求值，只允许基本数学运算符（+, -, *, /, **）和括号。

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"
        metadata: "$calc-ndvi.metadata"

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

9.6 raster.stats：栅格统计

参数

输出

示例

- 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 类型）。

参数

输出

示例

- 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"
        metadata: "$calc-ndvi.metadata"

    # 每 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 个栅格分析步骤：

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

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