GeoScene Pro像元统计：栅格最大值合成与对齐避坑

做遥感影像处理、地形分析或者多期监测的人，基本都会碰到同一个动作：手上叠着一堆同一区域、不同时相、不同传感器的栅格，要在每一个像元位置上把这些图层里的值挑出最大的那个，合成一张新的栅格。这件事在 GeoScene Pro 里有个非常直接的名字——相同位置栅格获取最大值，工具叫“像元统计”。听起来像是一键操作，但我见过太多人第一次点下去就得到一张全黑、全 NoData 或者整体错位半个像元的图，然后开始怀疑数据坏了。其实问题几乎都出在数据准备和环境设置上，工具本身反而是最简单的一环。这篇就按我自己的操作习惯，把这件事从头到尾拆开讲一遍：它到底在算什么、数据要先做哪些处理、工具参数怎么填、结果怎么验、出了问题怎么查。适合刚接触栅格分析的入门者照着抄，也适合用了一段时间但结果总对不上的朋友对照排查。

1. 先搞清楚“相同位置取最大值”到底在算什么

1.1 三种业务场景，需求其实不完全一样

“取最大值”这三个字在不同的业务里指向的逻辑差别挺大，动手前先对号入座，能省掉后面一半的返工。

第一种是多时相最大值合成，也是最常见的。比如一年十二个月的植被指数，每个像元取全年最大值，用来代表这个位置生长季最好那一天的状况，规避云、雪、阴影造成的低值污染。第二种是多源数据取最大值，比如同一区域有来自不同传感器的多幅影像，或者几个不同模型输出的强度栅格，取最大值相当于取“最不利情况”或“最强响应”，常见于风险叠加、信号强度合成。第三种是条件筛选型的最大值，比如只有满足某个阈值条件的图层才参与比较，本质是先做了掩膜再取最大。

这三种场景在工具层面是同一个操作，但数据预处理的重点完全不同：第一种最怕 NoData 处理不当，第二种最怕空间参考不一致，第三种最怕掩膜逻辑写反。我一般会先问自己一句——我要的是“所有图层的最大”，还是“有效图层中的最大”？这句话决定了后面那个勾选框怎么点。

顺带把一个特别容易混淆的工具拎出来对比。很多人搜“最大值”会搜到“焦点统计”，那里也有最大值选项，但它算的是同一张图上某个像元周围邻域内的最大值；而“像元统计”算的是多张图层在同一个像元位置上的最大值。前者是空间邻域，后者是图层叠加。这两个如果搞反了，结果看着还挺像那么回事，实际上完全不是一回事。

1.2 “相同位置”这四个字，比工具本身苛刻得多

真正让结果出问题的往往不是统计逻辑，而是“相同位置”这个前提根本没被满足。栅格数据要比大小，前提是每个图层在同一个行列号上代表地面上同一块地方。这要求三件事同时成立：坐标系一致、像元大小一致、网格原点对齐。

坐标系不一致是最好发现的，打开图层属性看一遍就知道。麻烦的是后两个：两个图层坐标系完全一样、像元大小都是 10 米，但网格原点差了 3 米，那么在叠加时每一个像元覆盖的地面范围都不一样，取出来的“最大值”其实是在两个错位的地块上比较，结果会有系统性偏差，而这种偏差在视觉上几乎看不出来，尤其在大范围、低对比度的影像上。

我见过一个典型情况：两批影像由不同渠道提供，一批是按整幅影像左上角对齐的，另一批是按标准格网切片的，两者分辨率相同但原点不同。直接跑像元统计，输出看起来没问题，但拿去做后续的时序分析时，早春的数值总是偏高，查了两天才定位到是网格错位导致的混合像元。

判断方法很土但很有效：把两个图层的显示范围都缩放到同一视图，用识别工具在几个特征明显的位置点一下，看两个图层的像元值变化边界是否严丝合缝地重合。另一种更严谨的做法是看两个图层的范围值（Extent），如果上下左右的坐标差不都是像元大小的整数倍，那就是没对齐。

