求解制图区域的地图图幅编号的算法研究

吴曜宏，乔俊军，胡冯伟,2

(1.武汉大学测绘学院，湖北 武汉 430079；2.浙江省测绘科学技术研究院，浙江 杭州 310012)

0 引 言

地图是国民经济建设、国防建设和科学研究的重要基础地理信息资料[1]，随着计算机科技的进步，计算机地图制图逐渐取代手工地图制图，成为地图生产的主要形式，国家基础地理信息数据库的建设更是极大地提高了地图制图的效率与精度。当我们进行地图制图时，仅需要根据制图区域的地理位置和形状大小，准确选择该区域的基础地理信息数据便可开展地图制图工作，因此，对基础地理信息数据范围的确定及编号的可视化查询就显得十分重要。

为了方便基础地理信息数据的管理查询和准确使用，我国基于地理格网制定了一套新的GBT13989-2012《国家基本比例尺地形图分幅和编号》[2]国家标准，该标准同样适用国家基础地理信息数据库的分幅和编号。为了实现地形图编号查询的自动化，国内许多学者做了相关研究，并取得了一定的研究成果[3]，但是，功能比较单一，有的没有实现点、线、面等多维度查询；有的查询结果是利用地理坐标系平面投影化后的环境生成的，在泛区域的情况存在图幅查询空缺。目前，还没有根据各种投影后的制图区域确定基础地理信息数据范围和查询基础地理信息数据幅编号等方面的研究，基础地理信息数据是否完全覆盖制图区域，直接影响着地图制图工作的进一步开展。

针对上述问题，本文基于投影后的制图区域，从点、线、面三个方面提出了归原法、斜率分段-同异侧判别法和投影反算-图廓内外判别法，这些方法可准确映射各种投影后制图区域所对应的基础地理信息数据范围，实现了基础地理信息数据编号的可视化查询，避免了图号缺失和冗余，为基础地理信息数据库驱动地图制图奠定了数据基础。

1 理论基础

1.1 基础地理信息数据范围的确定

地图是基础地理信息数据投影后的产物，由于投影方式各异，所以，经纬线形状及间隔也各不相同，这就造成了同一区域、同一比例尺且图廓大小相同的情况下，所对应的基础地理信息数据范围各不相同。下面以圆柱投影和兰勃特投影为例，不同投影下同一制图区域所对应的不同基础地理信息数据范围，如图1所示。

图1 不同投影下的基础地理信息数据范围比较Fig.1 Comparing the boundaries of diあerent projections

由于任何投影后的地图，都是由基础地理信息数据按照一定的数学法则从地理坐标转换而来，所以，无论投影后的制图区域如何变化，都可以将其反解到地理坐标上来，从而实现对基础地理信息数据范围的确定，具体确定方法详见算法部分。

1.2 基础地理信息数据的编号

基础地理信息数据的分幅和编号仍以GBT 13989-2012《国家基本比例尺地形图分幅和编号》为标准，采用国际上统一的“经纬线分幅，行列式编号”原则，编号由其所在的行号（字符码）和列号（数字码）组合而成。

1.2.1 1∶1 000 000基础地理信息数据的分幅和编号

国家1∶1 000 000基础地理信息数据分幅和编号的规定：

1）纬向成行：自0°纬线起算，每4°为一行，从赤道至南、北纬88°各有22行，用字母A，B，C，…，V表示。

2）经向成列：从180°经线起算，自西向东每6°为一列，全球分为60列，用阿拉伯数字1，2，3，…，60表示。

3）一个行号和一个列号组成一幅1∶1 000 000基础地理信息数据的编号。

1∶1000 000基础地理信息数据编号的计算公式如下：

式中，[ ]为取整符号，H为行号，L为列号，λ、ψ分别为某点的经度和纬度。

1.2.2 1∶5 000～1∶500 000基础地理信息数据的分幅和编号

1）1∶5 000～1∶500 000基础地理信息数据的分幅均以1∶1 000 000基础地理信息数据为基础，逐次加密划分而成，横为行，由上而下排序；竖为列，由左向右排序。

2）1∶5 000～1∶500 000基础地理信息数据的编号也均以1∶1 000 000基础地理信息数据为基础，由10位代码组成：前三位是1∶1 000 000基础地理信息数据的编号；第四位是比例尺代码，1∶500 000～1∶5 000分别由字母B～H表示；后六位分两组，前一组是图幅行号，后一组是图幅列号，不足三位补零，行号从上而下排序，列号从左到右排序。

1∶5 000～1∶500 000基础地理信息数据编号的计算公式如下：

式中，[ ]为取整符号，H、L分别为1∶1 000 000分幅的基础地理信息数据行号和列号，h、l分别为1∶5000～1∶500 000分幅的基础地理信息数据行号和列号，λ、ψ分别为某点的经度和纬度。△λ、△ψ分别为某比例尺分幅的基础地理信息数据经差和纬差。

