利用Python提取三维属性数据的路线剖面地形线绘制

潘黎黎，林龙斌，李细光

(1.广东省珠海工程勘察院，广东 珠海 519002；2.山西省第三地质工程勘察院，山西 晋中 030620；3.广西壮族自治区地震局，广西 南宁 530022)

0 引言

近年来，数字地质填图技术在地质行业，特别是区域地质调查填图、城市地质调查工作中发挥着重要作用[1-5]。国内地矿行业主要使用数字地质调查系统DGSS[6-9]进行填图等相关工作。在城市地质调查三维基础地质结构调查[10-11]野外工作中利用DGSS子系统数字地质填图系统RGMap[12-14](移动版)采集路线PRB(地质调查点、点间路线、地质界线)数据[15-17]，利用RGMap(PC版)对PRB野外数据进行整理、成图和管理等。移动版RGMap除能记录地质调查点上的描述信息，也能记录调查点的位置信息，如经纬度、投影平面XY坐标、高程等信息。

通常，为整体了解路线上地层、接触关系及构造特征，建立宏观认识，需要绘制路线信手地质剖面图。传统地质调查中，路线信手地质剖面图多是在野外与调查同步进行绘制，而在数字地质填图工作中往往利用数字地质填图系统绘制。RGMap系统(PC版)可以利用数字化地形图自动生成路线信手剖面框架。但通常在城市地质调查工作前期仅收集到大比例尺扫描版地形图，而没有数字化地形图，这使得RGMap系统无法自动生成路线信手剖面框架图中的地形线，而需从扫描版地形图上读取路线高程数据来绘制信手剖面地形线，极大地增加了工作量。因此，快速获取路线地形数据是绘制路线信手剖面首先需要解决的问题。随着人工智能、大数据等的快速发展，Python编程语言依靠其语言简单、入门快、功能强大、跨平台等优势，越来越受到大众的喜爱。本文拟介绍利用Python脚本批量获取地形线上地质点横纵坐标，并结合Section软件自动绘制路线信手剖面地形线的原理和方法等。

1 原理分析

信手剖面地形线常用展开法和投影法[18]绘制，RGMap系统默认的绘制方法为展开法，因此，本文也仅讨论展开法绘制的地形线。利用展开法绘制路线信手剖面地形线的一般步骤(图1)：① 确定信手剖面地形线方向和起止点；② 确定地形线比例尺；③ 确定各点间的水平距离和各点的高程。为完整的反映路线上地质概况，将路线首、尾地质点作为信手剖面地形线的起止点，方向一般由起点指向终点，若路线整体方向变化较大，可在剖面上标注多个方向。本文以比例尺为1∶1万的三维基础地质结构调查为例，设定信手剖面地形线比例尺为1∶1万，绘制时先按设定比例尺将第一个点绘制在方格纸上，计算第二点与第一点的水平距离和高程，确定第2点的位置，以此类推，在方格纸上按顺序先后绘出所有地质点，最后将各地质点用平滑曲线依次连接起来，即为路线信手剖面地形线。因此，绘制地形线最关键的是确定相邻地质点间的水平距离和各点高程。移动版RGMap系统记录了地质调查点的平面投影x、y坐标和高程值，在PC版RGMap系统中可以通过GPOINT.WT文件查看地质调查点的位置信息，利用Python脚本提取、处理地质点位置信息依次获取相邻地质点间的水平距离和高程，结合Section软件的批量点数据导入功能，即可自动绘制信手剖面地形线。

2 数据准备

PC版RGMap系统中GPOINT.WT文件记录了野外地质调查点的所有信息，并按属性字段分别存储，主要有ID(系统自增点记录)、ROUTECODE(自定义路线号)、GEOPOINT(自定义点号)、LONGITUDE(经度)、LATITUDE(纬度)、ALTITUDE(高程)、XX(投影平面X坐标)、YY(投影平面Y坐标)等(图2a)等属性字段。为快速提取地质点的ALTITUDE、XX、YY属性内容，计算各点在地形线上的横纵坐标，绘制信手剖面地形线，先利用RGMap系统中“属性数据导出Excel”的功能将这些属性内容导出为Excel文件。在RGMap系统中野外手图界面将GPOINT.WT文件设置为编辑状态，点击菜单栏“工具”—“图层属性导出Excel”(图2b)，在弹出的窗口“处理类型”默认选择点文件(图2c)，“是否导出PRB描述信息”的窗口选择“否”，系统自动调用MS-Office等工具打开导出的属性数据(图2d)，将其保存为Excel文件。

