EN中文

历史地理研究 >

2024 , Vol. 44 >Issue 2: 109 - 122

DOI: https://doi.org/10.20166/j.issn.2096-6822.L20220309

地图学史

基于OpenCV与ArcPy的民国大比例尺地形图批量配准方法——以汉珍数位民国时期五万分一地图集为例

  • 李爽
展开
  • 复旦大学历史地理研究中心,上海 200433

李爽,女,1992年生,陕西渭南人,博士,复旦大学历史地理研究中心青年副研究员,主要从事历史地理信息系统研究。

收稿日期: 2022-09-15

  网络出版日期: 2024-07-22

基金资助

教育部人文社会科学重点研究基地重大项目“数字时代的中国历史地名信息化研究”(22JJD770022)

收起

Batch Geo-registration Method for Large-scale Topographic Map of the Republic of China Based on OpenCV and ArcPy： A Case Study of the 1∶50，000 Scale Atlas Published by Transmission Books and Microinfo

  • LI Shuang
Expand

Received date: 2022-09-15

  Online published: 2024-07-22

Fold

摘要

民国时期,北洋政府及国民政府曾在全国范围组织测绘大比例尺地形图。对数千幅民国地图扫描归档、空间配准并与今天电子地图叠加实现古今对照,虽然繁重但意义重大。基于OpenCV与ArcPy方法库,发展一套针对民国地形图集的自动配准方法,保证精度的同时极大减少了人工,《汉珍数位地图集》配准的过程可资验证。最后,以《中国大陆五万分之一地图集成》为案例进行地图覆盖范围比较分析,以支撑其史料价值发挥与后续研究。

关键词： 地形图; 配准; OpenCV; ArcPy; 数位汉珍

本文引用格式

李爽 . 基于OpenCV与ArcPy的民国大比例尺地形图批量配准方法——以汉珍数位民国时期五万分一地图集为例[J]. 历史地理研究, 2024 , 44(2) : 109 -122 . DOI: 10.20166/j.issn.2096-6822.L20220309

一、 引言

