6.1 栅格与地图代数：把图层当数算

6.1 栅格与地图代数：把图层当数算

本节摘要：地图代数是栅格叠加的算术基础——同规格栅格逐像元对位运算，加减乘除、布尔逻辑、条件分支直接写成公式。本节讲清运算前提（像元大小、原点、范围、坐标系四要素一致）、重分类这一"统一刻度"的关键步骤、栅格计算器的语法要点，以及 NoData 在运算中的传播规则。

当叠加变成算术

第 5 章做相交，软件要计算两个多边形的交线交点，输出仍是多边形。栅格世界里没有这些几何谈判：两张 30 米分辨率的栅格，像元天然一一对齐，"叠加"就是第 (i,j) 格与第 (i,j) 格相加，遍历全图。三张因子图合成分数图，本质上是一个三维数组沿第三维求加权和——线性代数里最普通的运算。

这种朴素带来速度优势：百万像元的加权叠加毫秒级完成，同样范围的高精度矢量求交可能要分钟级。代价在第 3.1 节已经讲过：边界被量化到像元格。栅格叠加的适用判断因此很清晰——因子是连续变量（坡度、距离、温度）且需要全区域逐点评分时用栅格；因子是边界分明的法定对象（权属、规划红线）且输出要求边界精确时用矢量。

四要素一致是运算前提：两张栅格想逐像元对位，必须像元大小相同、左上角原点对齐、范围兼容、坐标系一致。四者有一不齐，软件要么报错，要么按最近邻规则强行取值——后者更危险，因为它静默地错位相加，输出的图看起来正常，数字全是错的。工程上的一致性做法是"锚栅格"制度：项目定一张基准栅格（通常是主 DEM），所有衍生栅格生成时都在环境里指定它为捕捉栅格，像元规格从源头统一。

图：逐像元运算的对应关系

重分类：统一刻度的关键一步

坡度单位是度，离水距离单位是米，土层厚度单位是厘米——三张图直接相加毫无意义。重分类（Reclassify）把每种量纲换算到统一的分数刻度：坡度 0 到 3 度记 9 分、3 到 8 度记 7 分……离水 0 到 200 米记 9 分、200 到 500 米记 7 分……分档与赋值由专业知识决定（种植手册、规划规范、专家经验），重分类工具只负责执行换算。

import arcpy
from arcpy.sa import Reclassify, RemapRange
arcpy.env.workspace = r"K:/gisdata/suit_tea.gdb"
arcpy.CheckOutExtension("Spatial")                      # 栅格运算需要空间分析许可
arcpy.env.snapRaster = "dem30"                          # 锚栅格：像元规格统一到 DEM
arcpy.env.cellSize = "dem30"
# 坡度重分类：度数换 1 到 9 分（分档依据种植技术规程）
slope_score = Reclassify("slope_deg", "VALUE",
    RemapRange([[0, 3, 9], [3, 8, 7], [8, 15, 5], [15, 25, 3], [25, 90, 1]]))
slope_score.save("slope_score")
print("坡度分栅格完成:", slope_score.minimum, "到", slope_score.maximum, "分")
# 输出: 坡度分栅格完成: 1.0 到 9.0 分
​

环境里的 snapRaster 与 cellSize 就是锚栅格制度的落地：此后本工程所有衍生栅格自动与 DEM 像元对齐，四要素一致性由环境保证而非人眼比对的运气。

栅格计算器：公式即分析

栅格计算器（Raster Calculator）是地图代数的交互前台，语法就是 Python 表达式：栅格用双引号包名，运算符直接写，函数从空间分析模块调用。几个高频模式：

import arcpy
from arcpy.sa import Raster
arcpy.CheckOutExtension("Spatial")
A, B, C = Raster("slope_score"), Raster("water_score"), Raster("soil_score")
# 一、加权求和：三因子合成为综合分
combo = 0.5 * A + 0.3 * B + 0.2 * C
combo.save("suit_combo")
# 二、布尔门槛：一票否决式（坡度分低于3的地块无论其他多好都出局）
gate = (A >= 3) & (B >= 3) & (C >= 3)
gate.save("suit_gate")          # 输出 1 与 0：可行 与 否决
# 三、条件合成：门槛内保留综合分 门槛外归零
from arcpy.sa import Con
final = Con(gate, combo, 0)
final.save("suit_final")
print("综合分范围:", final.minimum, "到", final.maximum)
# 输出: 综合分范围: 0 到 8.6
​

三段表达式演示了地图代数的三个层次：算术（加权）、逻辑（布尔与或非）、控制流（Con 条件函数，相当于逐像元的 if）。布尔栅格的 1 与 0 让"并且"用 &、"或者"用 |、"非"用 ~，与编程直觉完全一致。

界面上的栅格计算器写法相同：双击图层名入表达式、选运算符、设输出名。但脚本版有个不可替代的优势——表达式进版本管理，权重 0.5、0.3、0.2 的每次调整都有记录，这又是第 3 章谱系纪律的延续。

NoData：运算里的黑洞

NoData 是栅格里的特殊值，表示"这个格子没有信息"（云下区域、研究区外）。它在地图代数中有传染性：NoData + 5 = NoData——任何一格是 NoData，参与加法的结果就是 NoData。这常常符合预期（信息不全的格子确实不该有分），但也可能吞掉大片结果：DEM 边缘一圈 NoData 让最终适宜性图四周全部空白。

处理策略按意图选：希望传染（信息不完整即无结论）就顺其自然；希望把 NoData 当零处理（缺因子按最差计）用 Con(IsNull(ras), 0, ras) 先转换；希望在多因子中"跳过缺失因子"则需要更复杂的 Con 嵌套。关键是先弄清 NoData 从哪来（原始缺测？范围裁剪？重分类漏段？），再决定策略——尤其重分类的 RemapRange 漏掉某段值域时，那一段会静默变 NoData，这是新手最常踩的静默陷阱。

⚠️ 常见坑：两张"看起来一样"的栅格相加，输出莫名缺一块。九成是像元原点错半格或范围不同，软件按捕捉规则错位取值。运算前打印 Raster(...).extent 与 meanCellWidth 对比，十秒钟避免一场空欢喜。

💡 关键直觉：把地图代数当电子表格用。因子图是列，像元是行，公式写在计算器里向下填充到百万行——Excel 里你会检查两列是否对齐行，栅格运算前就检查像元是否对齐格。

要点回顾

• 对位即运算：栅格叠加是格子对格子的算术，四要素一致是前提，锚栅格制度是保障
• 重分类统一刻度：量纲各异先换 1 到 9 分，分档依据专业知识而非软件默认
• 三层语法：算术加权、布尔逻辑、Con 条件分支，公式即分析逻辑
• NoData 传染：加法中缺失即缺失，弄清来源再决定转换策略
• 脚本优于界面：表达式进版本管理，权重变更有谱系

算术备齐，下一节把方法论补全——因子怎么选、权重怎么定、结果怎么验，完成一场真正的适宜性合奏。