GeoScene Pro栅格叠加取最大值：单元统计与NoData实战

1. 需求拆解：多张栅格叠在一起，为什么要取最大值

手头攒了七八张同一区域的栅格，格式一样、范围差不多、像元大小也一样，领导一句"把每块地最旺的那期挑出来"，问题就来了：这些图叠在一起，同一个像元位置上落了好几个值，到底该听谁的？这就是相同位置栅格获取最大值要解决的事。说白了，就是把多张栅格的像元按空间位置一一对齐，然后每个位置只保留数值最大的那个，最后输出一张新的单波段栅格。听起来简单，但它背后牵扯的东西一点不少——坐标系对不对、像元格网有没有错位、NoData怎么算、统计范围取交集还是并集，任何一项没想清楚，结果就是错的，而且错得很隐蔽，你不做交叉验证根本看不出来。

这篇内容面向的是已经能打开软件、能找到工具箱，但对栅格叠加的"规矩"还不太有底的同行；也适合那些做过一两次、结果总是对不上、想搞清楚问题出在哪的人。我会拿GeoScene Pro里的**单元统计（Cell Statistics）**工具为主线，把参数一个个拆开讲清楚，再补上栅格计算器和脚本两条备选路线，最后把我在实际项目里踩过的坑整理成一张速查表。整篇不涉及具体业务数据，只讲方法、讲参数、讲判断逻辑，你换成自己的数据就能直接跑。

1.1 三个"取最大值"的典型判断依据

为什么偏偏是最大值，不是平均值、不是求和？这背后其实是业务语义在决定统计量。我把常见的判断逻辑归成三类，你对号入座就行。

第一类是**"取最优状态"**。多期植被指数是最典型的例子：一年十二个月的NDVI，云污染、阴影、雪覆盖会让某些月份的数值偏低甚至为负，但作物在生长旺季的那一期一定是全年最高的。这时候逐像元取最大值，等于自动把"最健康的那一天"挑出来，顺带把云和阴影的影响压下去。这就是遥感里常说的最大值合成法（MVC），是生成月度、年度植被指数产品的标准动作，不是什么野路子。

第二类是**"取最不利情景"**。听着跟上面相反，但逻辑一致——比如多套洪水淹没模拟方案叠在一起，你要的是每个位置淹得最深的那一套，用来定防洪堤高度；又比如多个气象情景下的最高温度，用来做热风险评估。这时候最大值代表的是"最坏情况"，是安全设计的依据。

第三类是**"取信号最强"**。多传感器、多时相的观测数据，同一个位置可能被反复扫到好几次，取最大值相当于选信噪比最好的那一次观测。

注意：这三类场景的共同点是"最大值本身有意义"。如果业务上问的是"总量"或者"平均水平"，那就该用求和或平均值，硬套最大值只会得出荒谬结论。工具选错比参数填错更致命。

1.2 不适合取最大值的场景，先排除掉

动手之前先花两分钟做个反向检查，能省掉后面大量返工。有几种情况看着像能用最大值，实际上不该用。

一种是量纲不一致的数据。比如一张是反射率（0到1之间），一张是DN值（0到255），你直接取最大值，结果必然是DN值全线胜出，这张图等于废了。要么先统一量纲，要么就别混着算。

另一种是物理意义冲突的数据。比如地表温度和 NDVI 放在一起取最大值，这在物理上讲不通，纯属数字游戏。统计的前提是这些栅格描述的是同一个可比的量。

还有一种是栅格分辨率差异过大。一张30米、一张10米，直接叠上去，工具会自动重采样到某个统一分辨率，但重采样本身就是一次信息损失和人为插值，取出来的"最大值"可能是插值插出来的，不是真实观测。这种时候应该先把粗分辨率数据做降尺度处理，或者干脆分开算，不要硬凑。

最后一种，时间跨度太长的数据。十年前和今年的同一区域栅格取最大值，土地利用可能早就变了，出来的图既不是十年前的也不是今年的，是个不存在的混合体。做时间序列分析前，先确认地物没有发生结构性变化，否则结论没有解释力。