1.3 为什么不推荐一开始就写栅格计算器

“取最大值”最直觉的写法是在栅格计算器里用条件判断或者调用统计函数。技术上完全可行，写一行表达式也能出结果，但我不建议把它作为默认方案，原因有三个。

一是NoData 的行为不直观。表达式的默认行为是 NoData 传播，也就是任何一个输入图层在那个位置是 NoData，整个结果就是 NoData。而最大值合成最讨厌的就是这个——某一期影像上有一小块云被掩掉了，那个位置明明还有其他十一期有效数据，结果却整片丢了。要在表达式里绕开它，得额外套条件函数，越写越长。

二是参数散落在表达式里，不好交接。一个工具对话框里填的东西，别人打开一看就知道设了什么；一串嵌套表达式，过两个月连自己都读不顺。

三是批处理不方便。像元统计作为标准工具，可以直接扔进模型构建器或者脚本里循环，输入图层列表是动态的；表达式里的图层名是硬编码的，换一批数据就得重写。

不是说计算器没用，它在做条件最大值、加权最大值这种“带业务逻辑”的场景里不可替代。但纯取最大值这种标准动作，交给标准工具，稳妥得多。

2. 数据准备：坐标系、像元对齐与 NoData 这三道坎

2.1 投影统一：宁可多花十分钟，也别事后返工

做叠加分析，我习惯先把所有输入的坐标系统一到同一个投影坐标系下，而不是地理坐标系。原因很实际：地理坐标系下单位是度，像元大小是零点几度，一来不好判断对齐没对齐，二来在不同纬度上“一度”对应的实际距离差很多，一个像元的面积会随纬度变化，取最大值这件事在物理意义上就不成立了。

统一的方式有两种：对原始数据重新投影，或者对原始数据重采样。区别在于前者是针对坐标系本身做变换，后者是在已有坐标系内改变像元大小。如果只是坐标系不同、像元大小刚好一致，做一次投影转换就够；如果连像元大小都不同，那必然要重采样，这时候就得做选择。

重采样方法的选择上，我自己的原则是这样：分类数据（土地利用、掩膜、等级栅格）一律用最近邻，保住原值不被插值污染；连续数据（温度、植被指数、高程）可以考虑双线性或三次卷积，但做最大值合成时我更倾向最近邻，因为最大值是一个排序统计量，插值会造出原始数据里不存在的数值，等于凭空多了一些“最大值”候选，结果可能比所有输入的真实最大值还大，这在数据严谨性上说不通。

分辨率统一时，取哪一个像元大小也有讲究。如果几个图层分辨率不同，我一般选中间值偏细的那个，别直接取最粗的。取最粗虽然能减少数据量，但会把精细图层的信息抹平；取最细则会造成粗分辨率图层被大量复制值填充，看起来精细其实是假精细。实际操作里，如果两个分辨率差距在两倍以内，统一到细的那个可以接受；差距超过三倍，我宁可先分开算再决策，也不做一次性的强行对齐。

注意：重采样会永久改变像元值，务必输出为新文件，不要覆盖原始数据。原始数据在后续核对时是唯一的参照物。

2.2 捕捉栅格与环境设置：一个勾选决定成败

GeoScene Pro 的栅格工具都受“环境设置”控制，这是最容易被忽略、也最容易导致结果错位的地方。三个设置项必须动：处理范围、捕捉栅格、像元大小。

处理范围我通常设为“与图层XXX相同”，选定其中一个输入图层作为基准，而不是用“所有图层的交集”或“并集”。用交集会把边缘一圈本来就该有值的地方裁掉；用并集则会造出一圈全是 NoData 的外框，后续再做统计时容易误判。

捕捉栅格是把结果网格的原点焊死在某个参照图层上，保证输出和参照图层逐像元对齐。这一项如果留空，系统会自己算一个范围，有可能和输入都错开。我一般就把第一个输入图层设为捕捉栅格。

