znlgis 博客

GIS开发与技术分享

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

9.1 概述

栅格分析步骤处理栅格（Raster）数据，如卫星影像、DEM（数字高程模型）、遥感影像等。GeoPipeAgent 基于 Rasterio 和 NumPy 实现栅格处理，所有栅格步骤的输入和输出都使用统一的栅格数据字典格式。

GeoPipeAgent 提供 5 个栅格分析步骤：

Step ID 名称 说明
raster.reproject 栅格投影转换 转换栅格 CRS
raster.clip 栅格裁剪 用矢量范围裁剪栅格
raster.calc 栅格计算 波段数学运算（如 NDVI）
raster.stats 栅格统计 计算波段统计值
raster.contour 等值线生成 从栅格生成等值线矢量

9.2 栅格数据格式

在 GeoPipeAgent 中，栅格数据以字典格式在步骤之间传递：

raster_info = {
    "data": numpy.ndarray,      # shape: (bands, height, width)
    "transform": Affine,        # 仿射变换矩阵
    "crs": CRS,                 # 坐标参考系统
    "profile": dict,            # rasterio 配置信息
}

这个字典由 io.read_raster 创建，并在各栅格步骤之间传递。

9.3 raster.reproject — 栅格投影转换

9.3.1 功能说明

将栅格数据从一个坐标参考系统（CRS）转换到另一个。使用 rasterio.warp.reproject 实现。

9.3.2 参数

参数 类型 必需 默认值 说明
input raster_info ✅ — 输入栅格数据
target_crs string ✅ — 目标 CRS（如 EPSG:3857）
resolution number ❌ — 输出分辨率（可选）

9.3.3 使用示例

- id: reproject-dem
  use: raster.reproject
  params:
    input: "$read-dem.output"
    target_crs: "EPSG:3857"

9.3.4 工作原理

栅格投影转换涉及像素重采样，框架使用 rasterio.warp 模块计算新的变换矩阵和栅格尺寸，然后执行投影转换和重采样。

9.4 raster.clip — 栅格裁剪

9.4.1 功能说明

使用矢量几何范围裁剪栅格数据。裁剪后，栅格范围缩小到矢量几何的外包矩形（或精确裁剪到几何边界）。

9.4.2 参数

参数 类型 必需 默认值 说明
input raster_info ✅ — 输入栅格数据
mask geodataframe ✅ — 裁剪范围（矢量数据）
crop boolean ❌ true 是否裁剪到范围外包矩形
nodata number ❌ 0 范围外的填充值

9.4.3 使用示例

pipeline:
  name: "DEM 裁剪"
  steps:
    - id: read-dem
      use: io.read_raster
      params:
        path: "data/dem.tif"
    - id: read-boundary
      use: io.read_vector
      params:
        path: "data/study_area.shp"
    - id: clip-dem
      use: raster.clip
      params:
        input: "$read-dem.output"
        mask: "$read-boundary.output"
        crop: true
        nodata: -9999
    - id: save
      use: io.write_raster
      params:
        input: "$clip-dem.output"
        path: "output/clipped_dem.tif"

9.5 raster.calc — 栅格计算

9.5.1 功能说明

对栅格波段执行数学运算，是遥感分析中最强大的步骤之一。支持使用波段变量（B1, B2, …）编写 NumPy 表达式。

9.5.2 参数

参数 类型 必需 默认值 说明
input raster_info ✅ — 输入栅格数据
expression string ✅ — 计算表达式

9.5.3 表达式语法

使用 B1, B2, B3, … 引用波段：

表达式 说明
(B4-B3)/(B4+B3) NDVI（归一化植被指数）
B1 + B2 波段加法
B1 * 0.5 波段缩放
np.sqrt(B1**2 + B2**2) 波段运算
np.where(B1 > 0, B1, 0) 条件运算

9.5.4 安全的 np 函数

表达式中允许使用以下 NumPy 函数：

函数 说明
np.abs() 绝对值
np.sqrt() 平方根
np.log() / np.log2() / np.log10() 对数
np.exp() 指数
np.sin() / np.cos() / np.tan() 三角函数
np.arcsin() / np.arccos() / np.arctan() 反三角函数
np.minimum() / np.maximum() 逐元素最小/最大值
np.where() 条件选择
np.clip() 值裁剪
np.nan_to_num() NaN 处理

9.5.5 安全验证

raster.calc 使用 AST 白名单验证 确保表达式安全：

  1. 解析表达式为 AST（抽象语法树）
  2. 遍历所有节点，只允许安全的节点类型
  3. 函数调用只允许 np.* 命名空间中的白名单函数
  4. 变量名只允许波段引用（B1, B2, …）和 np
  5. 拒绝任何不在白名单中的构造

9.5.6 使用示例

计算 NDVI：