3 数据提取

3.1 导入第三方模块

利用Python脚本从Excel文件中提取地质调查点的高程和平面投影XX、YY坐标并计算相邻地质点之间的水平距离，首先导入Python脚本处理数据所需的pandas、openpyxl、math、os、tkinter、subprocess等模块。其中pandas模块中的read_excel和ExcelWriter用于向Excel文件中读取和写入数据，openpyxl模块用于处理“.xlsx”格式的Excel文件，math模块中的sqrt用于计算两点间的水平距离，os模块中的path用于获取Excel文件路径相关参数，tkinter模块中的filedialog和Tk用于从磁盘选择要处理的Excel文件并返回文件的绝对路径。代码如下：

from pandas import read_excel，ExcelWriter

import openpyxl

from math import sqrt

from os import path

from tkinter import filedialog，Tk

from subprocess import Popen

结合tkinter模块filedialog函数返回的源Excel文件(Input_file)的绝对路径，利用os模块的path将结果文件(Output_file)保存于Input_file目录下，便于查看和管理。用pandas模块read_excel函数打开并读取Input_file的数据并保存到内存中。代码如下：

root = Tk()

root.withdraw()

Input_file = filedialog.askopenfilename()

图2 (a)地质点属性字段 (b)处理文件类型 (c)图层属性导出Excel (d)Excel属性数据

(c)Properties of map layer to excel file (d)Excel property data

if input_file：

Output_file = path.dirname(input_file) +'/地形线_'+ path.basename(input_file)

Worksheet_data = read_excel(input_file，'Sheet1'，index_col=None)

3.2 横坐标的提取

定义一个list列表保存地质调查点在信手剖面上的横坐标，并且设定第一个地质调查点的横坐标为0(即路线信手剖面的起点)。依次计算两相邻两地质调查点之间的水平距离，假定S1为第2点与起点之间的水平距离，S2为第3点与第2点之间的水平距离，以此类推，Sn为第n+1点与第n点之间的水平距离，那么第n+1点的横坐标应等于S1，S2，S3，...，Sn之和，依此原理计算各地质点的在信手剖面上的横坐标值。式(1)为平面投影直角坐标系下两点之间的水平距离计算公式，式(2)为地质点在信手剖面上的横坐标值计算公式。

(1)

(2)

式(1)中，Sn为第n+1地质点与第n地质点之间的水平距离，Xn和Yn分别为第n个地质点的平面投影坐标XX和YY；式(2)中Tn为第n地质点在信手剖面上的横坐标，默认T1为0。Python中可利用math模块内置的sqrt计算平面直角坐标系下两点间的水平距离，并进行循环累加求和计算各地质点的横坐标。野外RGMap系统采集的XX、YY坐标及高程ALTITUDE单位均为“m”，而路线信手剖面图的比例尺设定为1∶1万，RGMap系统中最小单位为“mm”，所以需要将以m为单位的Tn数据转换成1∶1万比例尺下以mm为单位的数据，即将Tn值除以10即得以mm为单位的横坐标x。代码如下：

#计算地质点的横坐标

List_x = [0，]

sum = 0

for i in range(1，len(Worksheet_data['XX']))：

sum+=sqrt((Worksheet_data['XX'][i]-Worksheet_data['XX'][i-1])**2+(Worksheet_data['YY'][i]-Worksheet_data['YY'][i-1])**2)/10

List_x.append(sum)

Worksheet_data['y'] = List_x

3.3 纵坐标的提取

通常信手剖面图横纵比例尺保持一致，路线信手剖面上地质点的纵坐标为地质点的高程，但需对RGMap系统采集的高程数据ALTITUDE进行变换，与XX、YY坐标变换方法一样，将ALTITUDE数据除以10即为1∶1万比例尺下以mm为单位的数据。一幅完整的路线信手剖面图需将比例尺及图例摆放于图下方，为便于制图，将地形线整体向上平移50 mm，所以需在转换后的高程数据(纵坐标)上加50mm得到最终纵坐标值y。代码如下：

#计算地质点的众坐标

List_y = []

for i in range(len(Worksheet_data['ALTITUDE']))：

item = Worksheet_data['ALTITUDE'][i]/10 + 50

List_y.append(item)

Worksheet_data['x'] = List_y

3.4 点编号及数据保存