2. 工具选型：单元统计、栅格计算器、脚本，我为什么选单元统计

在 GeoScene Pro 里实现"相同位置取最大值"，我知道的至少有三条路能走通：单元统计工具、栅格计算器写表达式、Python 脚本调用。三条路都能出结果，但适用场景差别很大，选错了不是跑不出来，而是跑得又慢又容易出错。

2.1 三条技术路线的横向对比

我先把三者的关键差异摆出来，你按自己的数据量和习惯挑。

栅格计算器那条路，写起来大概是Con(IsNull("r1"), "r2", Con("r1" > "r2", "r1", "r2"))这种层层嵌套的形式。两张图还行，三张图就开始眼晕，五张图基本就是噩梦，而且括号错一个就报错，排查起来极其痛苦。它适合的场景是"顺手算一下"或者需要跟其他运算混着来，不适合当主力方案。

脚本那条路适合批量和自动化，尤其是你要对几十个文件夹、上百个文件重复做同一件事的时候，写一次脚本顶手动干一天。但它有个门槛：环境配置、模块导入、路径处理，对不常写代码的人来说，光是让脚本跑起来就要折腾一阵。

单元统计工具是我在绝大多数情况下的首选，理由很实在：输入直接拖多值列表，统计类型下拉框里选"最大值"，NoData 用一个复选框控制，其余交给环境设置。参数少、语义清晰、结果可复现，把我的注意力全部留给"数据本身对不对"，而不是"表达式写没写对"。

一个容易忽略的前提：单元统计属于空间分析扩展模块的工具。如果你在系统工具箱里翻不到"单元统计"，先去看工程设置里的许可配置，确认空间分析模块已经勾选启用。这一步没做，后面所有操作都无从谈起。

2.2 单元统计工具到底做了什么（原理拆解）

把原理讲透，参数就好理解了。想象你手上有三张完全透明的方格纸，每张纸上每个格子里写了一个数字，三张纸严丝合缝地摞在一起。现在你从上往下看每一个格子，看到三个数字，规则是"只留最大的那个"，然后把留下的数字写在一张新纸上。这张新纸就是输出栅格。

关键在于"严丝合缝"这四个字。三张方格纸的格子必须完全对齐，误差为0，否则你看到的三个数字其实来自三个不同的地理位置，取最大值就失去了意义。这就是为什么工具运行前会自动做一次格网对齐检查——如果输入的栅格像元大小或原点位置不一致，工具要么报错，要么按环境设置里的规则强制对齐，而强制对齐的过程就是重采样，会引入误差。

所以单元统计真正的技术含量不在"取最大"这个动作本身，而在于输入栅格是否处在同一个空间格网上。这个前提满足了，工具就是傻瓜式的；不满足，你得到的是一张看起来很正常、实际上每个像元都错位的图。

工具内部对每个像元位置的处理逻辑是这样的：读取所有输入栅格在该位置的像元值，根据"NoData 处理"设置决定哪些值参与统计，然后执行最大值运算，写入输出。整个过程是逐像元独立进行的，像元之间互不影响，这也意味着它天然支持并行分块处理，数据量大的时候可以靠切块来提速。

3. 完整实操：从数据检查到出图的每一步

下面进入实操。我按"检查—配置—运行—验证"的顺序走一遍，每一步的做法和理由都写清楚，你照着做能少走弯路。

3.1 前置检查：坐标系、像元大小、格网对齐

正式打开工具之前，务必先做三项检查，这三项决定了后面会不会白干。

第一项，坐标系是否一致，且是否为投影坐标系。地理坐标系（经纬度）下做栅格运算是个经典陷阱：经线和纬线方向的像元实际地面尺寸不一样，越往高纬度差异越大，算出来的"最大值"位置可能偏移好几百米。正确做法是先统一投影到适合该区域的投影坐标系，再叠加。检查方法很简单，在内容列表里逐张看属性，或者在工具箱里跑一次批量投影，把结果作为新的输入。

