郭 霄 薛 睿 武慧琳 荀张媛

(1. 61363部队, 陕西 西安 710000; 2. 自然资源部第一大地测量队, 陕西 西安 710054;3. 自然资源部测绘标准化研究所, 陕西 西安 710054; 4. 中煤航测遥感集团公司, 陕西 西安 710199)

0 引言

标准经纬度分幅地形图接图表作为一种项目管理不可缺少的重要数据,在测绘生产领域有着十分重要的作用。当前地方数据格式管理不一,实际工作中可能遇见数据分别采用新旧编号规则管理的图幅,导致在数据汇交管理时,需重新修改接图表数据属性或重新生成接图表,现有的制作接图表方法操作复杂效率低。因此本文提出一种快速生成接图表的算法就具有十分重要的意义。

目前主流制作接图表方法依据目的分为两类:一是通过经纬度范围生成接图表,主要通过ArcGIS中ArcToolbox提供的“创建渔网”工具,生成接图表的骨架,使用者需要提供渔网的起始枢轴点,要生成对应比例尺的经差与纬差等参数共11项[1]。该方法处理不够自动化,图号等属性需要后续通过字段计算器添加,费时费力,对操作者专业水平要求较高,容易出现操作失误导致的渔网位置错位等问题。二是通过图号创建接图表,主要通过Arcpy实现,分为3步：①根据文件名获取图幅号列表;②按照图幅名计算坐标;③依照计算出的坐标生成接图表多边形,主要用到的函数有arcpy.Point()、arcpy. Array()、arcpy. Polygon()[2]。此方法一定程度的提高了生成接图表的效率,但后续依旧需要通过字段计算器添加图号信息,不能一步完成。

为更好地解决以上问题,本文提出了一种新的算法模型,该模型在部分以往的计算基础上,借鉴了图像处理深度优先搜索(depth first search,DFS)算法逻辑,减少了操作者干预接图表骨架的生成,提高了准确度,将图号信息在算法中直接添加,避免了后续操作,并将该算法与传统方法进行对比实验进行验证。

1 接图表生成算法

1.1 地形图经纬度分幅与编号

1.1.1地形图经纬度分幅

地形图经纬度分幅指地图的图廓由经纬线构成,我国的基本比例尺就是以经纬线分幅制作的[3]。它们是以1∶1 000 000地图为基础,按规定的经差和纬差划分图幅,行列数和图幅数成简单的倍数关系[4]。

1.1.2新旧图幅编号

新旧图幅编号规则都是在标准经纬度分幅的基础上设计,新图幅号由10位码组成,以1∶50 000图号为例,J50E001002中J代表该图幅所属1∶1 000 000图幅行号数字码、50代表其所属1∶1 000 000图幅列号数字码、E代表所属比例尺代码、001与002分表代表图幅行号与列号数字码[5]。旧图幅编号系统按图1编码。

所有比例尺都可实现新旧图号转换,以1∶50 000比例尺为例,其对应公式如表1所示,其中H代表新图幅中的行号,L代表新图幅中的列号,X代表旧图号中的地图代码值,[]代表取整,()代表取余。此为灵活生成不同编码规则的图号接图表奠定了算法基础。

图1 我国旧式地形图编号系统

表1 新旧图号转换

1.2 整体计算流程

本文提出算法将根据经纬度范围生成接图表与根据图号生成接图表结合,整体流程如图2所示。图号文件可根据需求输出为新编码规则图号与旧编码规则图号txt文件。

1.3 搜索计算图号

1.3.1深度优先搜索(DFS)[6]

DFS是计算机术语,是一种在开发爬虫早期使用较多的方法,是搜索算法的一种。它的目的是要达到被搜索结构的叶结点。

深度优先遍历图的算法是,假定给定图G的初始状态是所有顶点均未被访问过,在G中任选一个顶点i作为遍历的初始点,则深度优先搜索递归调用步骤[7-8]。①访问搜索到的未被访问的邻接点;②将此顶点标记为已访问节点;③搜索该顶点的未被访问的邻接点,若该邻接点存在,则从此邻接点开始进行同样的访问和搜索,反复进行直到所有节点都被访问为止[9]。

图2 整体算法流程

1.3.2计算图号

本文将经纬度范围围成的矩形四角顶点类比作图G的顶点i,这里称为基点,以基点为根结点搜索与其相邻的图幅,相邻图幅作为其子节点[10-11]。

当输入数据为经纬度范围时,根据下式计算四角其对应的新式图号编号,a为1∶1 000 000图幅所在经度带的字符所对应的数字码;b为1∶1 000 000图幅所在纬度带的字符所对应的数字码[12];c为所求比例尺地形图在1∶1 000 000地形图编号后的行号;d为所求比例尺地形图在1∶1 000 000地形图编号后的列号;λ为某点的经度;φ为某点的纬度;Δλ为所求比例尺的经差;Δφ为所求比例尺的纬差[13-14]。