像元大小设置在输入图层分辨率不一致时生效，填的时候要跟捕捉栅格保持一致，否则会出现“对齐了原点但对不齐格网”的尴尬情况。

这三项设完，跑出来的结果才有可能和输入严丝合缝。我的习惯是把这套环境设置存成一个自定义环境，后面一系列操作都复用它，避免每次手动填错。

2.3 NoData 的处理哲学：忽略还是传播

这是整个流程里最需要动脑子的一个决策。工具面板上的“忽略 NoData”勾选框，含义是：在某个像元位置上，如果部分输入图层是 NoData，是否要让剩下的有效值继续参与统计。

勾选（忽略 NoData）——只要有一个图层在该位置有有效值，就输出这些有效值的最大值。这是多时相最大值合成的标准做法，也是绝大多数场景该选的。云掩膜、条带缺失、传感器边缘无效区这些局部问题，都能被其它时相的有效数据补上。

不勾选（NoData 传播）——任何一个图层在该位置无效，整个输出就是 NoData。这种模式适合“必须所有输入都齐备才有意义”的场景，比如多条件同时满足的判断，或者做数据完整性检查。

我踩过的一个坑是这样的：做某区域一年的最大值合成，输入只有三期影像，其中一期在东北角有一块被掩掉。勾选了忽略 NoData 之后，那块区域的输出值来自剩下两期。结果看着没问题，但后来做同比分析时发现，那块区域的值整体偏低——因为剩下两期恰好都是低值期。这说明“忽略 NoData”虽然在数据完整性上更优，但会引入时间代表性偏差：参与统计的样本数在不同位置不一样，有的位置是十二期里取最大，有的是三期里取最大，两者不可比。

我的应对办法是额外生成一张“有效计数”栅格，也就是用像元统计的“有效数据个数”选项再跑一遍，把它作为质量掩膜。凡是计数低于某个阈值（比如少于总期数的一半）的位置，在后续分析里标记为低置信度或者不参与。这一步多花两分钟，能避免后面一堆解释不清的异常值。

3. 实操：用像元统计工具做最大值合成

3.1 工具入口与参数逐项说明

工具位置在空间分析工具箱下的局部分析组里，名字叫“像元统计”或者“单元格统计”，不同版本的中文译名可能略有差异，认准英文 Cell Statistics 就不会错。前提是空间分析扩展已经启用，否则工具会灰掉或者运行报错。

参数一共就几项，但每一项都有讲究：

统计类型的下拉里有十来个选项：最小值、最大值、范围、总和、平均值、标准差、变异度、多数、少数、中位数、有效个数。最大值只是其中一个，但它和“有效个数”配合使用价值最高——前者出结果，后者出质量。

3.2 一步一停的操作记录

我的习惯是分四步走，每一步停下来检查一次，别一口气跑完再看结果。

第一步，把参与计算的图层全部加载进来，逐个确认坐标系、像元大小、范围，把不一致的先用投影或重采样处理掉。这一步可以在图层属性里看，也可以在目录里看元数据，我一般用批量查看的方式，一次比对所有图层的关键参数，比一个个点开快得多。

第二步，打开环境设置，把三个关键项配好。处理范围选基准图层，捕捉栅格选同一个基准图层，像元大小填基准图层的分辨率。这三项配完可以先关掉对话框，它们会随工具一起生效。

第三步，填工具参数。输入列表全选所有图层，统计类型选最大值，忽略 NoData 按前面的判断勾或不勾，输出路径写到工作目录下的结果文件夹，文件名带上日期和参数标识，比如max_2024_ignore_nodata.tif。文件名规范这件事看着啰嗦，但一个项目做十几版结果是常态，不写清楚后面根本分不清哪版是哪版。