路线信手剖面上要求标注地质点号，所以需从属性字段GEOPOINT(点号)提取地质点号GEOPOINT(图2d)，将计算的横纵坐标x、y与地质点号ID一并写入Output_file，利用subprocess模块内置Popen自动打开Output_file，代码见图3a。Output_file包含三列数据，即x、y和ID(图3b)，分别代表地质点的横、纵坐标和地质点号。

图3 (a)代码段 (b)地质点的横纵坐标和编号

4 地形线绘制

Section V4.6是在MapGIS6.7基础上，利用Microsoft VC++6.0语言编写的地质图件制作软件。基于MapGIS 6.7输入编辑子系统强大的点、线、区编辑能力，添加专业的地质图件制作工具，大大提高了地质图件的制作效率。

4.1 绘制地形线

首先创建路线信手剖面图点、线、区文件。在PC版RGMap系统野外手图界面，选择菜单栏上“地质填图数据操作”—“信手剖面生成与浏览”—“显示”(图4a)，即可自动在相应野外路线的“素描图”文件夹里生成信手剖面图点、线、区和对应的工程文件。打开Section软件，找到“素描图”文件夹，加载自动生成的信手剖面工程文件，即可对信手剖面图点、线、区文件进行编辑。打开前述Output_file文件，点击菜单栏上“1辅助工具”—“表格数据投影”—“全部数据投影(EXCEL)”(图4b)。弹出“数据投影”窗口，去掉“线闭合”前面的勾(图4c)，分别设置“文字图元参数”、“子图图元参数”(图4d、4e)和“线图元参数”，为使信手剖面地形线圆滑，可以将线型设置为曲线。文字、子图、线图元等参数设置完成后，点击“确认”，即可自动生成地形线，并在地形线上标注了地质点号(图5a)。

4.2 调整点标注

自动标注的文字和地质点图元摆放并不美观，为便于整饰图面，在设置文字和子图参数时，需将文字和子图的图层号设置为不同的数字(图4d、4e)，地形线生成后，点击菜单栏上“A图层”—“改层开关”—“改点P”将地质点图层关闭(图4f)，移动地质点编号位置，使图面更美观(图5b)，至此地形线绘制完毕，补充完善信手剖面图地层、岩体界线和花纹、构造、图例、比例尺、方位等要素即可绘制完成信手剖面图。

5 结语

路线信手剖面图作为整体了解路线岩性、构造特征的重要图件，由于缺失数字化地形图，无法在数字地质填图系统中自动生成信手剖面图框架。经研究，利用路线上地质点的平面投影坐标计算相邻地质点间的水平距离，结合地质调查点的高程，可绘制出大致的路线信手剖面地形线。

图4 (a)自动生成信手剖面 (b)表格数据投影 (c)数据投影

图5 (a)自动生成的地形线和标注 (b)调整后的地形线和标注Fig.5 (a)Automatically created topographic line and annotation (b)Adjusted topographic line and annotation

1)利用Python脚本批量处理RGMap系统中地质点属性数据，提取地质点的平面投影x、y坐标和高程数据，采用平面坐标系中两点间的距离公式计算相邻地质点间的水平距离，并进行累加求和，即可计算出各地质点在设定剖面比例尺下的横坐标，而地质点的纵坐标为地质点高程属性数据在设定比例尺下的值。

2)利用Section软件的表格数据投影功能，将Python脚本计算所得的各地质点横、纵坐标和地质点号文件导入Section软件，同时设置地质点标注、地质点子图、线型等参数，即可自动绘制出路线信手剖面地形线。

3)利用Python提取、处理数据和Section批量导入数据绘制地形线能一定程度上提高工作效率，减少由于人为因素带来的误差，提高了数据的准确度。

本文主要利用野外采集的数据进行分析计算，获得地质点的横纵坐标。因此，通过计算获取的横纵坐标数据精度主要依赖于野外采集的x、y坐标及高程数据的精度。目前常用的野外数字填图设备的GPS定位系统获取的x、y坐标的精度较高，基本可以满足绘图的需要，而通常获取到的高程数据精度较低，这可能会影响到地形线的垂直精度，本文存在的不足之处在于直接使用了GPS采集的高程数据。因此，如需更准确的高程数据，可以在RGMap系统中根据电子扫描地形图对野外采集的地质点高程数据进行修正后，再利用python脚本进行处理。Python语言依据其独特的优势，不仅能批量处理地质点数据，还可以应用于城市地质调查工作中的钻孔资料整理分析、数据库建设数据准备等各阶段，能极大的提高批量数据分析处理的工作效率，相信在以后的地质工作中应用将会越来越广泛。