(一) 索引矢量文件的建立
首先根据地图覆盖范围绘制一个最小外包矩形(minimum bounding rectangle, MBR)。其次,根据索引图的最大行列数建立空间网格,这里运用了ArcGIS中的“create fishnet”(创建渔网)工具。之后根据空间选择要素的方法,选择与索引图重合的网格,即可生成完整覆盖索引图的网格图层。最后,考虑到图廓外还有很多属性信息,同时每张地图有编号,故在网格图层属性表增设ID列,通过标注每个地图的编号来关联其他属性信息。至此,每幅地图都会对应一个矩形范围,并且可以通过ID列与其他属性(如单独在Excel录入各图的属性数据)实现唯一关联。
(二) 地图抽样配准
对于单幅地图的配准,廖泫铭、张萍、韩昭庆等学者已有详细论述,在此不多做讨论。[20]本文主要考虑批量配准的情况。由于旧图廓没有明确的科学投影方式,只表述为平面坐标,无法在ArcGIS或QGIS等软件中按照投影规则来映射。根据索引图,选取地图总数10%左右且相对分散的地形图作为样本开展基于同名点古今对照的空间配准,利用城墙、道路交叉口、庙宇等明确的特征点与现代电子地图(如天地图)叠加,采用二次多项式等变换方法进行配准。又,索引图中只有地图的邻接情况,分布在地图内部的配准点难以直接映射到索引图中。因此,对于每张配准好的地图,需采集好其角点经纬度信息,建立若干样本地图的角点坐标与索引图中图幅角点的映射关系,即对索引图自身开展空间配准,推演出其他地图的角点坐标信息。根据相邻关系,最多配准50%的量即可得到所有地图的角点信息。得到映射关系后,还需对上一步骤中得到的索引矢量文件进行校正,使栅格图与矢量文件最终匹配。
(三) 利用ArcPy批量配准
在前两步对栅格接图表进行矢量化(未赋予真实坐标信息)并对地图抽样配准之后,笔者根据接图表的相邻关系和已配准地图的角点坐标推测出其他未配准地图的角点坐标信息,即整套图集的批量配准。这一过程使用了Python编程语言及其中ArcPy的站点包。Python是一种通用编程语言,也是一种支持动态输入的解释型语言,适用于交互操作以及一次性程序(即脚本)快速原型制作,同时具有编写大型应用程序的强大功能。ArcPy是一个继承“arcgisscripting”功能模块的Python站点包,可实用高效地通过Python执行地理数据分析、数据转换、数据管理和地图自动化。使用Python和ArcPy,只需要较小的成本便可开发出地理实用程序,用户可以使用由众多Python小群体开发的附加模块。[21]
本文所需的方法涉及ArcPy中多个函数,所编脚本具体包括以下步骤:
1. 确定每张地图四角点图像坐标
地图内图廓外的内容虽在OpenCV的图像处理阶段,经过了统一的透视变换、裁剪,但考虑到绘制、扫描等误差,地图图像仍未必为统一大小(横纵像素值有出入)。因此在配准之前需要计算每张地图四角点所对应的精确图像坐标,这里调用了“GetRasterProperties_management”(获取栅格属性)函数。[22]其中,因地图图像在处理后是一个规则矩形,因此只需计算图像的四至范围,即设置输入栅格的属性分别为TOP、LEFT、RIGHT、BOTTOM即可返回四至范围,从而获得四个角点的图像坐标。
2. 确定每张地图四角点真实地理坐标
通过已校正的矢量索引文件,可计算每个地图内图廓的真实地理坐标,这里使用ArcPy中的“SearchCursor”(查询游标)函数,通过查询对应id的系统自带编码FID,再提取出该FID对应四边形“geometry”类型中的四个角点真实坐标。SearchCursor用于建立从要素类或表中返回记录的访问权限,geometry.points可以返回空间点的位置坐标信息。
3. 批量配准
根据角点的原始图像坐标与真实地理坐标,即源控制点和目标控制点,使用“Warp_management”(扭曲)函数转换栅格数据集,实现循环配准。[23]在转换过程中,默认的变换方法即多项式阶数,将执行仿射变换;默认的重采样方法则为最邻近法,尽可能少更改像素值。事实上,当栅格数据需要根据多项式、矩阵变换来进行系统化校正之时,这一扭曲工具非常实用,它还可通过设置变换方法参数校正极为复杂的变形。在本研究中,笔者选择了默认的仿射变换[24],一则绝大多数地图是根据平面坐标系绘制而成,在绝大多数区域仿射变换可以使变换后的地图与原图的形状、范围相对一致;二则如采用多阶变换,地图或多或少会出现扭曲乃至极大变形,这会对地图阅读产生新的障碍。
(四) 配准方法分析
本文的研究对象是《汉珍数位地图集》,共1945幅地图。因大量地图为旧图廓,未有角点坐标信息,栅格索引表也无经纬度信息,不可避免要开展同名点匹配工作。假设每幅地图配准耗时10分钟,则完成该套图集中旧图廓地图的精确配准共需324.17小时。在本文算法中,处理一张地图的配准时间为28秒(该时间包含从物理磁盘中读取地图文件、写入地图文件的读写过程;处理环境即笔记本电脑配置为i7处理器,32G内存)。事实上,因缺乏角点坐标信息,整套图的批量配准无法直接开展,需按照一定比例抽样精确配准。如抽样比例设为10%,该部分地图配准耗时为32.42小时,之后90%的地图配准耗时为13.62小时。如设地图数量为N,抽样比例为x,总消耗时间(小时)应为T=(N×x×10+N×(1-x)×28/60)/60。至于抽样配准最多需要多少地图的问题,如为连续矩形图幅,最多仅需要1/4比例的抽样配准便可不经数学计算,而是通过边相邻或角相邻等邻接关系直接得到其他未配准地图的角点坐标信息,如图5a所示。因此,按照上面抽样精确配准、批量配准所耗单位时间及计算公式,理论上最多需要92.38小时,相较于整套图集的精确配准,效率也有了数倍提升。此时,角点经纬度均经过考订,精度未见损失,可等同于其他同名点配准方法。
关于新图廓地图能否直接按照角点处坐标进行批量配准的问题,则需对地图精度进一步讨论。栅格索引表虽未见标注,但可根据其中任意一幅地图的角点坐标信息及新图廓空间范围(15'×10')进行全部角点的迅速定位(图5b)。如按照该图集中共86幅新图廓地图来计算,不到1小时即可完成配准,其效率相较人工也有较大提升。但精度方面,CHMap平台发布的《集成》图集配准时的主要空间参考为索引图的经纬线[25];然而,因民国时期测绘技术受限,部分地图甚至由旧图廓直接改绘为新图廓,其角点坐标也有所偏差,根据潘威、万智巍等人的研究知其误差可达2—3千米。[26]如结合一定比例的抽样地图开展同名点配准,再对剩余地图进行批量配准,精度则有显著提升。在地图配准之时,采用二次多项式方法转换,最终地图均方根误差(RMS Error)在0.8像素以内,精度有较大提升。