利用Python提取三維屬性數(shù)據(jù)的路線剖面地形線繪制

潘黎黎，林龍斌，李細(xì)光

(1.廣東省珠海工程勘察院，廣東 珠海 519002；2.山西省第三地質(zhì)工程勘察院，山西 晉中 030620；3.廣西壯族自治區(qū)地震局，廣西 南寧 530022)

0 引言

近年來，數(shù)字地質(zhì)填圖技術(shù)在地質(zhì)行業(yè)，特別是區(qū)域地質(zhì)調(diào)查填圖、城市地質(zhì)調(diào)查工作中發(fā)揮著重要作用[1-5]。國內(nèi)地礦行業(yè)主要使用數(shù)字地質(zhì)調(diào)查系統(tǒng)DGSS[6-9]進(jìn)行填圖等相關(guān)工作。在城市地質(zhì)調(diào)查三維基礎(chǔ)地質(zhì)結(jié)構(gòu)調(diào)查[10-11]野外工作中利用DGSS子系統(tǒng)數(shù)字地質(zhì)填圖系統(tǒng)RGMap[12-14](移動版)采集路線PRB(地質(zhì)調(diào)查點、點間路線、地質(zhì)界線)數(shù)據(jù)[15-17]，利用RGMap(PC版)對PRB野外數(shù)據(jù)進(jìn)行整理、成圖和管理等。移動版RGMap除能記錄地質(zhì)調(diào)查點上的描述信息，也能記錄調(diào)查點的位置信息，如經(jīng)緯度、投影平面XY坐標(biāo)、高程等信息。

通常，為整體了解路線上地層、接觸關(guān)系及構(gòu)造特征，建立宏觀認(rèn)識，需要繪制路線信手地質(zhì)剖面圖。傳統(tǒng)地質(zhì)調(diào)查中，路線信手地質(zhì)剖面圖多是在野外與調(diào)查同步進(jìn)行繪制，而在數(shù)字地質(zhì)填圖工作中往往利用數(shù)字地質(zhì)填圖系統(tǒng)繪制。RGMap系統(tǒng)(PC版)可以利用數(shù)字化地形圖自動生成路線信手剖面框架。但通常在城市地質(zhì)調(diào)查工作前期僅收集到大比例尺掃描版地形圖，而沒有數(shù)字化地形圖，這使得RGMap系統(tǒng)無法自動生成路線信手剖面框架圖中的地形線，而需從掃描版地形圖上讀取路線高程數(shù)據(jù)來繪制信手剖面地形線，極大地增加了工作量。因此，快速獲取路線地形數(shù)據(jù)是繪制路線信手剖面首先需要解決的問題。隨著人工智能、大數(shù)據(jù)等的快速發(fā)展，Python編程語言依靠其語言簡單、入門快、功能強大、跨平臺等優(yōu)勢，越來越受到大眾的喜愛。本文擬介紹利用Python腳本批量獲取地形線上地質(zhì)點橫縱坐標(biāo)，并結(jié)合Section軟件自動繪制路線信手剖面地形線的原理和方法等。

1 原理分析

信手剖面地形線常用展開法和投影法[18]繪制，RGMap系統(tǒng)默認(rèn)的繪制方法為展開法，因此，本文也僅討論展開法繪制的地形線。利用展開法繪制路線信手剖面地形線的一般步驟(圖1)：① 確定信手剖面地形線方向和起止點；② 確定地形線比例尺；③ 確定各點間的水平距離和各點的高程。為完整的反映路線上地質(zhì)概況，將路線首、尾地質(zhì)點作為信手剖面地形線的起止點，方向一般由起點指向終點，若路線整體方向變化較大，可在剖面上標(biāo)注多個方向。本文以比例尺為1∶1萬的三維基礎(chǔ)地質(zhì)結(jié)構(gòu)調(diào)查為例，設(shè)定信手剖面地形線比例尺為1∶1萬，繪制時先按設(shè)定比例尺將第一個點繪制在方格紙上，計算第二點與第一點的水平距離和高程，確定第2點的位置，以此類推，在方格紙上按順序先后繪出所有地質(zhì)點，最后將各地質(zhì)點用平滑曲線依次連接起來，即為路線信手剖面地形線。因此，繪制地形線最關(guān)鍵的是確定相鄰地質(zhì)點間的水平距離和各點高程。移動版RGMap系統(tǒng)記錄了地質(zhì)調(diào)查點的平面投影x、y坐標(biāo)和高程值，在PC版RGMap系統(tǒng)中可以通過GPOINT.WT文件查看地質(zhì)調(diào)查點的位置信息，利用Python腳本提取、處理地質(zhì)點位置信息依次獲取相鄰地質(zhì)點間的水平距離和高程，結(jié)合Section軟件的批量點數(shù)據(jù)導(dǎo)入功能，即可自動繪制信手剖面地形線。

