一、概述
根据《四川省志·测绘志》,1927至1934年,四川省陆地测量局实测四川省1:10万地形图263幅。依据1913年制定的地形图图式,图幅为36厘米×46厘米,无经纬度和方里网。这批地图经历了多次制版印刷:四川省陆地测量局于1934至1940年所印图,四川省图书馆特藏部存189幅;日本参谋本部陆地测量部于1940至1942年依据1935年修正的地形图图式制版,1942年复制,共263幅;中国人民解放军东北军区司令部1949年翻印,四川省和陕西省测绘档案资料馆保存全套印刷图[1]。这些地图用密集而细腻的等高线刻画,县界、省界清晰可辨,详细记录了抗战前夕四川全省的行政建置、山川城镇、交通道路和地貌形态,是了解民国时期巴蜀大地直观的地理文献,具有较高的历史价值。然而,对于数百张分幅地图,研究者需要反复翻找、不断切换不同编号才能拼凑出完整地理信息。如何将263张地图图像,逐一矫正、裁切到统一尺寸后,再按正确的行列顺序排列成一幅完整的、更为直观的四川全境图?这正是本文要探讨的问题。通过一套完整的Python批量处理流程,实现从零散的分幅地图扫描件自动拼接为可供流畅整体浏览的全图[2]。
二、批量处理流程
关于类似地图拼接问题,前人已有丰富的研究成果。复旦大学历史地理研究中心李爽在汉珍数位民国时期五万分之一地图批量配准过程中,提出了一套由外框到内框,利用OpenCV库中findContours()函数在经过膨胀操作的二值化图像中查找闭合多边形,并以其中区域面积最大者再经过透视变换以完成图框识别和裁切的成熟方案[3]。但本文研究的这批地图中,除了常规的因绘制或印刷产生的边框空洞和缝隙,尚存在不少因省界边缘区域地图少量内容超出图框,或因纸张撕裂产生较大面积打断内外图框的情形(图1)。通过膨胀操作填充较大尺度的断裂区域以形成闭合多边形存在较大困难。在借鉴其从外框到内框逐层渐近识别检测并进行透视变换的思路后,本文优化了对不连续边框的识别方法:通过传统的霍夫变换(Hough Line Transform,见小贴士1)检测得出候选直线并进行分类以识别图框线,再两两相交求出四个角点后执行透视变换。
小贴士1 霍夫变换 (Hough Transform)
通俗理解: 图像中的“找线神器”。
作用: 哪怕老地图因为年代久远出现了断线、撕裂,或者扫描时角度歪了,这个算法也能从一堆杂乱的像素点中,精准地把地图的黑边框“找”出来,并计算出它的精确位置。
图1 地图边框不连续的情形
基于此,本文利用Python自动批量拼接上述地图的完整处理流程共分为五个步骤:步骤一,粗裁去边;步骤二,用霍夫变换识别外框线并矫正透视投影;步骤三,裁剪外框;步骤四,再次用霍夫变换识别内框线并矫正透视投影;最后步骤五,将全部263幅图按原始编号排布为连续大图,并提供两种输出方式:一是基于DeepZoom技术的Web端交互式地图,二是基于QGIS供下一步进行地理配准(见小贴士2)和叠加分析的VRT虚拟栅格(见小贴士3)(图2)。
小贴士2 地理配准 (Georeferencing)
通俗理解: 给老地图“定坐标”。
作用: 目前的拼接只是把图片拼在一起,但电脑不知道它对应地球上的哪个位置。地理配准就是寻找老地图上的特征点(比如老县城的位置)和现代卫星地图上的对应点,把老地图“拉伸”并贴合到地球上,这样才能进行古今对比分析。
小贴士3 VRT虚拟栅格 (Virtual Raster)通俗理解: 一个“虚拟的大拼图文件”。
作用: 在专业的GIS(地理信息系统)软件中,不需要真的把263张图片合并成一个巨大的文件(那样会占满硬盘),VRT只是一个“目录索引”,告诉软件这263张图应该摆在什么位置,软件读取时就像在读一张完整的地图。
图2 自动批量拼接流程图
步骤一:粗裁去边
部分地图原始扫描件的四个边缘留有空白区域,这些区域来自扫描仪的边缘留白,并非地图本身的实际内容。如果不先将其去除,后续的图框检测和透视校正将受到干扰——空白的纸面边缘本身也是高对比度的线条,在外框检测中容易与地图边框混淆。因此,需要首先裁剪掉这些区域。
利用Pillow(PIL)库的Image.crop()对所有地图执行简单的等距裁切,沿四个方向各切除120像素。依据对全批次扫描件边距的统计,绝大多数图幅的空白区域宽度在120像素内,因此基本可去掉大部分无效空白而不伤及地图主体(图3)。
图3 步骤一示意图
步骤二:识别外框线并矫正透视投影
在粗裁去边之后,地图主体部分被保留了,但其在扫描时几乎不可避免地存在角度倾斜与轻微透视畸变。本步骤目标为识别地图最外框粗线,并以此范围按固定比例的分辨率进行透视矫正。核心思路为在缩至一半的小图上用霍夫变换检测直线,然后把角点坐标映射回原图。具体流程如下(图4):
1. 原始图像等比缩小:由于直接在高分辨率的原图上运行霍夫变换会消耗大量内存并拖慢速度,因此先将图像缩小一半(使用INTER_AREA插值保持降采样质量),在较小空间内完成所有计算后,再将结果映射回原图坐标系。这样做的好处有三:一是显著降低内存占用;二是为每步骤的多进程并行处理加快效率;三是能天然过滤部分印刷噪点和纸张皱褶引起的虚假边缘线索。
2. 灰度化、高斯模糊和Canny边缘检测:对模糊化处理的灰度图进行Canny边缘提取。
3. HoughLines直线检测与分类:查找所有可视为边框的直线段,每一线的极坐标参数(ρ,θ)会被记录下来。对所有检测到的直线按倾斜角度分类,即水平近似线(75°≤θ≤105°)和垂直近似线(0°≤θ≤15°或θ≥165°),并按其中心位置分别归入顶部、底部、左侧、右侧四个队列中,且限定每条线的中心必须落在距边缘15%的区域内以防误判内部线条。
4. 计算外框线角点并映射回原图:每侧取坐标最靠外的一条:顶边取y最小值、底边取y最大值、左边取x最小值、右边取x最大值。四条线两两求交点得出外框的四个角点。由于所有计算均在缩图上进行,最后每个角点的坐标乘以2即可恢复至原始图像的像素精度。
5. 原图透视矫正:基于已计算出的外框角点坐标,使用cv2.getPerspective Transform()配合cv2.warpPerspective()以INTER_LANCZOS4(Lanczos高质量插值)执行透视变换,输出的图像均被统一矫正到6000×4750像素(宽高比为48:38)的标准矩形[4],同时输出边框角点识别结果标注,以便调试参数。
图4 步骤二示意图
步骤三:裁剪外框
经过第二步的透视变换后,每张图的黑色粗线外框已经变成了标准的矩形式样。但地图的内框,即地图上真正有等高线、地名和符号的区域,才是最终需要拼接的区域。为了便于进一步识别内框并矫正,需进行第二次裁剪。
策略与步骤一相同,再次利用Pillow(PIL)库的Image.crop()对所有地图执行简单的等距裁切,但改为对四边各裁去60像素。这一数值保证了能切掉外框,同时留出内框外一定空白区域,避免伤及任何有效地图区域(图5)。
图5 步骤三示意图
步骤四:识别内框线并矫正透视投影
通过前三步处理,每张地图获得一张内框区域外仍有一圈较窄空白区域的矩形图像。为了将地图沿内框进行裁切并矫正以实施最后的拼接,再次对其执行步骤二中的相同处理流程,输出的图像均被统一矫正到5750×4500像素(宽高比为46:36)的标准矩形。此时每张图像的边缘即为内框线所在位置,已具备下一步整体成图拼接的条件(图6)。
图6 步骤四示意图
理论上经过步骤二处理后,由于内外框间的宽窄相对统一(约125像素),通过统一裁剪图像四边该像素值后也可得到仅留存内框区域的图像。但考虑到制图与扫描等误差,本文仍采取裁去外框后再次用霍夫变换识别内框线并矫正透视投影,以保证准确性。
步骤五:成图拼接
最后一步将263幅完成裁切矫正的小图拼接成一张完整的大图。不同的场景对这一成果有不同的使用需求:一是仅需快速高效地对整幅拼接的地图进行简单的缩放、平移浏览;二是需要这幅图可加载到QGIS中作为标准地理图层进行进一步的地理配准或叠加分析。前者面向简单流畅的视觉浏览体验,后者面向专业的GIS工作流。
1. 面向网页浏览器:Deep Zoom瓦片金字塔
如果直接把263张图拼成一张约115000×99000像素的大图,则无法高效地在普通设备中流畅浏览。而DeepZoom Image(DZI,见小贴士4)格式提供了一种按需渐进式加载的方案。
小贴士4 DeepZoom (DZI)
通俗理解: 像“地图APP”一样看图的技术。
作用: 拼接后的全图像素极高(几亿像素),普通电脑无法直接打开。DeepZoom技术把大图切成无数个小方块(瓦片),你放大看哪里,它就只加载哪里的细节。这样即使在手机上,也能丝滑地缩放、平移查看这幅巨大的百年老地图。
根据地图索引图表(图7),整幅地图分为20列22行,以右上角为第1列第1行,将每张地图的行列位置保存在其文件名中,如01_06_xxx.jpg则表示其位于第1列第6行。将每张地图放在正确位置上构建虚拟画布(20列×22行=115000×99000像素),然后使用PyVips的dzsave生成Deep Zoom Image(DZI)瓦片金字塔,每一级包含不同分辨率的方块图片(瓦片)。在浏览器里放大到第N层时,只需加载当时视图范围内的瓦片。配合一个极简的HTML页面,即可支持在网页浏览器中通过OpenSeadragon JS库实现流畅缩放、任意平移的本地浏览体验(图8)。
图7 地图索引图表
图8 Web端成图浏览效果示意图
2. 面向QGIS:JGW文件+VRT拼接
对于需要在GIS系统中叠加分析的研究者而言,需为每张地图生成一个对应的JGW世界文件。在QGIS中加载这批裁切好的地图文件时,使用GDAL的gdalbuildvrt工具将263个带世界文件的地图合成为一个VRT虚拟栅格(见小贴士3),从而完成整图拼接。
本文仅在QGIS中实现了地图拼接。按照《四川省志·测绘志》原文描述,这批民国地形图本身没有经纬度或公里网参照系统,它们仅能表达相对位置关系,无法确定绝对坐标体系下的地理位置。如果想要进一步利用现代GIS平台进行叠加与时空动态分析,则未来必须完成下一步的地理配准工作:需要寻找可靠的古今对照控制点(GCPs),然后将图像拉伸扭曲至与实际地球曲面相匹配。同时考虑到地图本身的系统性误差和失真问题,可能还需要引入多项校准手段(图9)。
图9 GIS端成图浏览效果示意图
三、结语
本文介绍了一套基于Python脚本的自动化拼接方案,详细阐述了脚本原理与逻辑,将263幅四川省民国时期1:10万地图转化为可整体交互浏览的统一数据集。通过五步处理流程,实现历史地图从碎片化展示到整体化浏览,并提供Web端和GIS端两种输出形式。鉴于原图并无经纬度可参考,下一步仍需基于控制点进行地理配准,从而实现古今图层叠合对比分析等深层次应用。当前尚未进行配准的拼接地图,仍然对巴蜀地区历史地理、古道线路研究、考古调查、地名研究等具有积极的参考价值。
注释:
[1] 四川省地方志编纂委员会:《四川省志·测绘志》,成都地图出版社,1997年,第319-320页。
[2] 本文处理的地图扫描图像均来源于网络,仅用于学习研究。
[3] 李爽:《基于OpenCV与ArcPy的民国大比例尺地形图批量配准方法——以汉珍数位民国时期五万分一地图集为例》,《历史地理研究》2024年第2期,第109-122页。
[4] 原地图图幅宽高尺寸为46厘米×36厘米(内框),根据比例关系推算外框尺寸应为48厘米×38厘米,故考虑到原图分辨率尺寸,将外框裁切校正后的分辨率设为6000×4750像素(宽高比为48:38),同理步骤四中内框裁切校正后的分辨率设为5750×4500像素(宽高比为46:36)。
作者简介:孙锟,重庆市文物考古研究院建筑遗产研究所专业技术人员。