ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

ArcGIS Pro批量栅格重分类:Python脚本规则设计与工具封装实践

ArcGIS Pro批量栅格重分类:Python脚本规则设计与工具封装实践 1. 为什么要手动定义重分类规则批量之前的底层逻辑接触过栅格分析的人应该都不陌生重分类是栅格处理里几乎绕不开的一步。它本身不复杂就是把原有栅格数值按规则映射到新值上比如把DEM按高程区间分成5级把坡度图按临界阈值分成4级或者把土地覆盖类型重新编码。但真正做项目的时候麻烦不在“重分类”这一个动作上而在“需要重分类的栅格动不动就是几十个上百个而且分类规则还得完全一致”。ArcGIS Pro 里虽然有重分类工具面板操作一层一层点也不算麻烦麻烦的是重复。同一套断点值你今天手动填一次明天换一个测区再填一次后天换个年份的数据再填一次每次都要重新打开工具、重新输入规则、重新选输出路径。一旦数据量上来这套流程既耗时又容易手误输错一个小数点整批成果就废了。更现实的情况是项目验收要求可复现——你得说清楚你这套分类阈值依据是什么怎么应用的。手工点选的操作没有留痕事后补文档全是回忆录。所以我的选择是用 Python 写脚本把“手动定义规则”这件事做成半自动流程。熟悉 arcpy 的人会想到arcpy.sa.Reclassify但它默认的重映射方式其实并不完全适合手动规则真正的关键是自定义重映射表。这篇文章不会停在“调用一个工具”这种层面而是会把我自己在批量手动重分类里的完整实践拆开讲规则怎么组织、代码怎么写、参数怎么处理边界值、批量循环怎么避开常见坑、最后怎么封装成 Pro 里的脚本工具给同事用。如果你想在项目里真正落地而不是只看 API 文档这篇应该能帮你省下不少试错时间。1.1 重分类工具的核心机制在开始写代码之前先把arcpy.sa.Reclassify的工作机制理清楚。这个函数接收三个核心参数输入栅格、要重分类的字段、重映射表。重映射表定义了“原值或原值区间”到“新值”的对应关系。arcpy.sa模块里专门有RemapValue、RemapRange这组类用来构建重映射表。RemapValue适合处理离散类别的栅格比如土地利用类型编码一个原值对应一个新值。RemapRange适合处理连续值的栅格比如高程、坡度、NDVI需要定义区间范围。这两种类都要求范围边界是分段连续的简单说就是区间不能重叠且最好能覆盖整个有效值域。很多人第一次写重分类脚本最容易犯的错是觉得直接往Reclassify里塞一个列表就行其实它会报 TypeError。正确做法是先把边界值列表做成RemapRange对象再把对象作为重映射参数传进去。后面我会给出一个可以直接跑的代码骨架这里先说原理因为理解了这个机制批量循环时你才知道哪些参数需要变成变量哪些参数要保持固定。1.2 手动规则的价值可解释、可复用、可演替你可能会问ArcGIS 自带的重分类工具里有“分位数”“自然间断点”这些自动分类方法为什么还要手动定义规则我的理解是这样的自动分类适合数据探索阶段适合快速出一版结果看看趋势但真正进入生产流程、尤其是面向业务决策或成果交付时分类断点往往是有业务含义的。比如做洪水风险评价坡度分级阈值可能来自行业技术规范做生态敏感性分析植被覆盖度的分级阈值可能直接引用论文里已发表的参数。这些阈值不是数据自己长出来的是业务经验或者标准规定好的。手动重分类的核心价值就是把“这套业务规则”固化下来并且能批量复制到所有时空数据上。用 Python 写的好处更明显规则集中在代码里改一处全部重跑一遍成果自动更新。这正好引出一个判断到底什么场景值得写代码做批量手动重分类我的经验是只要满足两个条件就值得动手。第一输入栅格数量大于 5 个第二这套分类规则在可见的未来还会被重复使用。两个条件里满足任何一个写脚本的时间就能赚回来。2. 重分类映射表的规则设计从 Excel 到 Python 的迁移重分类规则的落地我习惯先做一张映射表。这张表是给人看的也是给代码用的底稿。你在 ArcGIS 里手动做重分类时面板上填的就是这个表在代码里它就是字典或列表。最常见的设计形态是一张两列的规则清单左边是原值区间右边是新值。假如我们要对一个研究区的 DEM 做高程分级原值区间含下界不含上界新值0 - 2001200 - 5002500 - 100031000 - 200042000 以上5这张表看起来简单但真到代码里就有讲究了。我下面分两个小节讲。2.1 断点设计的常见实践标准我接触过的栅格重分类项目里断点确定方式大体有三类来源。第一类是行业标准或技术规范这是最刚性的阈值清清楚楚写在条文里你不需要创造性发挥只需要严格对应。第二类是文献参考值适用于科研场景作者会在方法部分写明分类标准你复现时需要保持阈值一致。第三类是根据数据分布人为划定比如在满足可解释性的前提下让分类后的面积占比不至于过度失衡。不管来源是哪一类落到代码里之前都必须先确认区间边界是“左闭右开”还是“左开右闭”。ArcGIS 的重分类工具区间处理方式是带临界值语义的。如果你定义的是 100 到 200那 100 这个原值会进入这个区间200 则不进入。构建RemapRange时这个边界逻辑一定要和业务定义对齐。否则面积统计会莫名差出几十上百个栅格像元而且这种差异非常隐蔽。2.2 用字典管理映射关系我自己的项目习惯是先写一个映射字典把新值和区间列表关联起来再通过循环生成重映射对象。这样做的好处是规则一目了然改起来也方便。比如class_rules { 1: [0, 200], 2: [200, 500], 3: [500, 1000], 4: [1000, 2000], 5: [2000, float(inf)] }这里我用了float(inf)表示无上限避免因为数据最大值超出你预设末尾边界而整片变成 NoData。这是一个非常实用的细节后面会再展开。生成重映射表的代码可以抽成一个函数import arcpy from arcpy.sa import * def build_remap(rules): remap_ranges [] for new_val, (low, high) in rules.items(): remap_ranges.append([low, high, new_val]) return RemapRange(remap_ranges)这样规则和逻辑就分开了。规则调整的时候你只需要改字典不需要动循环和核心代码。这个习惯在我做多个年份、多个区域数据批处理时帮了大忙。3. 单栅格脚本跑通核心代码逐步拆解很多人在批量之前卡在单张栅格跑不通。所以我建议一定先把单张数据完整跑通再往循环里套。下面是我实际用过的一个最小可用版本我一行一行讲。3.1 环境与许可检查ArcGIS Pro 里跑 arcpy第一件一定要确认的事arcpy.sa模块能不能正常导入。这个模块对应 Spatial Analyst 扩展模块而且需要授权。如果没有在 Pro 界面里勾选启用脚本运行时会直接报错提示“模块不存在”或者“许可不可用”。最简单的检查方式是在代码开头强制设置import arcpy arcpy.CheckOutExtension(Spatial) from arcpy.sa import Reclassify, RemapRangeCheckOutExtension的返回值如果不等于字符串CheckedOut说明许可没有成功获取。严谨一点可以在正式流程前做判断。但在固定机器的项目环境里一般只要拿到过许可后面流程就稳定了。3.2 单栅格重分类的完整代码骨架假设现在要对D:/data/dem.tif做坡度分级规则就是前面那张表。完整流程是设置工作空间、构建重映射表、执行重分类、保存结果。import arcpy from arcpy.sa import Reclassify, RemapRange arcpy.CheckOutExtension(Spatial) arcpy.env.workspace rD:/data arcpy.env.overwriteOutput True in_raster rD:/data/dem.tif out_raster rD:/output/dem_reclass.tif rules { 1: [0, 200], 2: [200, 500], 3: [500, 1000], 4: [1000, 2000], 5: [2000, float(inf)] } remap_ranges [] for new_val, (low, high) in rules.items(): remap_ranges.append([low, high, new_val]) remap RemapRange(remap_ranges) reclassed Reclassify( in_raster, Value, remap, NODATA ) reclassed.save(out_raster) print(done:, out_raster)有几个参数我要单独说明。第二个参数Value指的是栅格值本身。如果你要重分类的栅格有属性表也可以填具体字段名。第三个参数是重映射对象。第四个参数NODATA的作用是定义范围内没有覆盖到的原值在输出里的取值方式。如果你填一个具体数值比如 99那么没人管的值都会变成 99如果填NODATA就变成 NoData。我在实际业务里建议优先用NODATA尤其是做分类统计时NoData 会被排除统计逻辑但一个异常值不会。万一某个栅格出现了你规则里没预料到的值比如 -9999 这种观测异常值它会直接进入 NoData 而非变成一个假分类这能让问题及时暴露。3.3 边界值处理的小技巧RemapRange的边界逻辑我再强调一遍因为项目里真的在这里翻过车。比如你定义了一条[100, 200, 1]语义是 100 到 199.999 这类小于 200 的值会变成 1200 本身不会。如果你还有另一条是[200, 300, 2]那 200 会落到第二个区间。所以规则表写 [0, 200]、[200, 500] 这种形式是没问题的相邻区间首尾相接。但如果你写的是 [0, 199]、[199, 500]那么在 199 即使是一个整型栅格里像 199.5 这种浮点情况也可能变得不可控。我的经验是边界尽量用统一规则并且把边界值放在规则表里让团队同行审阅比跑完发现统计对不上再回头排查要省事得多。4. 批量重分类的完整实现目录遍历与命名策略单张跑通之后接下来就进入正题批量。批量重分类说白了就是用一个循环把单张逻辑重复跑多遍。听起来简单但真正写好是不容易的。因为一旦数据量上到几十个你就不能一个文件一个文件手写输出了得有可靠的遍历逻辑和命名方案。4.1 批量输入的稳健遍历首先明确输入输出的组织方式。我习惯把待处理栅格统一放在一个输入目录下比如D:/data/dem/里面可能是多个年份的 DEM文件名带有年份信息格式可能是.tif也可能是一个 img、grid 文件夹栅格。区分格式会带来一个麻烦arcpy.ListRasters()是可以一次性列出来的但要注意通配符的写法。import os import arcpy input_dir rD:/data/dem arcpy.env.workspace input_dir raster_list arcpy.ListRasters(dem_*, TIF)这里的dem_*是通配符筛选TIF是格式过滤。我建议一定加上这两个参数理由很简单如果输入目录里混着一个不打算处理的某类辅助文件不加筛选直接全量遍历跑一半报错还要排查是谁混进来的。而且配合文件名规则你还可以只选中特定年份的数据比如dem_202*。另外要提到一个容易出问题的点如果目录里同时有tif和imgarcpy.ListRasters()返回的列表顺序并不保证和文件名排序一致。如果你在意输出顺序或者要和某个清单文件对应推荐先拿到列表后做一次sorted(raster_list)。这个细节在后续做多时相分析时非常重要。4.2 输出目录规划和命名策略输出目录千万别和输入目录混在一起。听起来是老生常谈但我确实见过有人把重分类结果直接保存到原目录结果下一次跑批时输入列表里混进了上一次的输出变成了重复重分类的雪球事故。我的方案是单独建一个reclass_output目录输出文件名统一在原名基础上加后缀比如dem_2000.tif变成dem_2000_reclass.tif。output_base rD:/data/reclass_output if not os.path.exists(output_base): os.makedirs(output_base) for ras in raster_list: base_name os.path.splitext(ras)[0] out_path os.path.join(output_base, base_name _reclass.tif) ...统一命名后缀的第二个好处是后续如果要按文件名排序、分组统计一眼就能区分原始数据和成果数据。尤其当一个项目跨了很多年成果文件越来越多时清晰的文件名能帮你省掉大量检索时间。还有一点是关于文件格式的。.tif是单文件栅格携带方便GRID 是文件夹栅格在批处理里路径操作会麻烦一些。我建议批量流程里统一输出为 TIFF格式过滤和路径解析都干净也不容易混进别的杂项文件。5. 实测中踩过的坑NoData、数据类型与性能代码写出来能跑跟跑得稳、跑得准之间中间隔着大量实际数据验证。下面是我在多个项目里反复踩过、最后沉淀下来的几个坑分享出来供参考。5.1 NoData 究竟是哪一个值不同来源的栅格数据对 NoData 的定义很不统一。有的用 -9999有的用 -3.402823e38有的只在栅格属性里标记但没有实际占用特定值。这就导致你定义重分类规则时很难提前知道是否有值会被漏掉。解决思路有两种。第一种是在重分类之前先用arcpy.Describe或arcpy.Raster的noDataValue属性查看。第二种更稳直接在重分类时把未覆盖值交给NODATA处理如前面代码所示。这样就算规则漏掉了一部分也不会被硬塞到一个不合理的类里后续检查统计结果时很容易发现。5.2 整型浮点型混用时的行为差异如果你的原始栅格是整型那么重分类后输出一般是整型。如果原始栅格是浮点型输出也往往是浮点型。但有些规则定义时用[low, high, new_val]new_val 是整数像 1、2、3如果原始是浮点型输出栅格的实际值可能是 1.0 而不是 1。这通常不影响分析但如果你要把它用于后续的分类统计或字段关联要注意类型问题。更麻烦的是有些工具生成的栅格虽然表面上值看起来是 1、2、3实际像素类型却是浮点导致你按整型字段去关联属性表时报错。所以在批量跑完后加一个小步骤用arcpy.Describe(out_path).pixelType抽查输出类型是很有必要的习惯。5.3 性能问题为什么批量跑得越来越慢栅格数据通常体积不小一个地区的高精度 DEM 动辄几百 MB。批量重分类时如果每跑一个文件都在读盘、写盘、刷新渲染性能会越来越差尤其当成果数据积累到一定规模磁盘碎片和缓存占用都会拖慢速度。我的性能优化手段集中在三个地方。第一设置arcpy.env.scratchWorkspace指向一个独立的临时目录避免中间数据和成果数据互相干扰。第二在代码开头设置arcpy.env.parallelProcessingFactor 75%让工具尽量利用多核。第三如果机器内存够大可以考虑在循环里把中间变量定期清空用del和gc.collect()配合避免栅格对象堆积占用内存。import gc for ras in raster_list: ... del reclassed gc.collect()这个小技巧在跑上百个文件时效果立竿见影内存占用能稳定在一个较低水平。6. 将脚本封装成 ArcGIS Pro 工具让批量流程变成可配置的工程资产单人单机跑脚本很方便但如果你在一个团队里工作或者这个流程需要交给不太熟悉 Python 的同事使用我建议把脚本封装成 ArcGIS Pro 的自定义脚本工具。这样别人打开工具面板填几个参数点确定就能批量处理不必碰代码。6.1 参数面板设计让规则变成输入参数封装工具的核心是定义参数。拿我这个重分类脚本来说有三个参数是必须暴露出来的输入栅格、输出目录、分类规则。前两个简单直接换成工具参数即可。分类规则参数稍微难一点因为规则是二维结构有区间和新值。我的方案是做一个多值字符串参数每一行用分号分隔每个字段之间用逗号分隔。比如0,200,1;200,500,2;500,1000,3;1000,2000,4;2000,inf,5然后在脚本里解析这个字符串重新构建规则字典。这样做的好处是同事不需要懂 Python也能在工具面板里清晰看到规则。坏处是格式必须约定好写错一个分隔符整个脚本就报错了。所以我在解析时做了一个辅助函数能检查每段是否有三个值少了就报出明确定位信息def parse_rules(rule_string): rules {} parts rule_string.split(;) for i, part in enumerate(parts): items part.split(,) if len(items) ! 3: raise ValueError(f第{i1}段规则格式错误: {part}) low float(items[0]) high float(items[1]) if items[1].lower() ! inf else float(inf) new_val int(items[2]) rules[new_val] [low, high] return rules这个函数让工具面板成了一个低门槛入口规则还是那个人能看懂和修改的规则。6.2 封装脚本工具的注意事项如果你在 ArcGIS Pro 的“目录”里右键工具箱选“添加脚本”会进入脚本工具设置。有几个容易被忽略的设置点我说一下第一参数的数据类型必须选择“栅格图层”或“工作空间”不能选字符串不然传参路径容易有歧义第二在“语法”里要勾选“从脚本工具获取参数”也就是脚本里用arcpy.GetParameterAsText(0)来取值而不是硬编码路径第三工具的“运行”模式建议选择“在后台运行”因为批量重分类耗时较长前台运行容易把 Pro 界面锁死。我封装完之后同事只需要在面板里选好输入文件夹、输出文件夹、粘贴规则字符串点运行。整个流程不需要任何人接触代码。考虑到批量重分类可能是整个栅格分析链里的一个环节这个自定义工具还可以继续作为子模块嵌套进更大的 ModelBuilder 流程里。7. 从批量重分类到完整工作流的扩展跑完批量重分类工作一般还没结束。重分类往往只是分析链路的中间步骤之后还要统计各分类的面积、生成专题图、甚至做多期数据的变化检测。我不会在这里展开全套流程但我想讲两个我实际用过的扩展方向给你一点参考。第一个方向是批量统计分类面积。重分类结果栅格往往有属性表表里每一行对应一个类别Count字段是像元个数乘以像元面积就是类别面积。可以在重分类后的下一个循环里直接用arcpy.Raster转属性表或者直接访问.attributeTable。不过要注意属性表要用arcpy.sa.ZonalGeometryAsTable或者简单地用arcpy.management.BuildRasterAttributeTable处理后才稳定。第二个方向是多期重分类结果的变化分析。比如你有 2005 年、2010 年、2015 年三期地表覆盖重分类成果可以直接用栅格计算器做差值快速找出类别发生变化的位置。这里重分类的“手动规则一致性”就体现出来了三期数据采用的断点完全一致输出类别编码也完全一致差值图才能直接解释否则三套规则的分类编号互相错位整个变化分析根本无法开展。我在实际项目里把重分类脚本放到一个更大的工作流里以后真正体会到一件事写代码不是为了让事情“自动化”而是为了让决策规则“留住”。你可以不断迭代规则可以对比新旧规则对成果的影响可以在项目交接时把规则和数据一起交付。这套代码时至今日还在被我复用不同项目的规则内容不同但骨架几乎没有变过。如果你也经常和大量栅格数据打交道希望这篇内容能帮你跨过批量重分类里的一些隐形门槛。遇到问题欢迎多交流。
返回列表