2 數(shù)據(jù)準(zhǔn)備

PC版RGMap系統(tǒng)中GPOINT.WT文件記錄了野外地質(zhì)調(diào)查點的所有信息，并按屬性字段分別存儲，主要有ID(系統(tǒng)自增點記錄)、ROUTECODE(自定義路線號)、GEOPOINT(自定義點號)、LONGITUDE(經(jīng)度)、LATITUDE(緯度)、ALTITUDE(高程)、XX(投影平面X坐標(biāo))、YY(投影平面Y坐標(biāo))等(圖2a)等屬性字段。為快速提取地質(zhì)點的ALTITUDE、XX、YY屬性內(nèi)容，計算各點在地形線上的橫縱坐標(biāo)，繪制信手剖面地形線，先利用RGMap系統(tǒng)中“屬性數(shù)據(jù)導(dǎo)出Excel”的功能將這些屬性內(nèi)容導(dǎo)出為Excel文件。在RGMap系統(tǒng)中野外手圖界面將GPOINT.WT文件設(shè)置為編輯狀態(tài)，點擊菜單欄“工具”—“圖層屬性導(dǎo)出Excel”(圖2b)，在彈出的窗口“處理類型”默認(rèn)選擇點文件(圖2c)，“是否導(dǎo)出PRB描述信息”的窗口選擇“否”，系統(tǒng)自動調(diào)用MS-Office等工具打開導(dǎo)出的屬性數(shù)據(jù)(圖2d)，將其保存為Excel文件。

3 數(shù)據(jù)提取

3.1 導(dǎo)入第三方模塊

利用Python腳本從Excel文件中提取地質(zhì)調(diào)查點的高程和平面投影XX、YY坐標(biāo)并計算相鄰地質(zhì)點之間的水平距離，首先導(dǎo)入Python腳本處理數(shù)據(jù)所需的pandas、openpyxl、math、os、tkinter、subprocess等模塊。其中pandas模塊中的read_excel和ExcelWriter用于向Excel文件中讀取和寫入數(shù)據(jù)，openpyxl模塊用于處理“.xlsx”格式的Excel文件，math模塊中的sqrt用于計算兩點間的水平距離，os模塊中的path用于獲取Excel文件路徑相關(guān)參數(shù)，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

結(jié)合tkinter模塊filedialog函數(shù)返回的源Excel文件(Input_file)的絕對路徑，利用os模塊的path將結(jié)果文件(Output_file)保存于Input_file目錄下，便于查看和管理。用pandas模塊read_excel函數(shù)打開并讀取Input_file的數(shù)據(jù)并保存到內(nèi)存中。代碼如下：

root = Tk()

root.withdraw()

Input_file = filedialog.askopenfilename()

圖2 (a)地質(zhì)點屬性字段 (b)處理文件類型 (c)圖層屬性導(dǎo)出Excel (d)Excel屬性數(shù)據(jù)