2 数据基础与功能设计

本文以中国及周边1∶1 000 000基础地理信息数据库为基础数据，研究基于投影后的制图区域，确定基础地理信息数据的范围，查询基础地理信息数据的编号等算法，下面以兰勃特投影为例，要求算法严密，运算快捷，查询灵活，结果准确。

本文基于投影后的制图区域，从点、线、面三个方面设计了基础地理信息数据编号的可视化查询窗口，编写了点、线、面查询算法模块，并提供了八种基本比例尺选项，当用户选择其中一种比例尺时，可从点、线、面三个方面进行查询。查询结束后释放内存，没有数据冗余。程序界面如图2所示。

图2 编号查询程序设计界面Fig.2 Programmatic interface of querying map sheet numbers

设计流程如图3所示。

图3 设计流程图Fig.3 Design flow-process diagram

2.1 点查询功能

点查询能自动计算某点所在的基础地理信息数据编号，并将该编号所对应的基础地理信息数据范围准确地显示在地图上。本文提供了两种点查询的方法：图上拾取点查询和输入经纬度查询。

2.1.1 图上拾取点

首先，通过双击鼠标左键拾取图面上的一点（或多点），并显示在地图上；然后，点击“查询图幅”按钮，将该编号所对应的基础地理信息数据范围准确地显示在地图上，以此来完成点的查询。

2.1.2 输入经纬度

首先，通过键盘输入某点的经纬度，其形式可以是小数度（D.DD）或者是度分秒（DMS）任意一种；然后，点击“查询图幅”按钮，将该编号所对应的基础地理信息数据范围准确地显示在地图上，以此来完成点的查询。

2.2 线查询功能

线查询能自动计算该线所经过的基础地理信息数据编号，并将编号所对应的基础地理信息数据范围一并显示在地图上。

首先，通过双击鼠标左键获取线（或折线）的节点，并显示在地图上，再双击鼠标右键结束画线。点击“查询图幅”按钮，该线所经过的基础地理信息数据范围及编号将准确地显示在地图上，以此来完成点的查询。

2.3 面查询功能

面查询能自动计算该面所覆盖的基础地理信息数据编号，并将编号所对应的基础地理信息数据范围一并显示在地图上。

首先，通过要素选择工具选择制图主区的边界（如行政区划界线等）；然后，点击主界面的“查询图幅”按钮，弹出“初设比例尺”界面，完成比例尺设计；最后，点击“初设比例尺”界面的“查询图幅”按钮，该面所覆盖的基础地理信息数据范围及编号将准确地显示在地图上，以此来完成面的查询。

3 算法设计与查询结果

本文以兰勃特投影为例，基于投影后的制图区域，从点、线、面三个方面进行探讨，通过算法设计，程序解算，实现了基础地理信息数据范围的确定和编号的可视化查询，以下为解算过程和程序查询结果。

3.1 点查询

在进行点查询时，首先，通过投影转换获取点的经纬度，利用公式（1）、（2）计算该点所在的基础地理信息数据编号。为了将该编号所对应的基础地理信息数据范围准确地显示在地图上，本文提出了归原法，即以图幅的左下角点为图幅原点，推算整个图幅。计算图幅原点的公式如下：

式中，λ1、ψ1分别为该图幅左下角点的经度和纬度，H、L分别为1∶1000000分幅的基础地理信息数据行号和列号，h、l分别为1∶5000～1∶500000分幅的基础地理信息数据行号和列号，△λ、△ψ分别为某比例尺分幅的基础地理信息数据经差和纬差。

由于是在投影后的制图区域上进行分幅，所以，图幅形状不同于投影前的矩形基础地理信息数据，有一定的投影变化，如图4所示。因此，在进行图幅可视化时，不同于投影前只计算矩形基础地理信息数据的四个角点，而是对矩形基础地理信息数据边上的点进行加密（以每边加密100点为例），通过公式（3）获得矩形基础地理信息数据左下角点的经纬度，在此基础上，利用基本比例尺确定的经差和纬差获得每条边上加密点的经纬度，并通过投影变换[4]得到投影后的加密点坐标，连线得到如图4b所示的图幅。

图4 投影前后图幅对比Fig.4 The comparison of map sheet line before and after map projection

在进行点查询时，还要判断以下两种特殊情况：点在两幅图的公共边上时，要显示该点所在共边的两幅图；点在四幅图的公共角点上时，要显示该点所在共点的四幅图。

3.1.1 图上拾取点

图上拾取点首先获取该点的图上坐标，通过投影反解获得该点的经纬度，然后再进行计算，程序运行结果如图5所示。