现代遥感研究大多以1972年(美国国家航空航天局的陆地卫星Landsat计划发射元年)为上限,因此古旧地图对于扩充遥感之前时段的长时序研究有重要意义,其自身也可以作为一种重要史料。①(①Michael A. Wulder, David P. Roy, Volker C. Radeloff et al., Fifty years of Landsat science and impacts, Remote Sensing of Environment, 2022, Vol.280, pp.113-195.)就这一问题,满志敏曾指出研究区域人地关系的空间和时间变化,一定要考虑多源数据融合问题,尤其重建100年以上尺度的土地利用时空过程时,应充分利用人文社会资料(如历史文献资料、统计资料)。②(②满志敏: 《小区域研究的信息化: 数据架构及模型》,《中国历史地理论丛》2008年第2辑。)
民国时期,北洋政府和国民政府先后组织了全国范围的大比例尺地形图测绘工作。北洋政府于1914年制定了第一个“十年速测计划”,该计划将测制全国军用地图的任务划归各省执行。在这一阶段,我国测绘的统一基准尚未建立,仍沿用清光绪年间测量的各地经纬点作为平面控制点,继而采用各省设置的假定高程起算点测量,无法构成有效的坐标、高程计算系统,省与省之间也难以拼接。1929年,南京国民政府在原“十年速测计划”的基础上,由参谋本部陆地测量局制定了新的《全国陆地测量十年计划》;1933年,增改为《全国军用、地籍图测量计划纲要》。据统计,1∶5万比例尺地图最终完成了计划的77.2%,有约一半的国土面积未能进行任何比例尺的地形测绘。①(①《中国测绘史》编辑委员会编: 《中国测绘史》第2卷《明代—民国》,测绘出版社1995年版,第533—534页。)这一阶段,机构采用大三角测量法,辅以若干天文测量、基线测量和水准测量等,并对地图投影、图廓范围进行了科学规定,如采用兰伯特正形割圆锥投影,分幅方式从旧图廓分幅(36cm×46cm矩形图廓)转为经纬度分幅(15'×10'),地图内容则采用符号规范表示。这批地图对当时地形、河流、城市进行了精确描绘,可直接用于城市、水文、地形等专题研究,具有相当重要的学术和人文价值。民国时期,由于测量主体变更、分省组织、时局动荡等多重原因,地图馆藏分散,加之时间久远,地图纸张也难免老化、变黄甚至霉变、破碎,因此对其开展抢救保护、数字化存档十分必要。
国内外对于民国时期地形图集的整理与出版,亦不乏代表性成果。2018年,台湾汉珍数位图书股份有限公司出版了五套民国时期五万分一地图集,包括《民国时期五万分一: 河北地图集》(3册)、《民国时期五万分一: 河南地图集》(4册)、《民国时期五万分一: 山东地图集》(4册)、《民国时期五万分一: 山西地图集》(4册)、《民国时期五万分一: 陕甘地图集》(2册),共17册(以下统称“《汉珍数位地图集》”)。②(②侯巧芬主编: 《民国时期五万分一河北地图集》,汉珍数位图书股份有限公司2018年版;侯巧芬主编: 《民国时期五万分一河南地图集》,汉珍数位图书股份有限公司2018年版;侯巧芬主编: 《民国时期五万分一山东地图集》,汉珍数位图书股份有限公司2018年版;侯巧芬主编: 《民国时期五万分一山西地图集》,汉珍数位图书股份有限公司2018年版;侯巧芬主编: 《民国时期五万分一陕甘地图集》,汉珍数位图书股份有限公司2018年版。)地形图反映了不同时期的测绘标准、技术情况,例如该套地图集中的《民国时期五万分一: 陕甘地图集》,既包括了旧图廓地图,也包括了带有经纬度信息的新图廓地图。两种不同分幅方式的地图可能同时出现在同一册地图集,在同一张索引图上不同图廓也可被清晰区分(图1)。通过对地图的图幅情况整理汇总可知,在全部1945幅地图中,有86幅地图在内图廓角点处标有经纬度,即新图廓,约占总地图数的4.42%,具体统计情况见表1。此前,日本科学书院曾于1986—1998年间先后出版发行了八册《中国大陆五万分之一地图集成》(以下简称“《集成》”),包括分布在20个省份的4088幅民国时期地形图。③(③[日] 科学书院: 《中国大陆五万分之一地图集成》,科学书院1986—1998年版。)任玉雪与邓发晖在《<中国大陆五万分之一地图集成>所收地图来源分析》一文中,对《集成》中地形图的测图、制版、发行、注记、等高线等信息和地形图图表开展了系统分析,认为《集成》地图主要为中国政府测绘,被日本窃取之后加以改制和发行。④(④任玉雪、邓发晖: 《<中国大陆五万分之一地图集成>所收地图来源分析》,《历史地理研究》2020年第3期。)车群等对《集成》的数字化过程、CHMap作为地形图的载体以及作为互联的连结进行了介绍,并从空间知识论的视角对未来进行了思考。⑤(⑤车群、林农尧、陈诗沛: 《以空间连接时间: 以陆地测量军事地形图为基础的地图(图像)共享与互操作》,《数字人文研究》2021年第4期。)此外,杭州市档案馆于2013年整理出版了《民国浙江地形图》,收录了民国时期绘制的浙江五万分之一比例尺地形图345幅。①(①杭州市档案馆编: 《民国浙江地形图》,浙江古籍出版社2013年版。)
图1 《民国时期五万分一陕甘地图集》的封面、检索地图与全幅简图②(②《民国时期五万分一陕甘地图集》检索地图的左上部分为甘肃,右侧为陕西中部区域,其中甘肃部分为新图廓,陕西为旧图廓,新旧图廓之间存在相交情况。)

资料来源: 侯巧芬主编《民国时期五万分一陕甘地图集》(汉珍数位图书股份有限公司2018年版)。