第四步，运行完成后不要急着关，先看输出的范围、像元大小、波段数对不对，再看直方图分布是否落在输入数据的合理区间内。如果最大值出现在输入数据的极端上边界，或者整幅图只有一个值，基本可以判断出问题了。

3.3 图层一多，就该上批处理了

三期五期手动跑还行，几十期就没法看了。这时候有两个选择：模型构建器或者 Python 脚本。模型构建器适合不太写代码的人，把像元统计工具拖进去，把输入设成参数，用迭代器遍历文件夹里的所有栅格，一次性跑完。好处是可视化、能存成工具箱给同事复用。

脚本适合需要精细控制的情况，核心逻辑就是把文件列表收集起来，循环调用统计函数。大致骨架是这样：

import os import arcpy from arcpy.sa import CellStatistics arcpy.CheckOutExtension("Spatial") arcpy.env.workspace = r"D:\work\input" arcpy.env.overwriteOutput = True arcpy.env.snapRaster = r"D:\work\base\ref_grid.tif" arcpy.env.extent = r"D:\work\base\ref_grid.tif" arcpy.env.cellSize = r"D:\work\base\ref_grid.tif" tif_list = [os.path.join(arcpy.env.workspace, f) for f in os.listdir(arcpy.env.workspace) if f.lower().endswith(".tif")] out_max = CellStatistics(tif_list, "MAXIMUM", "DATA") out_max.save(r"D:\work\output\max_result.tif") out_cnt = CellStatistics(tif_list, "COUNT", "DATA") out_cnt.save(r"D:\work\output\valid_count.tif")

这段里有两个细节值得说。"DATA"这个参数就是“忽略 NoData”，与之对应的是"NODATA"，别写成布尔值，会报错。另外，环境变量里的 snapRaster、extent、cellSize 三个一定要在调用工具之前设好，设晚了不生效——这个顺序问题坑过我不止一次，脚本跑完发现结果网格和基准图层差了半个像元，回头查了半天才发现是环境变量写在调用之后了。

脚本运行前记得确认 Python 环境的扩展许可能正常检出，不同版本的模块命名可能略有差别，用之前先在交互环境里 import 试一下，比跑到一半报错强。

4. 结果校验与常见问题排查

4.1 结果对不对，四个动作就能判断

跑完不看结果就往下走的，早晚要返工。我固定做四件事。

范围与网格核对。看输出的范围、像元大小是否和基准图层完全一致，行列数是否匹配。这一步用图层属性就能看，三十秒的事。

数值区间核对。输出栅格的最小值不应该小于所有输入的最小值中的最小者，最大值同理。如果输出的最大值比任何一个输入的最大值都大，说明中间做过插值或者统计类型选错了。

逐点抽查。挑几个特征位置，用识别工具同时看输入各层和输出层的值，确认输出值确实等于输入中的最大者。我一般挑高值区、低值区、边缘 NoData 区各一个点，三类都过一遍，基本能覆盖大部分逻辑错误。

直方图形态比对。单个输入的直方图和最大值合成后的直方图，形态上应该有明显差别——最大值合成会把分布整体往右推，低值端被削掉。如果两张直方图几乎一样，八成是统计类型没改，还是平均值或者第一个输入被原样输出了。

4.2 常见问题速查

4.3 几个只有踩过才知道的细节

输入数量多的时候不要一次全加。二三十个图层一次性扔进去，内存占用会飙得很快，尤其分辨率细的时候。我的做法是分两轮：先把前一半合成一个中间结果，再和后半合，结果完全等价，内存压力小很多。这种“归并式”叠加对最终最大值没有任何影响，因为最大值运算满足结合律。

输出格式优先选 .tif。有些格式对 NoData 的支持不一致，输出再读进来时 NoData 可能变成 0，而 0 在最大值比较里会拉低结果，尤其当真实数据全是负值的时候（比如某些温差栅格）。用 .tif 并明确指定 NoData 值，能避开这类问题。