第二项，像元大小是否一致。不一致的时候，工具默认会按环境设置里的像元大小规则处理，通常是取"最大输入"或"最小输入"。这里有个取舍：取最小像元意味着粗分辨率数据被双线性或最邻近插值放大，图看着细了，但细节是编的；取最大像元意味着细分辨率数据被聚合，细节被抹掉。我的建议是，如果分辨率差异在两倍以内，可以接受重采样；超过两倍，先考虑降尺度或分区域处理。

第三项，格网原点是否对齐。这一项最隐蔽。两张图像元大小都是30米，但一个原点在 (0,0)，另一个在 (15,15)，它们的格子在空间上是错开半个像元的。这时候工具会自动重采样，结果偏差不会立刻显现，但做逐像元对比时就会发现问题。

实操心得：在环境设置里指定一个"捕捉栅格"（Snap Raster），把它设为你所有输入里最可靠的那一张，其余栅格都会向它的格网对齐。这一个设置能解决八成的错位问题，比事后排查省事得多。

处理范围我一般设为"输入的并集"，这样不会因为某张图边界稍小就把边缘数据裁掉；如果业务只需要共同覆盖区域，那就设"输入的交集"。这两者的差别在后续做面积统计时非常明显，一定要事先想清楚。

3.2 工具参数逐个拆解

工具位置在系统工具箱的空间分析工具下、局部分析类别里，名称就是单元统计。打开后能看到四个主要参数和一个环境设置入口，我一个一个说。

输入栅格：这是一个多值参数，点右侧的加号可以一次选多个文件，也可以直接把内容列表里的图层拖进来。顺序不影响结果，最大值运算跟先后无关。但有个细节要注意：如果输入里混了多波段栅格，默认行为是只处理第一个波段，除非你勾选了"处理为多波段"。这一项漏掉，输出可能少了好几个波段，你还以为是工具出错了。

统计类型：下拉框里有平均值、多数、最大值、中位数、最小值、少数、范围、标准差、总和、变异度十种。选"最大值"就行。这里顺便提一句，"范围"是最大值减最小值，"变异度"是不同值的个数，这几个统计量在做数据质量检查时特别有用，后面会讲到。

忽略 NoData：这是个复选框,勾选和不勾选的结果差别巨大。勾选表示"该位置只要有一个输入有值就参与统计"；不勾选表示"该位置只要有一个输入是 NoData，输出就是 NoData"。选哪个取决于你的数据特性：如果 NoData 代表"没有观测"而不是"观测值为零"，那必须勾选，否则有效区域会被一点点啃掉；如果 NoData 代表"该位置物理上不存在"，那可以不勾选，保持严格一致。

输出栅格：指定路径和文件名，格式建议选 TIFF 或者写入文件地理数据库。格式这块在后面性能部分还会展开说。

环境设置里有几个参数值得单独调：处理范围、像元大小、捕捉栅格、输出坐标系统。前三个前面已经讲过，输出坐标系统建议显式指定，别依赖默认值，避免输出跟输入坐标不一致导致后续流程出错。

参数就是这么几个，没有花哨的东西。但参数的语义必须跟数据的语义对上，这是我在无数次返工里学到的最贵的一课。

3.3 结果验证：怎么确认最大值得对了

工具跑完弹出"成功"两个字，不代表结果是对的。我一般做两步验证。

第一步，抽样点检。用识别工具在输出图上点几个位置，记下值；然后切回每张输入图，在同一点位分别读值，手动比对最大的那个是否跟输出一致。点位要挑得有讲究：至少一个在数据重叠核心区，一个在边缘地带，一个在 NoData 交界处。边缘和交界处是问题高发区，核心区通常没事。

第二步，差值检查。用栅格计算器做一次输出 - 某一张输入，结果里如果出现负数，说明输出在该位置小于这张输入，逻辑上就矛盾了，一定有问题。反过来也可以用输出 >= 每张输入的组合判断，逐张验证。这个检查跑起来很快，但能抓住几乎所有逻辑错误。

第三步，看统计信息。打开输出栅格的最小值、最大值、平均值，跟输入图层对比。输出的最大值应该等于所有输入最大值中的最大者，输出的最小值应该不小于所有输入最小值中的最小者。这两个不等式不成立，说明统计过程有问题。

