摘要
矿业权核实坐标因其直观可读性广泛应用于矿产压覆评估和地质勘查项目查重等工作。然而,矿业权核实坐标存储在字符串中,空间化时需手动拆解为经纬度列,难度较大。为提高效率,本文基于ArcGIS Pro Python提出了矿业权核实坐标与 ArcGIS要素类批量互转的思路,并用 Python实现,封装为脚本工具。通过生成2000个随机坐标要素,验证了工具的转换精度和可靠性。文章还探讨了转换过程和Python开发中的关键问题,为矿业权核实工作提供了参考和依据。
关键词
Abstract
The coordinates of mining rights verification are widely used in mineral overburden assessment and geological exploration project duplication checking due to their intuitive readability. However, the coordinates of mining rights verification are stored in strings, and they need to be manually disassembled into longitude and latitude columns when spatialized, which is difficult. To improve efficiency, this paper proposes the idea of batch conversion between mining rights verification coordinates and ArcGIS feature classes based on ArcGIS Pro Python, and implements it in Python and encapsulates it as a script tool. The conversion accuracy and reliability of the tool were verified by generating 2000 random coordinate features. The article also discusses the key issues in the conversion process and Python development, providing a reference and basis for mining rights verification.
0 引言
随着市场经济发展,因建设项目压覆矿产资源涉及的补偿事务越来越多,矿业权核实坐标广泛的应用于建设项目矿产压覆评估、地质勘查项目查重等工作中(尚晓鹏等,2017;董少波等,2021;石文君等,2024)。矿业权坐标是常用的坐标传递形式(宋子东,2021;和吉豫,2023),相关学者针对此坐标形式做了许多研究,朱树叶和杨波(2017)在其研究中详细介绍了宁夏回族自治区矿业权数据整合及坐标转换的项目,展示了如何通过技术手段实现数据的统一管理和更新,从而提升矿业权管理效率;杨文府(2012)聚焦于山西省矿业权实地核查工作,讨论了如何通过建立属性数据库和空间数据库来加强矿政管理,并提供了数据成果综合应用的实例; 侯莉(2020)的研究提出了基于GIS二次开发组件和空间数据存储技术的矿业权坐标数据GIS管理系统的设计与实现方法,为矿业权坐标数据的一体化、智能化、可视化管理提供了解决方案;丁晓晓 (2020)通过实例验证,探讨了不同坐标系之间的转换方法,为矿业权坐标转换工作提供了参考。
相比经纬度列形式存储的结构化坐标点数据,矿业权核实坐标可以直接判读出矿业权的部件个数与挖洞情况,在直观性方面具有显著优势(汪凯明,2023)。但将其空间化为要素类较为困难,需要手动拆解为经度和纬度列形式,然后使用 GIS 类软件空间化(王利东,2021),拆解过程中容易出现错误,导致空间化偏差或失败。本文基于 ArcGIS Pro Python,提出一种矿业权核实坐标与 ArcGIS 要素类批量互转思路,并使用 Python 实现,为了提高代码使用率,将Python代码封装为ArcGIS Pro脚本工具,同时探讨转换过程和 Python 语言开发中遇到的基础层面问题,为矿业权核实工作提供一定的参考和依据。
1 基本概念
1.1 矿业权核实坐标
矿业权核实坐标分为部件数和部件实体两部分,其中部件实体分为点数和坐标对两部分(图1)。各部分之间分割符均为英文逗号(,),经纬度的坐标格式为(DDD.MMSSssssss),度与分之间的点仅做分割作用,非小数点意义,这与十进制度坐标极易混淆。秒与秒的小数之间没有小数点做分割,是按 2 位截取计算,小数部分可以自定义长度,通常在0-6 位之间选择,部件结尾有填充区(0,0,0)或挖洞区 (-1,0,0)两种形式。
图1矿业权核实坐标基本结构
例如:1,3,121.5001,42.3102,121.5546,42.3102,121.5546,42.2632,0,0,0
示例坐标为最简坐标形态,仅具备1个部件,且部件只有3个点。
1.2 ArcGIS Pro Python
Python 是一种解释型语言,可以在源代码执行时直接翻译成机器码,这使得 Python代码易于调试和修改。在 ArcGIS Desktop 中,使用 Python2.7 版本 (张宇冉等,2021;赵军鹏和刘军,2021;孙亚珍等, 2023),对中文支持较差,对GIS开发人员不够友好,往往在解决业务问题的同时还要考虑到中文兼容性,影响开发体验。从 ArcGIS Pro 版本开始,Python 使用 3.9及以上版本(孙跃龙和鞠立华,2024),完美解决中文问题,并提供多种新库支持,为GIS开发提供便捷。ArcGIS Pro Python具备如下特点和功能:
(1)地理数据处理和分析:用户可以使用 Arc‐ GIS Pro Python编写脚本来读取、写入和操作地理数据。这包括筛选、裁剪、合并等操作,使用户能够灵活地处理地理数据。
(2)地理空间分析和建模:ArcGIS Pro提供了强大的地理空间分析工具,用户可以通过 Python脚本来调用这些工具,完成复杂的地理空间分析和建模任务。例如,缓冲区分析、网络分析、地形分析等都可以通过Python脚本来实现。
(3)自动化地理数据处理:通过编写 Python 脚本,用户可以自动化地处理地理数据,提高工作效率。这包括批量导入、导出、转换地理数据等操作。
(4)与其他软件集成:Python是一种通用的编程语言,因此在 ArcGIS Pro 中可以方便地与其他软件进行集成。这使得用户能够更灵活地处理和分析地理数据,并将其与其他类型的数据进行结合。
2 实现思路
2.1 矿业权核实坐标转要素类
在ArcGIS Pro中,常用的是点、线、面、多面体等 4种要素类,矿业权核实坐标表示的是封闭区域,因此需要将其转换为面要素类。转换要素类过程中,笔者认为首先应该解决单个坐标串转要素的情况,然后再实现批量转换。
1)单个矿业权核实坐标转要素类
具体分为参数输入与部件解析、坐标点格式转换、几何体构建和要素输出等步骤。
(1)参数输入与部件解析
输入参数即为矿业权核实坐标字符串,截取字符串第一个逗号后半部分(cut_partVal_str),使用代码1,将其分割为部件列表,其中re为Python正则对象。
为正则表达式,同时考虑到了填充区结尾、挖洞区结尾、结尾逗号的不确定性、空字符串和“正则贪婪”情况。
为正则表达式,同时考虑到了填充区结尾、挖洞区结尾、结尾逗号的不确定性、空字符串和“正则贪婪”情况。
(2)坐标点格式转换
由于ArcGIS中保存的坐标为十进制度格式,故需要将每个点的坐标由 DDD.MMSSssssss 转为十进制度,转换思路如图2,先将小数点后移四位,再按照固定宽度分割字符串,最后使用公式 1 得到十进制度。
图2DDD.MMSSssssss转为十进制度
(1)
式(1)中,D10为十进制度,DDD为度整数部分, MM为分整数部分,SS.ss为秒部分。
(3)几何体构建
ArcGIS 要素中,存放的是几何实体,使用代码 2,将坐标矩阵转换为ArcGIS的面对象(Polygon),其中坐标矩阵 feature_info 是三维数组,第一维代表几何,第二维代表部件,第三维代表坐标。arcpy.Poly‐ gon 是构建面对象的核心方法,其中 WKID 为 4490 (ArcGIS 官网,2024a)的投影方式对应国家 2000 大地坐标系(GCS_China_Geodetic_Coordinate_System_2000)。
(4)要素输出
将内存中的面几何要素列表转换为数据库存储的要素类(图3),使用拷贝要素类功能即可实现 (代码3)。其中fc为要素类保存路径。
2)批量矿业权核实坐标转要素类
本次使用ArcGIS Pro的模型构建器来将单一任务转为批量任务,ArcGIS Pro 模型构建器是一种可视化编程语言,用于构建地理处理工作流。使用模型构建器,可以便捷地制作工作流,降低开发难度。
(2)批量模型构建
批量执行的整体思路是遍历所有矿业权坐标串,以此调用单个转换脚本,并将转换结果使用收集值功能暂存在内存中,最终使用合并功能导出结果,模型详见图4。
①本次将输入参数设置为表类型,这样既可以从 Excel 等第三方表格中读取数据,也可以从 Arc‐ GIS要素类或表格中获取数据。
②在遍历操作后增加“计算值”功能,该功能没有业务功效,代码块中编写如下代码:return "正在进行第"+str(a)+"行",可以提示正在进行转换的数据行数,在数据量较大的情况下,出现问题便于定位问题位置。
③获取字段值功能是将遍历到的行中坐标值取出,坐标字段列作为模型第二个参数,由用户选择坐标列。
④单个字符串转要素类脚本中,输出路径设置为 memory\YS%值%,是将此临时变量暂存在内存中,%值%为 ArcGIS Pro 提供的行内变量,用来获取该模型中正在进行的行ID。
⑤收集值配合合并工具,将内存中的要素类合并为一个导出,值得一提的是,按照图4的模型写法,每循环一次都会执行一遍合并功能,虽然对最终结果无影响,但会降低效率。若要实现最后一次循环结束后再执行合并,需要将模型拆分为主子模型,本文不再具体探讨。经测试,在要素超过 2000 个时,目前模型会有明显速度降低现象。
图3单个转换方法封装
图4批量矿业权核实坐标转要素类模型构建
2.2 要素类转矿业权核实坐标
将 ArcGIS Pro 要素类(GDB、MDB 或 shp 格式) 批量转换为矿业权核实坐标格式,具体可划分为要素投影、游标遍历、部件解析、坐标转换、字符串输出等步骤。
(1)要素投影
由于矿业权核实坐标格式为国家 2000 经纬度投影,故首先应将要素类投影为该类型坐标系。
(2)游标遍历
使用 Python 读取 ArcGIS 要素类的方式有很多,最灵活高效的是使用 SearchCursor游标读取(ArcGIS 官网,2024b),如代码4,fc变量为要素类路径,该游标可直接读取要素的GeoJson字符串,从中提取拐点坐标,Json中拐点坐标位于rings属性内,该属性为多维数组,其格式与代码2中feature_info变量格式一致。
(3)部件解析
feature_info 数组为三维数组,其中第一维度即为要素、第二维度为部件,第三维度是拐点。考虑到部件分为填充部件和挖洞部件,不同部件结尾字符串不一样,因此需要对所有部件进行类型判断。首先应该确定一个主部件,主部件往往面积最大,其余部件与主部件做包含判断,当包含在主部件内时,为挖洞部件,不包含在主部件内时,为填充部件。
(4)坐标转换
ArcGIS Pro 要素类中,拐点坐标使用十进制度形式存储,此时需要转换为矿业权核实坐标形式,转换思路如图5,先将度小数部分乘以 60,得到分,再将分小数部分乘以 60,得到秒。但需要注意的是,当分和秒小于10时,需要右对齐左补零占位,否则会出现坐标错乱。秒小数部分需左对齐向右补零,满足用户保留小数位数个数。最终将各部分坐标拼接,得到矿业权核实坐标格式。
(5)字符串输出
在每次要素遍历结束时使用 arcpy.AddMessage()方法输出矿业权核实坐标字符串,使用 AddMes‐ sage 方法的好处是无需编写结果保存方面的代码,直接将矿业权核实坐标输出至消息台,用户可在工具执行消息中复制转换结果。
图5十进制度格式转矿业权核实坐标格式
3 关键问题探讨
在实际使用过程中,会遇到各种各样的问题,笔者将最重要且会影响坐标精度的问题列举如下。
3.1 float精度受损问题
在计算机内部,浮点数(float)是使用二进制来表示的。但是,有些十进制小数无法精确地转换为二进制小数(方星星和吕永强,2019;王兆华, 2020)。例如,0.1是一个无限循环的二进制小数,由于计算机在内部使用有限的二进制位来表示这些数字,所以它们只能近似地表示这些十进制小数,导致精度丢失。在 Python 或其他开发语言中 0.3-0.1 并不等于 0.2,而是 0.19999999999999998,如图6。这种情况被称为float精度受损。
图6Python中0.3-0.1结果
float精度受损问题是计算机开发中常常遇到的问题,本次大量的用到了 Python 代码,且涉及大量浮点小数运算后转为字符串做截断的操作,难免会因为精度受损而造成坐标大范围偏移的问题,例如正确应该得到的是 00 秒,使用字符串截断的方式,在某种巧合下会得到 99 秒,造成坐标发生大量偏移。
针对此类问题,推荐使用Decimal类型而非float 类型做浮点数的运算,Decimal类型在各编程语言中是金钱常用计算类型,精度比较可靠,Decimal 使用时,要包裹字符串类型,不能直接包裹 float,否则没有效果,如图7a,包裹float下,120.3-120不等于0.3,而图7b中则等于0.3。
图7Decimal中包裹float(a)与包裹字符串(b)使用
3.2 补零问题
前文已述,在十进制度转度分秒时,分和秒整数部分为个位数时,应该在左方向补零,否则坐标会由于错位导致偏差;秒的小数部分也应该右方向补零,满足用户的小数长度。
除此之外。还有一种情况未提及,秒小数部分左方向补零问题,例如,该坐标为十进制度形式120.5983427,转换为秒时为 54.03372 秒,假设用户需要秒级小数保留 3 位,若使用字符串分割整数和小数部分的方法,再对小数部分做四舍五入操作,得到的秒小数部分将变为 337,会丢失左侧 0,而非正确值 034,从而造成坐标偏移。此类问题的解决方式应该尽量使用减法减去整数部分,而非字符串分割,同时在使用减法时,需要用Decimal类型。
3.3 科学计数问题
在 Python 中,当小数位数过多时,把该数字转为字符串会变成科学计数法,例如:0.00001转为 str 类型,会变成'1e-05'。
本次工作中,假设用户要求秒小数点保留6位,度=120,分=00,秒=00,秒小数=0.000001 时,使用 Python 字符串拼接形成的坐标会变为 120.00001e-05,而非120.0000000001,从而造成坐标错误。
此类问题的解决方案是,对于精度较高的小数转字符串,使用"{:f}".format(xxx)方式,而非 str(xxx),这样可以避免出现科学计数法。
3.4 挖洞部件问题
前文中阐述了简单的与主子部件识别算法与挖洞部件认定方式,但还有一种形式会造成挖洞部件认定错误,主部件为右侧面积较大部件,遍历两个子部件与主部件做比较,均为非包含情况,则会认定子部件1和2均为填充部件,与事实不符。
此类问题的解决方案是,需要将所有部件均假设为主部件,做一遍包含遍历,一旦有一个部件存在包含关系,那该部件应该属于挖洞部件,反之亦然。
笔者本次并没有采用该方式,主要原因是实现过于繁琐,且此类情况极少出现。即使出现导致子部件 2 结尾标识符号反向,但对导出的矿业权核实坐标字符串使用并无影响。在矿业权核实坐标字符串转为要素类使用的过程中,不需要告知 GIS 软件部件的性质(代码 2),仅需要构造好要素数组即可,GIS 软件会自动根据叠复情况判断是否为挖洞部件。
3.5 秒小数保留位数
若矿业权核实坐标中秒不保留小数,则与实际最大可能偏差 30 m 左右,随着小数保留位数的增多,由保留小数造成的坐标误差也随之减少,经测试,秒级小数保留6位时,矿业权核实坐标与要素类互转将不存在偏差。这是由于ArcGIS中,经纬度投影的要素类坐标精度为 0.00000001度,转换为秒为 0.000036秒,此时秒级小数位数刚好6位。
4 正确性核验
本次使用 ArcGIS Pro 随机生成了坐标不等、大小各异、部件复杂的2000个要素并存放至一个要素类中,通过要素类转矿业权核实坐标工具将其全部转换为字符串坐标后,再使用矿业权核实坐标转要素类工具,将字符串坐标批量转换为要素类。
将前后两个要素类放在同一个要素数据集中,建立拓扑规则(a要素必须被b要素覆盖),做高精度 (9位小数容差精度)拓扑检验,验证拓扑后结果0错误;0异常(图8),说明两个要素类位置完全一样,没有出现坐标偏差。
图8验证拓扑结果
5 结论
(1)基于 ArcGIS Pro Python 可以实现快速批量的矿业权核实坐标与ArcGIS要素类批量互转,减少人工重复,提高工作效率,本文提供了核心思路和算法代码。
(2)本文为日常矿业权核查工作提供了一定的参考,并促进该领域的技术交流和进步。
(3)本文所遇到的关键问题及解决方案,不仅在此项目中可用,在其他 Python项目中仍然有重要的借鉴意义。
(4)使用拓扑验证功能可快速的检查要素类正确性。