图5 图上拾取点查询结果Fig.5 The query result of pick up points

3.1.2 输入经纬度

经纬度的输入有两种形式，如图4所示，当选择任意一种形式输入时，程序会自动弹出填写范例，无论哪种形式，都是以小数度的形式参与运算。度分秒（DMS）转化为小数度（D.DD）的公式如下：

式中，dec为小数度形式的经纬度，dms为度分秒形式的经纬度，[ ]为取整符号，△为一个极小偏差量。这两种形式在计算机中转换时，由于计算机运算是二进制，与十进制转换存在极小误差，不会影响点的精度。但是，当十进制与六十进制转换时，误差就会比较大，直接影响到点的精度，所以在进行转换时要人为添加一个极小偏差量。图6为输入的经纬度恰好为四幅图交点的情况。

图6 输入四幅图交点的查询结果Fig.6 The query result of the four corners

3.2 线查询

在进行线查询时，首先，根据线的外包络矩形获得四个极值点坐标，然后，根据选择的比例尺计算图幅的经差和纬差，进行初步查询，其结果如图7所示。

图7 线查询初步分幅结果Fig.7 The tentative query result of lines

由图7可以看出，初步查询的结果存在多余的图幅，为了清除多余的图幅，本文提出了斜率分段-同异侧判别法。首先，对折线进行分段，如图8所示，并计算出每一段的斜率，公式如下：

然后，根据线段的端点坐标、斜率以及图幅的包络线，对每一幅图进行判断，判断该图幅的边是否位于线段的同侧，位于同侧则删除，位于异侧则保留。对图7中所有图幅进行遍历，得到最终结果如图8所示。

图8 线查询结果Fig.10 The query result of lines

3.3 面查询

在进行面查询时，首先，要根据制图区域确定适宜的比例尺和矩形图廓，再根据矩形图廓反解出其所覆盖的基础地理信息数据范围和编号。

当选取制图区域时，点击主界面上的“查询图幅”按钮，程序弹出“初设比例尺”窗口，当用户输入内图廓（宽度与高度）尺寸，程序则根据制图区域的形状和大小自动绘出外接矩形，即内图廓线。同时，程序根据输入的内图廓尺寸和制图区域的实际大小，按照“长边比长边，短边比短边”的比例尺计算法则，计算出该制图区域的横向比例尺与纵向比例尺。用户再根据横向和纵向比例尺，经取整后，给出合适的设计比例尺，程序再根据用户的设计比例尺重新绘制内图廓线，如图9中的矩形区域，经适当偏移和避开微小区域后，该区域即为制图区域所需要的基础地理信息数据范围。

图9 设计比例尺Fig.9 Designing the scale

由于采用的是投影后数据，矩形制图区域投影反解后会变成扇形，所以，根据矩形区域获得的四个极值与该矩形区域投影反解后所求得的四个极值并不相同。

根据矩形区域的四个极值点，通过投影反解获得对应的经纬度，进而求得经差和纬差以及图幅数量。最终，绘制查询图幅范围如图10a黑线所示。从图中可以看出，矩形的四个角点并不对应经纬度的最大最小值，故而查询结果会出现某些地方缺失，而另一些地方却是多余的情况。

图10 投影反算前后结果对比Fig.10 The comparison before and after projection inversion

为了解决这个问题，本文提出了投影反算-图廓内外判别法。首先对矩形边上的点加密并进行投影反算，本文以矩形宽度（或高度）与该比例尺图幅经差（或纬差）的比值为加密的单元数，每个单元加密100个点，加密后的矩形区域如图10b中的绿线所示。加密过程中比较每个加密点的经纬度，得出经纬度的最大最小值，进而求得经差和纬差以及图幅数量，绘制查询图幅范围如图10b黑线所示，从图中可以看出，查询的图幅不存在缺失现象，但却存在着多余现象，所以本文对每个图幅的范围与图廓范围进行比较，判别图幅的位置，位于图廓外的删除，位于图廓内的保留，最后加上图幅编号，最终查询结果如图11所示。

图11 初设比例尺查询结果Fig.11 The query result of initial design scale

4 结束语

本文以兰勃特投影为例，基于投影后的制图区域，从点、线、面三个方面提出了归原法、斜率分段-同异侧判别法和投影反算-图廓内外判别法，这些算法可准确映射各种投影后制图区域所对应的基础地理信息数据范围，实现了基础地理信息数据编号的可视化查询。特别是在编制《中华人民共和国普通地图集》前期的基础地理信息数据库建设与扩建项目中，发挥了精准确定基础地理信息数据范围的作用，提高了基础地理信息数据编号的查询效率和准确度，避免了数据缺失和数据冗余，为基础地理信息数据库驱动《中华人民共和国普通地图集》的编制奠定了数据基础。