(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 橫坐標(biāo)的提取

定義一個list列表保存地質(zhì)調(diào)查點在信手剖面上的橫坐標(biāo)，并且設(shè)定第一個地質(zhì)調(diào)查點的橫坐標(biāo)為0(即路線信手剖面的起點)。依次計算兩相鄰兩地質(zhì)調(diào)查點之間的水平距離，假定S1為第2點與起點之間的水平距離，S2為第3點與第2點之間的水平距離，以此類推，Sn為第n+1點與第n點之間的水平距離，那么第n+1點的橫坐標(biāo)應(yīng)等于S1，S2，S3，...，Sn之和，依此原理計算各地質(zhì)點的在信手剖面上的橫坐標(biāo)值。式(1)為平面投影直角坐標(biāo)系下兩點之間的水平距離計算公式，式(2)為地質(zhì)點在信手剖面上的橫坐標(biāo)值計算公式。

(1)

(2)

式(1)中，Sn為第n+1地質(zhì)點與第n地質(zhì)點之間的水平距離，Xn和Yn分別為第n個地質(zhì)點的平面投影坐標(biāo)XX和YY；式(2)中Tn為第n地質(zhì)點在信手剖面上的橫坐標(biāo)，默認(rèn)T1為0。Python中可利用math模塊內(nèi)置的sqrt計算平面直角坐標(biāo)系下兩點間的水平距離，并進(jìn)行循環(huán)累加求和計算各地質(zhì)點的橫坐標(biāo)。野外RGMap系統(tǒng)采集的XX、YY坐標(biāo)及高程ALTITUDE單位均為“m”，而路線信手剖面圖的比例尺設(shè)定為1∶1萬，RGMap系統(tǒng)中最小單位為“mm”，所以需要將以m為單位的Tn數(shù)據(jù)轉(zhuǎn)換成1∶1萬比例尺下以mm為單位的數(shù)據(jù)，即將Tn值除以10即得以mm為單位的橫坐標(biāo)x。代碼如下：

#計算地質(zhì)點的橫坐標(biāo)

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 縱坐標(biāo)的提取

通常信手剖面圖橫縱比例尺保持一致，路線信手剖面上地質(zhì)點的縱坐標(biāo)為地質(zhì)點的高程，但需對RGMap系統(tǒng)采集的高程數(shù)據(jù)ALTITUDE進(jìn)行變換，與XX、YY坐標(biāo)變換方法一樣，將ALTITUDE數(shù)據(jù)除以10即為1∶1萬比例尺下以mm為單位的數(shù)據(jù)。一幅完整的路線信手剖面圖需將比例尺及圖例擺放于圖下方，為便于制圖，將地形線整體向上平移50 mm，所以需在轉(zhuǎn)換后的高程數(shù)據(jù)(縱坐標(biāo))上加50mm得到最終縱坐標(biāo)值y。代碼如下：

#計算地質(zhì)點的眾坐標(biāo)

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 點編號及數(shù)據(jù)保存

路線信手剖面上要求標(biāo)注地質(zhì)點號，所以需從屬性字段GEOPOINT(點號)提取地質(zhì)點號GEOPOINT(圖2d)，將計算的橫縱坐標(biāo)x、y與地質(zhì)點號ID一并寫入Output_file，利用subprocess模塊內(nèi)置Popen自動打開Output_file，代碼見圖3a。Output_file包含三列數(shù)據(jù)，即x、y和ID(圖3b)，分別代表地質(zhì)點的橫、縱坐標(biāo)和地質(zhì)點號。

圖3 (a)代碼段 (b)地質(zhì)點的橫縱坐標(biāo)和編號

4 地形線繪制

Section V4.6是在MapGIS6.7基礎(chǔ)上，利用Microsoft VC++6.0語言編寫的地質(zhì)圖件制作軟件?；贛apGIS 6.7輸入編輯子系統(tǒng)強大的點、線、區(qū)編輯能力，添加專業(yè)的地質(zhì)圖件制作工具，大大提高了地質(zhì)圖件的制作效率。

4.1 繪制地形線

首先創(chuàng)建路線信手剖面圖點、線、區(qū)文件。在PC版RGMap系統(tǒng)野外手圖界面，選擇菜單欄上“地質(zhì)填圖數(shù)據(jù)操作”—“信手剖面生成與瀏覽”—“顯示”(圖4a)，即可自動在相應(yīng)野外路線的“素描圖”文件夾里生成信手剖面圖點、線、區(qū)和對應(yīng)的工程文件。打開Section軟件，找到“素描圖”文件夾，加載自動生成的信手剖面工程文件，即可對信手剖面圖點、線、區(qū)文件進(jìn)行編輯。打開前述Output_file文件，點擊菜單欄上“1輔助工具”—“表格數(shù)據(jù)投影”—“全部數(shù)據(jù)投影(EXCEL)”(圖4b)。彈出“數(shù)據(jù)投影”窗口，去掉“線閉合”前面的勾(圖4c)，分別設(shè)置“文字圖元參數(shù)”、“子圖圖元參數(shù)”(圖4d、4e)和“線圖元參數(shù)”，為使信手剖面地形線圓滑，可以將線型設(shè)置為曲線。文字、子圖、線圖元等參數(shù)設(shè)置完成后，點擊“確認(rèn)”，即可自動生成地形線，并在地形線上標(biāo)注了地質(zhì)點號(圖5a)。

4.2 調(diào)整點標(biāo)注

自動標(biāo)注的文字和地質(zhì)點圖元擺放并不美觀，為便于整飾圖面，在設(shè)置文字和子圖參數(shù)時，需將文字和子圖的圖層號設(shè)置為不同的數(shù)字(圖4d、4e)，地形線生成后，點擊菜單欄上“A圖層”—“改層開關(guān)”—“改點P”將地質(zhì)點圖層關(guān)閉(圖4f)，移動地質(zhì)點編號位置，使圖面更美觀(圖5b)，至此地形線繪制完畢，補充完善信手剖面圖地層、巖體界線和花紋、構(gòu)造、圖例、比例尺、方位等要素即可繪制完成信手剖面圖。

5 結(jié)語

路線信手剖面圖作為整體了解路線巖性、構(gòu)造特征的重要圖件，由于缺失數(shù)字化地形圖，無法在數(shù)字地質(zhì)填圖系統(tǒng)中自動生成信手剖面圖框架。經(jīng)研究，利用路線上地質(zhì)點的平面投影坐標(biāo)計算相鄰地質(zhì)點間的水平距離，結(jié)合地質(zhì)調(diào)查點的高程，可繪制出大致的路線信手剖面地形線。

圖4 (a)自動生成信手剖面 (b)表格數(shù)據(jù)投影 (c)數(shù)據(jù)投影

圖5 (a)自動生成的地形線和標(biāo)注 (b)調(diào)整后的地形線和標(biāo)注Fig.5 (a)Automatically created topographic line and annotation (b)Adjusted topographic line and annotation

1)利用Python腳本批量處理RGMap系統(tǒng)中地質(zhì)點屬性數(shù)據(jù)，提取地質(zhì)點的平面投影x、y坐標(biāo)和高程數(shù)據(jù)，采用平面坐標(biāo)系中兩點間的距離公式計算相鄰地質(zhì)點間的水平距離，并進(jìn)行累加求和，即可計算出各地質(zhì)點在設(shè)定剖面比例尺下的橫坐標(biāo)，而地質(zhì)點的縱坐標(biāo)為地質(zhì)點高程屬性數(shù)據(jù)在設(shè)定比例尺下的值。