表1 《汉珍数位地图集》地图编号及经纬度信息统计
图 集 名 称 地 图 编 号 含有经纬度信息的地图编号
《民国时期五万分一河北地图集》第1册 1—118 —
《民国时期五万分一河北地图集》第2册 119—237 —
《民国时期五万分一河北地图集》第3册 238—356 —
《民国时期五万分一河南地图集》第1册 1—121 —
《民国时期五万分一河南地图集》第2册 122—242 —
《民国时期五万分一河南地图集》第3册 243—363 —
《民国时期五万分一河南地图集》第4册 364—484 —
《民国时期五万分一山东地图集》第1册 1—109 —
《民国时期五万分一山东地图集》第2册 110—218 —
《民国时期五万分一山东地图集》第3册 219—327 —
《民国时期五万分一山东地图集》第4册 328—436 429—436
《民国时期五万分一山西地图集》第1册 1—112 —
《民国时期五万分一山西地图集》第2册 113—224 —
《民国时期五万分一山西地图集》第3册 225—336 —
《民国时期五万分一山西地图集》第4册 337—448 —
《民国时期五万分一陕西地图集》第1册 1—111 —
《民国时期五万分一陕西地图集》第2册 112—221 144—221
空间信息技术为传统纸质地图的存档与研究带来新的活力,尤其为历史时期的纸图提供了存档、检索、共享、古今对照等功能,进而为挖掘历史地理信息、提升资料利用率提供了技术支撑。国内外多个大学、图书馆等存有一定数量的民国时期地形图,可线上获取,在此仅举例一二。陕西省测绘档案资料馆馆藏民国时期地形图采用地形图叙录形式编制,为民国时期地形图档案的整理奠定了基础。该馆还建成、开放了陕西省测绘档案资料馆内藏陕西省民国时期地形图叙录系统,其中,五万分一比例尺地形图有两批,分别为关中、陕南地区五万分一地形图161幅(制图时间为1915—1948年),陕甘宁边区五万分一地形图共82幅(制图时间为1948—1949年)。①(①施小溪、王立言: 《民国时期地形图叙录编纂与展示方法研究》,《兰台世界》2019年第2期;董延红、焦惊眉、高瑛等: 《民国时期地形图保护与抢救的必要性与可行性分析》,《测绘技术装备》2013年第1期。)在数字存档的基础上,台北“中研院”将所藏地图进行拼接、配准,可以加载至GIS软件中调用。②(②廖泫铭、范毅军: 《中华文明时空基础架构: 历史学与信息化结合的设计理念及技术应用》,《科研信息化技术与应用》2012年第4期;中华文明之时空基础架构,http://ccts.ascc.net。)上海交通大学历史系联合德国马普斯·普朗克科学史研究所,以CHMap为主要数据开展配准及可视化工作,同时提供WMTS服务,其中,配准时的主要空间参考为索引图中的经纬线。③(③CHMap, https://chmap.mpiwg-berlin.mpg.de/lgtu-new/;车群、林农尧、陈诗沛: 《以空间连接时间: 以陆地测量军事地形图为基础的地图(图像)共享与互操作》,《数字人文研究》2021年第4期。)丝绸之路历史地理信息开放平台上发布有苏南浙北、丝绸之路沿线的民国地形图,可直接加载于现代电子底图之上。④(④丝绸之路历史地理信息开放平台,http://www.srhgis.com/。)复旦大学历史地理研究中心发布的中国历史地理信息平台,在古旧地图模块上线了《汉珍数位地图集》以及《集成》配准的成果,在时空框架模块可直接加载,并通过设置透明度开展古今对照。⑤(⑤中国历史地理信息平台,https://timespace-china.fudan.edu.cn/。)在地图配准功能方面,国内外一些古旧地图在线平台也提供了相关支撑,如David Rumsey地图中心提供Georeferencer模块实现在线配准,丝绸之路历史地理信息开放平台、中国历史地理信息平台也均支持地图的在线配准。上述平台均已实现对单张地图的配准功能,并可与现代电子地图进行叠加,开展古今对照。而当配准对象是一套卷帙浩繁、规范完整的地图集时,部分配准流程就具备了进行自动化处理的基础。
民国时期地形图的精度情况是其史料价值发挥的前置问题。地图学家曾世英在《陆地测量局制图工作的检讨》一文中指出,民国地形图旧图廓所涉及的平面坐标系在制图区域狭小时尚可采用,区域稍大时平面坐标就会发生错误。①(①曾世英: 《陆地测量局制图工作的检讨》,《地质论评》1938年第1期。)民国时期旧图廓地图采用的是各省自设的平面坐标系,角点处无经纬坐标信息,且直接由各省自行制定起点、标高,缺乏科学的投影说明,不同省之间的地图难以拼合,难以在GIS软件中进行地图投影换算。即便是采用了经纬度分幅的新图廓地图,角点处标明的经纬度与实际经纬度也有一定偏差,根据潘威、万智巍等人的研究知其误差范围可达2—3千米。②(②潘威、满志敏: 《大河三角洲历史河网密度格网化重建方法——以上海市青浦区1918—1978年为研究范围》,《中国历史地理论丛》2010年第2辑;万智巍、贾玉连、洪祎君等: 《基于历史地图与遥感影像的近百年来长江荆江段河道演变》,《地理科学》2019年第4期。)因此,对不含经纬坐标信息的旧图廓配准,或对新图廓地图精确配准,必须要结合现代电子地图、遥感影像,确认特定地物的古今位置信息,如城墙、道路交叉口、文庙等古建筑或重要的地理要素,再进行空间配准。
这般精细的比对、考证适用于地图数量较少的场景;当地图数量达到数千张时,手动考订各张地图同名点将成为一项十分繁重的工作。经过笔者的研究,发现《汉珍数位地图集》具有地形图样式相对固定、绘制清晰的特征,特别是分省索引图的接图信息明确(边相邻的地图可以共用一条边的两个角点坐标,角相邻的地图可以共用一个角点坐标),具有批处理的资料基础。该图集如果角点经纬度信息精确度允许,就可以不进行古今同名点对照,直接按照转换关系批量配准。令人遗憾的是,该套地图集的索引图缺乏《集成》中索引图的经纬度线,即接图表仅是对原地图之间相对位置的说明,未标明明确的地理坐标,需要利用接图表进行配准点换算。
综上,民国地形图存在旧图廓缺乏地图投影、新图廓测绘误差等问题,导致难以构建准确的测绘图幅框架,进而限制了大量地形图的高效配准及应用。本文以科学性与实用性兼备为原则,先选取《汉珍数位地图集》为样例图集,以该图集的配准、空间范围考证为主要研究目的,引入Python在ArcGIS中的支持、OpenCV的图像学等多种方法,批量裁剪地图的边框以避免叠盖。再根据抽样思想,提取一定比例地形图进行同名点配准,获取矩阵变换公式求取参数。在此基础上,充分利用索引表中的接图信息即相对位置关系,推演得到所有地形图的角点坐标,极大提升了工作效率。

二、 基于OpenCV的地图拼接与图廓检测

在地图集印刷过程中,由于版面限制,一幅地图被分为了两张,因此在地图配准前首先要将被拆分为两部分的地图拼接成一张完整的地图。同时,为了解决图幅间相互叠盖的问题,需要识别出地图的图廓并裁剪。由于上述步骤均不涉及地理坐标,可视为纯粹的计算机图像学问题。鉴此,笔者尝试引入OpenCV(Open Source Computer Vision Library)库以解决上述问题,技术路线参见图2。OpenCV是一个开源计算机视觉库,由英特尔公司发起研发,用于开发实时的图像处理、计算机视觉、模式识别等程序。③(③OpenCV, https://opencv.org/。)OpenCV轻量级且高效,底层由C++语言撰写,同时提供了Python、Java、Matlab等接口,在工业领域与学术研究中都有着广泛的应用。
图2 基于OpenCV的地图拼接与图廓检测技术路线

(一) 基于关键点匹配的地图拼接

为了阅读的连贯性,被拆分的地图通常在内容上有所重叠(图3a)。又因为在印刷出版、扫描过程中均可能存在扭曲、交叠、倾斜等情况,每张地图的拼接范围或有区别,不能根据重叠区域的像素范围直接拼接,需要匹配两张地图中同名特征点进行像素级的拼接。图像拼接的方法在计算机视觉、图像处理方面是重要的研究方向。本文所用方法在图像层面(不涉及地理坐标)调用相关图像拼接函数实现批处理,将一幅或多幅存在重复区域、相同特征点的图像配准、投影、融合,再进行组合,即建立待拼接图像与原图之间的对应关系。最后找出相应的数学转换模型,将其转换至同一参考系,形成一幅新的、更大的图像。
图3 膨胀操作与轮廓识别原理示意
图像配准是其中最核心的内容,目前主要有基于灰度和基于特征的配准算法,后者在准确率、效率方面表现更好。具体来看,图像同名点配准过程通过尺度不变特征转换算法(scale-invariant feature transform,SIFT)等来构建尺度空间,提取局部不变描述符,即特征点;定位特征点,确定主方向以确保转换前的鲁棒性(robustness);计算图像之间的单应性矩阵,进行投影变换;最终返回布尔值(boolean),拼接成功时为true,否则为false,输出生成的拼接图。由于不能保证两张图色彩一致,在坐标变换过程中,应当对两幅图的同名点像素取平均值,以消除与两张图中其他区域亮度的不稳定度。本文采用的方法相较于重叠区域的直接拼接,具有更好的适用性和稳定度。这里具体可调用OpenCV的“stitch()”函数,实现批处理。

(二) 图廓的识别与裁剪

一张地图可分为外图廓外、内外图廓之间、内图廓内三个部分。其中,外图廓外含有图名、制版时间、制图机构、图示及标高等信息,内外图廓之间可能存在经纬度(位于新图廓角点处)、等高线高程数值、道路延伸方向等信息,内图廓内则是地图的主体内容。以《民国时期五万分一河北地图集》中的北平一图为例(地图编号为88),图廓外包括以下信息: 图名“北平”,图示“据民国三年所定之地形原图图式”,标高“自秦皇岛之中等潮位以米远起算”,制图单位“参谋本部陆地测量总局”,时间“中华民国十八年十一月缩制,二十一年十二月制版”。
在空间处理时,只需对内图廓内的地图内容开展配准(外图廓外为属性信息,需要单独采集),因此第一步就要裁掉内图廓之外的内容避免覆盖。图幅边框的准确识别是裁剪时最为重要的工作,这一工作涉及图像形态学多个算法,具体如下。

1. 基于膨胀操作的边缘增强

由于地图在印刷、扫描过程中难免存在畸变,特别是图廓的线条可能有间断、不闭合的现象,所以要结合形态学滤波方法,先对图像进行二值化,在此基础上进行膨胀操作,使边缘闭合。膨胀操作是指通过核算子,寻找出原图像中的极大值区域,从而填充一些遗漏的像素值,表现在图像上则是“扩大”,填充原图的孔洞、缝隙,使图像更为连贯。①(①王树文、闫成新、张天序等: 《数学形态学在图像处理中的应用》,《计算机工程与应用》2004年第32期。)与之相对应的操作即对偶运算则是腐蚀。经过膨胀操作,图廓的线条会更为闭合,从而有助于后续对图廓的整体识别。

2. 外图廓识别

在图像层面上,地图图廓是由一系列点环绕形成的闭合多边形(由于绘制、制版、扫描等影响,大多不是规则矩形)。在图像二值化的基础上,从有数值的点出发,沿着顺时针或逆时针方向,便可识别出由一组连续点构成的集合,即图中的轮廓,这一过程被称为边界跟踪,原理见图3b。②(②Satoshi Suzuki, Keiichi Abe., Topological structural analysis of digitized binary images by border following, Computer Vision, Graphics, and Image Processing, 1985, Vol. 3, p.396.)在OpenCV库中,有用于寻找二值图像中轮廓的函数,即“findContours()”函数。笔者进一步对轮廓围绕的区域面积开展数值计算,将数值由大到小排序,其中面积最大区域对应的轮廓便是地图的外图廓。

3. 透视变换

考虑到栅格化地图经过了绘制、制版、扫描等步骤,外图廓极有可能是不规则矩形,因此需要根据四角点定位透视变换,即转为上下两条边、左右两条边像素范围一致的规则矩形。透视变换是一种常用的图像变换方法,是对图像投影与原投影之间的转换,即将一个图像根据特定的映射关系,投至一个新的规则视平面。事实上,外图廓虽然不是规则的矩形,但倾斜、变形程度相较于图幅的范围并不大,这里对可能引起的细小偏离暂不讨论。

4. 外图廓裁剪

重新透视变换后的地图,外图廓成了一个长、宽像素一致的规则矩形。因此,可直接选取四个外图廓角点的新像素值,根据一个像素范围(矩形形状),裁剪掉与地图空间信息无关的外图廓外区域。

5. 内图廓的识别与裁剪

在裁剪掉外图廓外的内容后,需要对内图廓再次识别与裁剪。识别内图廓与外图廓的原理及方法并无二致,所以在后续识别、裁剪内图廓时,只需要重复步骤2、3、4,再识别出内图廓的边框、透视变换以及裁剪。
至此,内图廓以外的地图内容均被裁剪完毕。上述方法将扫描后的栅格地图视为一种图像,从图像学视角尝试开展地图拼接、边框识别、透视变换、裁剪等数据处理工作,并通过Python调用OpenCV,实现了4小时内完成1945幅地图的批处理(笔记本电脑配置为i7处理器,32G内存)。这一流程如使用手动处理,按照上述流程在Photoshop等专业图像处理软件中操作,采用人工操作(按熟练工每幅图10分钟计算),至少超过300小时。因此可以认为该方法在保证精度的同时极大降低了人工的工作量。

三、 基于ArcPy的快速配准

裁剪好内图廓外的区域后便要对这些地图进行配准。各省地图数量虽不一致,但都有数百幅,河南、山东、山西等省更分别有400余幅。对于这类缺乏明确地图投影和经纬度信息的地图,手动配准难度较高。因此,笔者利用样本抽样思想以及接图表的相对位置信息,通过引入ArcPy站点包,设计了一套快速配准的方法。该方法的具体思路是以同名点匹配的考订方法配准一定比例的地形图,再以索引图作为桥梁,通过索引图中反映的接图关系推演出其他地图的角点坐标信息,具体技术路线如图4所示。
图4 基于ArcPy的快速配准方法

(一) 索引矢量文件的建立

首先根据地图覆盖范围绘制一个最小外包矩形(minimum bounding rectangle, MBR)。其次,根据索引图的最大行列数建立空间网格,这里运用了ArcGIS中的“create fishnet”(创建渔网)工具。之后根据空间选择要素的方法,选择与索引图重合的网格,即可生成完整覆盖索引图的网格图层。最后,考虑到图廓外还有很多属性信息,同时每张地图有编号,故在网格图层属性表增设ID列,通过标注每个地图的编号来关联其他属性信息。至此,每幅地图都会对应一个矩形范围,并且可以通过ID列与其他属性(如单独在Excel录入各图的属性数据)实现唯一关联。

(二) 地图抽样配准

对于单幅地图的配准,廖泫铭、张萍、韩昭庆等学者已有详细论述,在此不多做讨论。①(①廖泫铭、范毅军: 《中华文明时空基础架构: 历史学与信息化结合的设计理念及技术应用》,《科研信息化技术与应用》2012年第4期;张萍、双静如: 《民国<西京市区>图的数字化、数据提取和应用价值》,《图书馆论坛》2021年第11期;韩昭庆、杨霄、刘敏等: 《康熙<皇舆全览图>长城以南地区绘制精度的空间分异》,《清华大学学报(哲学社会科学版)》2021年第3期。)本文主要考虑批量配准的情况。由于旧图廓没有明确的科学投影方式,只表述为平面坐标,无法在ArcGIS或QGIS等软件中按照投影规则来映射。根据索引图,选取地图总数10%左右且相对分散的地形图作为样本开展基于同名点古今对照的空间配准,利用城墙、道路交叉口、庙宇等明确的特征点与现代电子地图(如天地图)叠加,采用二次多项式等变换方法进行配准。又,索引图中只有地图的邻接情况,分布在地图内部的配准点难以直接映射到索引图中。因此,对于每张配准好的地图,需采集好其角点经纬度信息,建立若干样本地图的角点坐标与索引图中图幅角点的映射关系,即对索引图自身开展空间配准,推演出其他地图的角点坐标信息。根据相邻关系,最多配准50%的量即可得到所有地图的角点信息。得到映射关系后,还需对上一步骤中得到的索引矢量文件进行校正,使栅格图与矢量文件最终匹配。

(三) 利用ArcPy批量配准

在前两步对栅格接图表进行矢量化(未赋予真实坐标信息)并对地图抽样配准之后,笔者根据接图表的相邻关系和已配准地图的角点坐标推测出其他未配准地图的角点坐标信息,即整套图集的批量配准。这一过程使用了Python编程语言及其中ArcPy的站点包。Python是一种通用编程语言,也是一种支持动态输入的解释型语言,适用于交互操作以及一次性程序(即脚本)快速原型制作,同时具有编写大型应用程序的强大功能。ArcPy是一个继承“arcgisscripting”功能模块的Python站点包,可实用高效地通过Python执行地理数据分析、数据转换、数据管理和地图自动化。使用Python和ArcPy,只需要较小的成本便可开发出地理实用程序,用户可以使用由众多Python小群体开发的附加模块。②(②什么是ArcPy?https://desktop.arcgis.com/zh-cn/arcmap/latest/analyze/arcpy/what-is-arcpy-.htm;杨柳: 《一种基于ArcPy的地图数据快速分幅导出方法》,《测绘与空间地理信息》2022年第1期。)
本文所需的方法涉及ArcPy中多个函数,所编脚本具体包括以下步骤:

1. 确定每张地图四角点图像坐标

地图内图廓外的内容虽在OpenCV的图像处理阶段,经过了统一的透视变换、裁剪,但考虑到绘制、扫描等误差,地图图像仍未必为统一大小(横纵像素值有出入)。因此在配准之前需要计算每张地图四角点所对应的精确图像坐标,这里调用了“GetRasterProperties_management”(获取栅格属性)函数。③(③获取栅格属性,https://desktop.arcgis.com/zh-cn/arcmap/latest/tools/data-management-toolbox/get-raster-properties.htm;GetRasterProperties_management函数是ArcPy提供的从元数据和栅格数据集的相关描述性统计数据中检索信息的函数,具体调用方法为arcpy.GetRasterProperties_management(in_raster,{property_type},{band_index}),in_raster为待检索的,property_type为要从输入栅格获取的属性,band_index指从哪个波段获取属性,默认为第一个波段。)其中,因地图图像在处理后是一个规则矩形,因此只需计算图像的四至范围,即设置输入栅格的属性分别为TOP、LEFT、RIGHT、BOTTOM即可返回四至范围,从而获得四个角点的图像坐标。

2. 确定每张地图四角点真实地理坐标

通过已校正的矢量索引文件,可计算每个地图内图廓的真实地理坐标,这里使用ArcPy中的“SearchCursor”(查询游标)函数,通过查询对应id的系统自带编码FID,再提取出该FID对应四边形“geometry”类型中的四个角点真实坐标。SearchCursor用于建立从要素类或表中返回记录的访问权限,geometry.points可以返回空间点的位置坐标信息。

3. 批量配准

根据角点的原始图像坐标与真实地理坐标,即源控制点和目标控制点,使用“Warp_management”(扭曲)函数转换栅格数据集,实现循环配准。①(①扭曲,https://desktop.arcgis.com/zh-cn/arcmap/latest/tools/data-management-toolbox/warp.htm;Warp_management函数的功能是使用源控制点和目标控制点转换栅格数据集,调用方法为arcpy.management.Warp(in_raster, source_control_points, target_control_points, out_raster, {transformation_type}, {resampling_type}),其中参数分别为要转换的栅格、要扭曲的栅格坐标、待扭曲的源栅格坐标系、要创建的栅格、变换方法、重采样方法。)在转换过程中,默认的变换方法即多项式阶数,将执行仿射变换;默认的重采样方法则为最邻近法,尽可能少更改像素值。事实上,当栅格数据需要根据多项式、矩阵变换来进行系统化校正之时,这一扭曲工具非常实用,它还可通过设置变换方法参数校正极为复杂的变形。在本研究中,笔者选择了默认的仿射变换②(②仿射变换,在图像映射中,是指一个图像经过一次线性变换、一个平移,转换为另一个图像,经过该变换前后的图像会保持共线性和比例特性。),一则绝大多数地图是根据平面坐标系绘制而成,在绝大多数区域仿射变换可以使变换后的地图与原图的形状、范围相对一致;二则如采用多阶变换,地图或多或少会出现扭曲乃至极大变形,这会对地图阅读产生新的障碍。

(四) 配准方法分析

本文的研究对象是《汉珍数位地图集》,共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小时,相较于整套图集的精确配准,效率也有了数倍提升。此时,角点经纬度均经过考订,精度未见损失,可等同于其他同名点配准方法。
图5 配准过程分析
关于新图廓地图能否直接按照角点处坐标进行批量配准的问题,则需对地图精度进一步讨论。栅格索引表虽未见标注,但可根据其中任意一幅地图的角点坐标信息及新图廓空间范围(15'×10')进行全部角点的迅速定位(图5b)。如按照该图集中共86幅新图廓地图来计算,不到1小时即可完成配准,其效率相较人工也有较大提升。但精度方面,CHMap平台发布的《集成》图集配准时的主要空间参考为索引图的经纬线①(①车群、林农尧、陈诗沛: 《以空间连接时间: 以陆地测量军事地形图为基础的地图(图像)共享与互操作》,《数字人文研究》2021年第4期。);然而,因民国时期测绘技术受限,部分地图甚至由旧图廓直接改绘为新图廓,其角点坐标也有所偏差,根据潘威、万智巍等人的研究知其误差可达2—3千米。②(②潘威、满志敏: 《大河三角洲历史河网密度格网化重建方法——以上海市青浦区1918—1978年为研究范围》,《中国历史地理论丛》2010年第2辑;万智巍、贾玉连、洪祎君等: 《基于历史地图与遥感影像的近百年来长江荆江段河道演变》,《地理科学》2019年第4期。)如结合一定比例的抽样地图开展同名点配准,再对剩余地图进行批量配准,精度则有显著提升。在地图配准之时,采用二次多项式方法转换,最终地图均方根误差(RMS Error)在0.8像素以内,精度有较大提升。
总体来看,本方法在保证精度的同时,可以极大提升配准效率。其中,抽样比例与效率成反比关系,与精度成正比关系。但如需进一步提升效率,可适当降低抽样比例,此时应根据实际需求予以平衡。

四、 方法应用与讨论——与日本科学书院地图覆盖范围的对比

在批量配准后,笔者对《汉珍数位地图集》共1945幅地图的覆盖范围进行可视化,得知共覆盖河北、河南、山东、山西、陕西、甘肃6个省区(图6a)。《中国大陆五万分之一地图集成》(八卷本)中共包含4 088幅五万分一比例尺的地形图,笔者对《集成》的覆盖范围可视化后,发现其覆盖20个省区(图6b)。从覆盖范围看,《集成》覆盖的省份远多于《汉珍数位地图集》,但对一些特定区域,如山东、山西绝大多数区域以及河北南部、河南西南局部以及甘肃部分区域,后者的地图可以进行有效补充。相较CHMap平台的配准方法,本文重点讨论的方法采纳了更多图中重要地理要素的坐标,并未全然依赖索引图的经纬线。
图6 地形图覆盖范围与今北京故宫区域
上述讨论仅围绕两套地图的空间范围,若要全面、完整论述相同图名或相同位置地图之间的具体区别,则需要对地图空间信息、图廓外的文字信息逐一考订、比对。例如,民国时期,日本对于窃取的大陆地图进行更新、再版,地图的制图单位、版式、绘制方法等皆因此有所出入。笔者通过检索两套图集中的现北京故宫区域,发现二者在绘制细节上出入较大,地图的图幅范围、出版单位、时间亦难统一。其中,《汉珍数位地图集》的北京地图(图6c)缩制于民国十八年十一月(1929年11月),民国二十一年十二月(1932年12月)制版;而《集成》的北平地图(图6d)测于明治四十三年七月(1910年7月),次年制版、发行。可以看出,对于近代大比例尺地形图的测绘、脉络、地图比对等话题的研究,尚有广阔的空间。

五、 结语

民国地形图因存在旧图廓缺乏科学的地图投影、新图廓测绘误差等问题,难以在GIS软件中直接投影换算,造成无法高效配准大量地形图的现实瓶颈,进而制约其研究价值、史料价值的发挥。如需要空间配准,应结合现代地图、遥感影像,确认特定地物的古今位置信息,如城墙、道路交叉口、文庙等古建筑或重要的地理要素,进行空间校正。但当地图数量达到数千幅时,工作量十分繁重。基于上述背景,本文以抽样及推演作为指导思想,通过配准特定比例的地图,再充分利用索引图的接图信息,推演得到其他地图的图幅范围,最终形成一套基于OpenCV与ArcPy方法库针对民国地形图集的批量配准方法。该方法相较逐一手动配准极大降低了工作量,相较直接基于索引图的配准(针对新图廓)也增加了很多考证的坐标位置信息,因此精度更为可靠。
笔者选择汉珍数位图书出版的民国时期五万分一比例尺地形图进行了方法验证,并与日本科学书院出版的《中国大陆五万分之一地图集成》覆盖范围以及特定地图开展比较分析。分析认为,《汉珍数位地图集》在空间范围上可以有效补充《集成》的覆盖范围;而在地图的图幅、属性、绘制细节上,两者之间仍存在考订、分析的空间。在省级或更大尺度民国地形图的配准过程中,该配准方法能够充分利用接图表信息、利用特定比例地图进行批量配准,原则上具备可推广性。但考虑到民国时期地图繁杂、质量不一,存在文字与图框叠盖、地图破损图廓不全的现象,不少地图因保存不善,字迹不清,甚至存在纸张老化、变黄、霉变、破损等问题,在配准之前需要检查地图质量,对存在问题的地图进行图像层面的预处理工作。因此,增强该算法在特殊情况下的适用性、推动其工作效率的进一步优化,将成为未来的重要研究方向。
Options
文章导航

/

〈 〉