输入是整型还是浮点，会影响结果的可解释性。整型栅格输出的最大值也是整型，看着干净；浮点栅格容易出现 0.30000000000000004 这种值，做阈值判断时要留意。我一般不对浮点结果做等值比较，改成区间判断。

同一个位置最大值来自哪一期，其实也能查。像元统计旁边有个“最高位置”工具，输出的是最大值所在的输入序号。做时序分析时这个信息特别值钱——不仅知道最大是多少，还知道什么时候达到最大，一张图就能讲清楚物候节律。这一步很多人不知道，白白丢掉了信息。

提示：输出文件名里带上“max”“ignore_nodata”“count”这类关键词，半年后回来找结果时，你会感谢现在的自己。

5. 进阶：把最大值统计接进真实的分析流程

5.1 和其他统计量搭配，才是一套完整方案

单独一张最大值栅格，很多时候不够用。我一般会同时产出三张：最大值、有效计数、最大值出现的位置。三张放一起，才能回答“最大是多少、样本够不够、什么时候最大”这三个问题。

举个实际用法：做某个区域的生长季峰值分析，最大值栅格给出峰值强度，位置栅格给出峰值出现的时序序号，有效计数栅格用来剔除样本过少的边缘区域。三张叠在一起做分区统计，每个子区域能得到峰值均值、峰值时间的众数、有效样本比例，这套指标比单一的均值有说服力得多。

如果业务上还关心“最不利情况”，比如风险类分析里要看多个情景中的最坏值，那逻辑是一样的，只是把最大值换成最小值或者取分位数。工具还是那个工具，统计类型换一下。

5.2 分区最大值和逐像元最大值是两件事

有个需求经常被混在一起：算每个行政区内的最大值。注意这和逐像元最大值完全不同——前者输出的是一个表，每个区一行；后者输出的是一张栅格，每个像元一个值。

分区最大值要用分区统计工具，输入一个是分区矢量或分区栅格，另一个是要统计的值栅格，统计类型选最大值。这个工具在分区边界跨越像元的时候会有归属判断问题，默认是像元中心落在哪个区就算哪个区，边界附近的值可能和直觉有出入。如果对边界精度敏感，得先把分区转成栅格再做，或者提高分区栅格的精度。

还有个小细节：分区统计默认只统计完全落在分区内的像元，或者按像元中心判断，不同设置结果会有差异。同一份数据换一种设置，某些小区域的统计值可能整个变掉。做对比分析时，务必保证两次用的是同一套设置。

5.3 性能和数据量管理

栅格最大值合成是个 IO 密集加内存密集的活。几个能明显提速的习惯：把输入统一转成整型（如果业务允许），能省一半以上空间；给输入栅格建好金字塔，虽然对计算本身没帮助，但中间检查时速度快很多；输出到本地固态盘而不是网络位置，几十 GB 的中间结果走网络会慢到怀疑人生；尽量用工具而不是在计算器里做多层嵌套，标准工具的内部实现对分块处理优化得更好。

数据量特别大的时候，我倾向于先做空间分块，把研究区切成若干块分别跑，再拼接。拼接时注意块与块之间的接边要重叠一两个像元，或者用基准格网严格对齐，避免接缝处出现一条细缝或者重复计算。

另外提醒一句，做这类叠加分析，中间结果会迅速堆积。我给自己定的规矩是：每完成一个阶段，把确定不再需要的中间栅格删掉，只保留最终成果和关键中间产物，并在一个文本文件里记录每张图的生成参数。这个习惯救过我很多次——三个月后被问“这张图是怎么来的”，翻记录比翻记忆靠谱得多。

最后分享一个我常用的偷懒办法：把整套流程（环境设置 + 像元统计 + 有效计数 + 位置输出）做成一个模型工具，输入只留一个文件夹路径和一个基准图层，其余全封装好。以后再来同样的活儿，双击运行，去泡杯茶回来就完事了。前面花两小时搭的东西，后面能省几十个小时，这笔账怎么算都划算。