一个小技巧：如果输入栅格数量很多，逐个比对太累，可以先用"范围"统计类型跑一遍，得到最大值减最小值的分布图，异常区域会非常显眼——差值异常大或者异常小的位置，往往就是数据质量问题所在。

4. 问题排查实录：我踩过的六个坑

这一节是我最想写的部分，因为工具本身不难，难的是数据总有各种意想不到的状况。下面这几类问题，我几乎每个项目都会遇到至少一种。

4.1 NoData 是最大的变量

NoData 在栅格里的地位很特殊，它不是一个数值，而是"此处无数值"的标记。工具怎么对待它，直接决定输出的形状。

第一个坑，边界被啃掉。假设三张图覆盖范围略有差异，A图是完整矩形，B图缺了右下角，C图缺了左上角。如果不勾选"忽略 NoData"，输出图只保留三张图都有的区域，面积比你预期小一圈，而你如果没注意，拿去做面积统计，结果就全错了。勾选之后，缺角位置会用剩下两张图的值来算，覆盖范围才符合直觉。

第二个坑，0 和 NoData 混为一谈。有些数据源把"无观测"直接写成0，这是数据生产环节的错误，但到了分析环节就变成灾难：0 会被当成真实数值参与最大值运算，虽然它很难胜出（除非所有值都是负数），但如果你做的是最小值或者平均值，影响就大了。遇到这种情况必须先做一次条件替换，把 0 转成 NoData，或者在读取时用栅格计算器做掩膜。

第三个坑，输出全是 NoData。这种情况通常是坐标系没对上，或者栅格根本没重叠——你以为它们同区域，实际上一个在另一个的几千公里外。排查方法是在地图里同时加载所有输入，看一眼是不是真的叠在一起。这个错误听起来很蠢，但在批量处理时，路径配错、文件张冠李戴，真的很常见。

4.2 像元错位：最隐蔽的错误

像元错位不会报错，不会警告，输出图看着漂漂亮亮，但每个位置的值都可能来自错误的地理坐标。我把它称为"沉默的错误"。

产生原因有三：一是输入栅格像元大小不同但差距很小，工具自动重采样时原点偏移；二是某张图做过裁剪或投影转换，输出格网跟原始格网错开半个像元；三是数据来源不同，一个是国产卫星产品，一个是国际产品，两者的格网定义规则本来就不一样。

排查方法我总结了一个笨但有效的：在两张图上各画一条相同的线段，读出线上所有像元值，逐点比对。如果同一个地理位置上，两张图的值呈现规律性偏移（比如总是差一个像元），那就是错位了。更快的办法是用栅格计算器做一次减法图A - 图B，如果两图本来高度相关，正确的减法结果应该是平滑的、有空间结构的；错位的结果则会出现大量高频噪声和边缘伪影，一眼能看出来。

解决手段就一个：在环境设置里明确指定捕捉栅格，强制所有输入向同一个格网对齐。宁可多花时间做一次规范的重采样，也不要把错位的数据带进下游分析。

4.3 性能与数据组织

栅格运算的耗时跟数据量基本是线性关系，但有几个因素会让它变成指数级的灾难。

文件格式。同样是 TIFF，带金字塔和不带金字塔的读取速度差好几倍。大数据量时，事先构建金字塔和统计信息，能让工具运行时间砍掉一大半。这个操作在图层属性或者工具箱里有现成命令，跑一次也就几十秒。

中间文件堆积。每次运行都在工作目录留下一堆临时文件，磁盘满了之后工具会莫名其妙失败。我的习惯是每个项目单独建一个工作目录，用完就清，别让临时文件满世界乱放。

并行处理。单元统计是逐像元的独立运算，天然适合分块并行。软件的环境设置里通常有并行处理因子这个选项，默认值不一定吃满你的机器，可以手动调高。不过要注意，并行度高了对内存要求也高，机器配置一般的话，调到跟物理核心数持平就行，别贪多。