以矩形4个交点图号为控制搜索生成边框图号,当搜索结束后,计算四角点图号内对角图号,再以内对角图号为基准再次进行搜索。当搜索进行到最后一个循环时会有6种情况,分别为(a)田字分布(b)竖向一字分布(c)横向一字分布(d)单幅分布(e)竖向两行分布(f)横向两行分布。

针对(a)情况,当四角图号的对角图号为其本身,则跳出搜索;针对(b)、(c)情况以两边图幅为基准向下或向右搜索图号至完毕;(d)中当四角图号的对角图号为同一图号时跳出循环,结束搜索;(e)(f)判断向下或向右搜索是否右图幅如果没有,结束搜索。

1.4 通过图号生成接图表

在已知新图号的前提下,通过(5)(6)式可计算出该图幅在西南图廓点的经纬度坐标[4]。

通过该比例尺所属经差与纬差可计算出其他3个角点的经纬度坐标,获取的XY坐标,写出文件时将四角点坐标与其对应图号一起写出,为后续直接在接图表面数据中增加图号信息做准备。根据实际需求转换图号编码实现需要图号的直接输入。

将XY坐标转为点要素并通过spatRef = arcpy.SpatialReference(4326)设置其空间参考,再将点要素生成线,线转面时将之前写出的图号写入面要素属性中。完成接图表制作,具体流程如图3所示。

图3 通过图号计算接图表算法

2 实验与分析

2.1 算法的实现

2.1.1计算对角图号实现

搜索计算四角图号由4个主要函数组成,分别是ldown(string number)计算基准图幅左下角图号;lup(string number)计算基准图幅左上角图号;rdown(string number)计算基准图幅右下角图号;rup(string number)计算基准图幅右上角图号。每一个计算方法都有4种情况：(1)a变b不变；(2)a不变b变；(3)a、b同时变化；(4)a、b都不变,最后返回计算图号。

2.1.2循环获取边框图号

循环获取图号主要借助toright(string num, string destination);todown(string num,string destination),采用while(true)死循环,当搜索图号num与目标图号destination相等,或num图号排在destination下方和左方时循环跳出,搜索结束。

最后一个循环借助hashset结构,在预先存储了搜索基准图号的前提下,hashset的长度为4,hashset主要特点为无序的不重复,因此借助hashset可以了解最后一个循环属于上文6种情况中的哪一种,利用while(true)死循环,当为(a)情况时添加其对角图号后hashset长度为仍为4,跳出循环;(b)(c)情况时hashset长度为6,借助todown(string num,string destination)、toright(string num, string destination)完成搜索后跳出;(d)情况时长度为5,break跳出循环。

所有图号最终都存储在hashset中通过foreach循环遍历借助stream类写出。通过python实现角点图号与边框图号函数编写,并将脚本插入Arctoolbox中,方便加入modelbuilder使用。

2.1.3图号生成接图表实现

借助modelbuilder,实现图3流程如图4所示。其中“计算xy坐标”与“xy转为点要素”工具为插入的脚本工具。

图4 Modelbuilder实现算法

2.2 与传统方法对比分析

2.2.1参数对比

本文计算经纬度范围内的图号,只需输入上下纬度与左右经度值,与保存路径。也可以将该脚本拖入的Model中,一步实现生成。

依据图号生成接图表界面参数包括存储图号文件的路径,判断是否为一条线上的点的字段name,以及输出路径。

而传统方法所需参数繁杂,后续还需要在属性表中创建字段,使用字段计算器,计算几何还有投影设置等设置接图表图号、投影属性[15]。

2.2.2实验结果对比

以生成1∶50 000数据接图表为例,假定使用传统方法熟练的情况下,与本文方法对比如图5所示，其中粗线为本文方法,细线为传统方法。

从图5可以看出传统方法人工操作执行时间受操作者的影响较大,且基本运行时间是本文方法的1～2倍。

图5 整体对比分析

3 结束语

本文借鉴深度优先搜索算法逻辑,类比接图表格网为栅格数据的像素格网,以接图表四顶点开始搜索经纬度范围内的所有图号,并在通过图号生成接图表的过程中将图号属性写入接图表矢量数据中,与现有制作接图表的方法相比减少了处理步骤与人工操作,克服了因人工误操作导致的误差问题。通过多次对比实验分析证明,该方法大大提高了接图表生成的效率,有效地降低了制作的复杂度。但该方法还存在识别图号格式不够灵活的限制,例如在图号中如果出现其他连接符需更改为标准图号格式算法才可使用,未来将对算法进行改进使其更加灵活适应图号的不同格式,扩大其适用领域。