2)利用Section軟件的表格數(shù)據(jù)投影功能，將Python腳本計算所得的各地質(zhì)點橫、縱坐標(biāo)和地質(zhì)點號文件導(dǎo)入Section軟件，同時設(shè)置地質(zhì)點標(biāo)注、地質(zhì)點子圖、線型等參數(shù)，即可自動繪制出路線信手剖面地形線。

3)利用Python提取、處理數(shù)據(jù)和Section批量導(dǎo)入數(shù)據(jù)繪制地形線能一定程度上提高工作效率，減少由于人為因素帶來的誤差，提高了數(shù)據(jù)的準(zhǔn)確度。

本文主要利用野外采集的數(shù)據(jù)進(jìn)行分析計算，獲得地質(zhì)點的橫縱坐標(biāo)。因此，通過計算獲取的橫縱坐標(biāo)數(shù)據(jù)精度主要依賴于野外采集的x、y坐標(biāo)及高程數(shù)據(jù)的精度。目前常用的野外數(shù)字填圖設(shè)備的GPS定位系統(tǒng)獲取的x、y坐標(biāo)的精度較高，基本可以滿足繪圖的需要，而通常獲取到的高程數(shù)據(jù)精度較低，這可能會影響到地形線的垂直精度，本文存在的不足之處在于直接使用了GPS采集的高程數(shù)據(jù)。因此，如需更準(zhǔn)確的高程數(shù)據(jù)，可以在RGMap系統(tǒng)中根據(jù)電子掃描地形圖對野外采集的地質(zhì)點高程數(shù)據(jù)進(jìn)行修正后，再利用python腳本進(jìn)行處理。Python語言依據(jù)其獨特的優(yōu)勢，不僅能批量處理地質(zhì)點數(shù)據(jù)，還可以應(yīng)用于城市地質(zhì)調(diào)查工作中的鉆孔資料整理分析、數(shù)據(jù)庫建設(shè)數(shù)據(jù)準(zhǔn)備等各階段，能極大的提高批量數(shù)據(jù)分析處理的工作效率，相信在以后的地質(zhì)工作中應(yīng)用將會越來越廣泛。