分块处理。如果单张栅格就超过了内存能承载的规模，可以把研究区分成若干块，逐块统计再拼接。这个方法比较土，但很稳，尤其是在跑几十张高分影像的时候。

4.4 常见问题速查表

我把上面这些问题的表现和处置办法整理成表，遇到对不上的情况可以直接查。

5. 进阶玩法：把这件事做成一键可复用

单次操作学会了，接下来该考虑的是"下次怎么更快"。毕竟业务里这种统计往往不是做一次，而是每个月、每个季度都要来一遍。

5.1 用模型构建器封装成工具

模型构建器是把操作步骤固化下来的最省事方式。思路很简单：把单元统计工具拖进画布，把输入栅格设为模型参数，把输出路径也设为参数，然后保存。以后每次用，就在模型工具对话框里填两个参数，点运行，完事。

模型里还能加前置步骤，比如自动做投影转换、自动建金字塔、自动把 0 转成 NoData。这一串动作封装起来，新人拿到模型也能直接跑，不会漏步骤。我在团队里推进标准化的时候，这一招比写文档管用得多——文档没人看，工具是天天要用的。

有个细节要留意：模型里的输入参数最好设计成"文件夹"或者"栅格列表"类型，这样一次能传多个文件，而不是一个一个填。参数类型设置不对，用起来会很难受。

5.2 用 Python 脚本批量跑

脚本适合真正的大批量场景。核心代码就几行，逻辑是把输入文件列表丢进去，指定统计类型，输出保存。

import arcpy from arcpy.sa import CellStatistics arcpy.env.workspace = r"D:\project\input" arcpy.env.overwriteOutput = True arcpy.env.snapRaster = r"D:\project\input\ref_grid.tif" inputs = arcpy.ListRasters("*", "TIF") result = CellStatistics(inputs, "MAXIMUM", "DATA") result.save(r"D:\project\output\max_result.tif")

这段代码里有两个参数值得说明。第一个是统计类型"MAXIMUM"，换成"MINIMUM"、"MEAN"、"SUM"就是别的统计量，改一个字符串的事。第二个是"DATA"，对应工具里的"忽略 NoData"勾选状态；如果想严格不忽略，改成"NODATA"。这两个位置我见过不少人写错，结果跑出来跟界面操作对不上还找不着原因。

环境设置那两行也别省。snapRaster指定捕捉栅格，能避免格网错位；overwriteOutput打开覆盖，省得每次报文件已存在的错。至于模块名，不同版本的工具包命名可能略有差异，如果你的环境里导入报错，去软件的帮助文档里确认一下当前的模块名和函数签名，以官方文档为准。

5.3 统计量组合的延伸思路

单元统计工具里有十个统计类型，最大值只是其中一个。既然流程都跑通了，顺手把其他几个也用起来，能挖出不少信息。

比如最大值加最小值，两者相减就是"范围"，能直观反映一个区域内数值的波动幅度。波动大的地方，往往是数据质量有问题，或者地物本身变化剧烈，值得单独去看看。

再比如变异度，它统计的是该位置不同取值的个数。假设你有十二期数据，某个像元的变异度是 1，说明十二期里这个位置的值完全没变——这要么是稳定的水体、建筑，要么是数据本身有问题（比如被同一张底图填充了）。变异度异常低或者异常高的区域，都是数据质检的重点目标。

还有一个组合思路：先取最大值，再取最小值，两张图做成假彩色合成，能一眼看出"高值区"和"低值区"的空间分布格局。这个图拿去做汇报比单张最大值图有说服力得多。

我在实际项目里的体会是，工具的价值往往不在它最直接的那个功能上，而在于你把它的输出当成下一步的输入，一层层往下推。最大值只是这条链的第一环，后面还能接分区统计、重分类、变化检测，路很长。

最后分享一个小习惯：每次做完这类逐像元统计，我都会把参数配置截图存档，跟输出文件放在一起。过半年再回头查，不用重新推演当初到底勾没勾那个复选框，看一眼图就明白了。这个习惯帮我省下的返工时间，比学任何技巧都值。