pipeline:
  name: "NDVI 计算"
  steps:
    - id: read
      use: io.read_raster
      params:
        path: "data/landsat.tif"
    - id: ndvi
      use: raster.calc
      params:
        input: "$read.output"
        expression: "(B4-B3)/(B4+B3)"
    - id: save
      use: io.write_raster
      params:
        input: "$ndvi.output"
        path: "output/ndvi.tif"
  outputs:
    ndvi_stats: "$ndvi.stats"

波段比值运算：

- id: ratio
  use: raster.calc
  params:
    input: "$read.output"
    expression: "B4 / np.maximum(B3, 0.001)"    # 避免除以零

9.5.7 输出统计

raster.calc 自动计算结果的统计信息：

{
  "expression": "(B4-B3)/(B4+B3)",
  "min": -0.95,
  "max": 0.98,
  "mean": 0.35
}

9.6 raster.stats — 栅格统计

9.6.1 功能说明

计算栅格各波段的统计信息，包括最小值、最大值、均值和标准差。

9.6.2 参数

参数 类型 必需 默认值 说明
input raster_info ✅ — 输入栅格数据

9.6.3 输出

属性 说明
stats.bands 各波段统计列表
stats.bands[n].min 最小值
stats.bands[n].max 最大值
stats.bands[n].mean 均值
stats.bands[n].std 标准差

9.6.4 使用示例

pipeline:
  name: "DEM 统计分析"
  steps:
    - id: read
      use: io.read_raster
      params:
        path: "data/dem.tif"
    - id: stats
      use: raster.stats
      params:
        input: "$read.output"
  outputs:
    statistics: "$stats.stats"

9.7 raster.contour — 等值线生成

9.7.1 功能说明

从栅格数据生成等值线，将栅格中连续的值转换为矢量线要素。常用于：

  • 从 DEM 生成等高线
  • 从温度栅格生成等温线
  • 从降雨数据生成等降水线

9.7.2 参数

参数 类型 必需 默认值 说明
input raster_info ✅ — 输入栅格数据
interval number ✅ — 等值线间隔
band number ❌ 1 使用的波段号

9.7.3 使用示例

从 DEM 生成等高线：

pipeline:
  name: "等高线生成"
  steps:
    - id: read-dem
      use: io.read_raster
      params:
        path: "data/dem.tif"
    - id: contour
      use: raster.contour
      params:
        input: "$read-dem.output"
        interval: 50          # 50 米等高线间隔
    - id: save
      use: io.write_vector
      params:
        input: "$contour.output"
        path: "output/contour_lines.geojson"

9.7.4 实现原理

raster.contour 使用 matplotlib.contour 生成等值线，然后将等值线转换为 Shapely LineString 几何对象，最终封装为 GeoDataFrame。

9.8 栅格分析综合案例

9.8.1 遥感影像 NDVI 分析流水线

pipeline:
  name: "遥感影像 NDVI 分析"
  description: "从 Landsat 影像计算 NDVI 并生成统计报告"

  variables:
    satellite_path: "data/landsat_8.tif"
    boundary_path: "data/study_area.shp"

  steps:
    # 1. 读取卫星影像
    - id: read-image
      use: io.read_raster
      params:
        path: "${satellite_path}"

    # 2. 读取研究区范围
    - id: read-boundary
      use: io.read_vector
      params:
        path: "${boundary_path}"

    # 3. 裁剪到研究区
    - id: clip-image
      use: raster.clip
      params:
        input: "$read-image.output"
        mask: "$read-boundary.output"
        crop: true

    # 4. 计算 NDVI
    - id: calc-ndvi
      use: raster.calc
      params:
        input: "$clip-image.output"
        expression: "(B4-B3)/(B4+B3)"

    # 5. 统计 NDVI
    - id: ndvi-stats
      use: raster.stats
      params:
        input: "$calc-ndvi.output"

    # 6. 保存 NDVI 结果
    - id: save-ndvi
      use: io.write_raster
      params:
        input: "$calc-ndvi.output"
        path: "output/ndvi.tif"

  outputs:
    ndvi_result: "$save-ndvi.output"
    ndvi_statistics: "$ndvi-stats.stats"
    image_info: "$read-image.stats"

9.8.2 DEM 分析与等高线生成

pipeline:
  name: "DEM 分析"
  steps:
    - id: read-dem
      use: io.read_raster
      params:
        path: "data/dem.tif"

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

    - id: contour-50m
      use: raster.contour
      params:
        input: "$read-dem.output"
        interval: 50

    - id: save-contour
      use: io.write_vector
      params:
        input: "$contour-50m.output"
        path: "output/contour_50m.geojson"

  outputs:
    elevation_stats: "$dem-stats.stats"
    contour_features: "$contour-50m.feature_count"
↑