張偉，邱崇濤，謝明宏，趙叢

(1.核工業(yè)航測遙感中心，河北 石家莊 050002；2.中核集團(tuán) 鈾資源地球物理勘查技術(shù)中心(重點(diǎn)實(shí)驗(yàn)室)，河北 石家莊 050002)

0 引言

在地質(zhì)、物探工作中，自始至終都會應(yīng)用不同形式的專題地圖來直觀地表現(xiàn)其研究成果[1]。為獲取最佳的地質(zhì)效果，技術(shù)人員需不斷地對地質(zhì)情況深入研究，反復(fù)斟酌，經(jīng)多次易稿后才能獲取最終滿意的成果，而伴隨其中的專題圖件也在不斷的修改中完善。

在物探工作中，地質(zhì)推斷解釋圖是反映物探成果的必備專題圖件之一?？傮w而言，地質(zhì)推斷解釋圖中的圖元要素并非十分復(fù)雜；但由于修改次數(shù)較多，又受到專業(yè)制圖軟件功能的限制，工作效率難以滿足快節(jié)奏項(xiàng)目需求。以目前國內(nèi)應(yīng)用范圍最廣的國產(chǎn)軟件MapGIS為例，雖有強(qiáng)大編輯功能，但憑手工繪圖方式，其效率依然難有明顯提升。如電法勘探中的反演電阻率斷面的地質(zhì)推斷解釋圖，不整合線是一種常見的地質(zhì)現(xiàn)象，圖中使用波浪線表示，而在MapGIS軟件中波浪線則由一個(gè)特殊線型呈現(xiàn)，實(shí)際坐標(biāo)位置處于曲線的中軸線上；若直接成區(qū)，線與區(qū)間的分界面無法滿足制圖要求，仍需徒手繪制或借助于第三方軟件生成，操作繁瑣且費(fèi)時(shí)費(fèi)力。筆者針對該類項(xiàng)目制圖特點(diǎn)和需求，利用Python語言及其強(qiáng)大的第三方庫，設(shè)計(jì)了自動(dòng)制圖程序，主要功能包括波浪線(不整合界線)、自動(dòng)剪斷線、微短線識別與刪除以及多邊形化成區(qū)等功能，目的是提高工作效率，減少制圖過程中不必要的錯(cuò)誤發(fā)生。

1 Python語言簡介

Python語言是一個(gè)純凈、簡潔明了的編程語言，特別適合于初學(xué)者，并且還可幫助專業(yè)人員快速地解決復(fù)雜的問題，它可以輕松快速地對地理空間數(shù)據(jù)進(jìn)行可視化處理[2]。Python具有強(qiáng)大的科學(xué)計(jì)算、繪圖功能，同時(shí)在地圖繪制、地理空間數(shù)據(jù)的處理與轉(zhuǎn)換方面具有豐富的類庫[3]。本文主要涉及了NumPy(N維數(shù)組和矩陣計(jì)算)、Shapely(二維圖形處理庫)、Pandas(數(shù)據(jù)清洗庫)和Scipy(數(shù)學(xué)工具包和算法庫)。另外，在程序調(diào)試中，還使用了Geopandas(地圖繪制)和GeoJson(空間數(shù)據(jù)格式)庫，用于預(yù)覽繪制效果。

Shapely為Python中二維圖形軟件包，基于廣為流傳的GEOS庫中的函數(shù)，進(jìn)行幾何體的操作和分析[4]。本文主要涉及該庫中的幾何體裁剪、多邊形化、節(jié)點(diǎn)抽稀和幾何體空間位置判定等函數(shù)[5]。

Scipy是構(gòu)建在Python擴(kuò)展庫Numpy上的數(shù)學(xué)算法和便捷函數(shù)的集合。增加了強(qiáng)大的交互式Python任務(wù)，向用戶提供了高級別的指令或類用于數(shù)據(jù)操作和可視化[6-7]。本文主要利用該庫提供的插值(interpolate)模塊。

2 程序設(shè)計(jì)及實(shí)現(xiàn)

2.1 總體思路

一般情況下，當(dāng)物探工程師在繪制好地質(zhì)解釋圖后，交付給制圖人員進(jìn)行線編輯、顏色填充、圖面整飾等工作。利用MapGIS的制圖流程大體包括自動(dòng)裁剪線、線轉(zhuǎn)弧段、拓?fù)渲亟?、修改地?巖體)顏色、手工去除微短線等。筆者基于Python語言及其Shapely、NumPy、Scipy和Pandas等擴(kuò)展庫，參照上述MapGIS繪制流程，同時(shí)也補(bǔ)充了不整合線(波浪線)的繪制，達(dá)到了自動(dòng)制圖工作的目的。程序?qū)崿F(xiàn)的基本流程見圖1。

圖1 程序自動(dòng)化繪制流程Fig.1 The flow of automatic drawingby programming

2.2 程序的實(shí)現(xiàn)

MapGIS圖形按點(diǎn)、線、區(qū)分類存儲，與之對應(yīng)的文件分別為點(diǎn)文件、線文件、區(qū)文件，明碼格式文件名后綴分別為.wat、.wal和.wap[8]，詳細(xì)格式可查閱MapGIS相關(guān)文獻(xiàn)[9]。為脫離MapGIS運(yùn)行環(huán)境，擴(kuò)展程序的適用性，此次調(diào)用了MapGIS的明碼格式。

2.2.1 不整合界線的繪制

在上、下地層之間存有一個(gè)沉積間斷面，叫不整合面，在地面的出露線稱作不整合線[10]。在斷面圖中，地層不整合界線一般用波浪線表示。由于MapGIS程序中提供波浪線線型無法滿足區(qū)域填充的要求，首先通過程序生成一條適合區(qū)域填充的波浪線，從圖形學(xué)上，可以一條曲線或折線為軸心，由多個(gè)連續(xù)的正弦曲線或余弦曲線繪制成波浪線。繪制過程主要包括曲線插值平滑曲線，等路徑插值提取節(jié)點(diǎn)坐標(biāo)，根據(jù)節(jié)點(diǎn)位置繪制正弦曲線或余弦曲線形成波浪線，以及節(jié)點(diǎn)抽稀減輕運(yùn)算負(fù)荷等4部分(圖2)。

圖2 波浪線繪制流程Fig.2 Drawing of wave line

1) 曲線插值

曲線插值的目的是光滑曲線，使圖面更為美觀，同時(shí)又要避免曲線過度變形。在Scipy庫中，scipy.interpolate.splprep是適合于多維曲線的B樣條插值法，經(jīng)多次試驗(yàn)，s(光滑系數(shù))取0，k(曲線擬合度)取2較為合適。Python實(shí)現(xiàn)代碼如下：

from scipy.interpolate import splprep, splev

tck, u = splprep([xr, yr], s=0, k=2)

u = np.arange(0, 1., 0.0002)

#獲取XY坐標(biāo)值

xySet = splev(u, tck)

2) 等路徑提取節(jié)點(diǎn)坐標(biāo)

在Matlab中提供了沿路徑插值法(interpolate along the path)，而Scipy并未提供類似的函數(shù)。在Shapely庫中interpolate函數(shù)功能是沿線性幾何體獲得指定距離的坐標(biāo)值，通過等距離迭代計(jì)算可獲得類似效果。

#獲得線的長度

myLine = LineString(xy0)

lineLen = myLine.length

# sizeUnit為單個(gè)正弦曲線的波長

segNum = int(lineLen/sizeUnit)+2

dataSet3=[]

for i in range(segNum):

#獲取等路徑的坐標(biāo)值，存入數(shù)組中

ip = myLine.interpolate(i*sizeUnit, normalized=False)

dataSet3.append(list(ip.coords)[0])

3) 繪制波浪線

為了防止正弦(余弦)曲線變形，將預(yù)先生成的正弦曲線或余弦曲線，通過旋轉(zhuǎn)、平移操作移至曲線的相應(yīng)位置，旋轉(zhuǎn)平移的代碼如下：

#空間幾何體的旋轉(zhuǎn)、平移函數(shù)

def rotateMoving(xunit, ang, yunit, x0, y0):

x1 = xunit * np.cos(ang) + yunit * np.sin(ang) + x0

y1 = yunit * np.cos(ang) - xunit * np.sin(ang) + y0

return x1, y1

為盡量減少曲線誤差，正弦曲線或余弦曲線的范圍各為0～π/2，交替繪制。此時(shí)，需要考慮節(jié)點(diǎn)的矢量方向，確保正弦(余弦)曲線銜接連續(xù)。代碼如下：

if i % 2 == 0: #偶數(shù)點(diǎn)

if x3[i] < x3[i + 1]:

#當(dāng)后一點(diǎn)坐標(biāo)值X大于前一點(diǎn)時(shí)，

xx1, yy1 = rotateMoving(xunit, ang, yunit0, x3[i], y3[i])

else:

xx1, yy1 = rotateMoving(xunit, ang, -1*yunit0, x3[i+1], y3[i+1])

#翻轉(zhuǎn)數(shù)組

xx1 = np.flipud(xx1)

yy1 = np.flipud(yy1)

else：#奇數(shù)點(diǎn)

(余弦曲線與正弦曲線平移類似，略)

#存儲節(jié)點(diǎn)位置

allXPoints.append(xx1a)

allYPoints.append(yy1a)

4) 節(jié)點(diǎn)抽稀

此時(shí)波浪線的節(jié)點(diǎn)十分致密，將會加大計(jì)算機(jī)運(yùn)行負(fù)荷，尤其是長剖面，線節(jié)點(diǎn)過多會導(dǎo)致MapGIS無法正常運(yùn)行。線節(jié)點(diǎn)的抽稀方法采用道格拉斯—普克算法(douglas-peucker algorithm)，在Shapely中提供了此功能(simplify函數(shù))。經(jīng)試驗(yàn)對比，抽稀系數(shù)為0.01較為適合。代碼如下：

xySet = linTem.simplify(0.01, preserve_topology=False)

2.2.2 自動(dòng)剪斷線與微短線識別

從數(shù)學(xué)上，一個(gè)多邊形(ploygon)由同一平面上閉合曲線構(gòu)成，填充后可形成區(qū)。因此，繪制圖件時(shí)，為形成一個(gè)閉合多邊形，構(gòu)成多邊形輪廓的曲線必須有“搭接”。此環(huán)節(jié)可使用shapely.opt庫中的split函數(shù)，通過迭代計(jì)算，進(jìn)行不同幾何體間相互分割，實(shí)現(xiàn)要求：① 被分割幾何體為單個(gè)線幾何體(linestring)；② 分割幾何體為被分割體外的其他所有線幾何體的集合(multilinestring)，確保分割幾何體必須一次性完成剪斷任務(wù)。關(guān)鍵代碼如下：

#導(dǎo)入shapely庫

from shapely.ops import polygonize, split

#遍歷所有曲線數(shù)組

for i in range (len(lineSets)):

#生成第i條線(被裁剪)

oneLine = LineString(lineSets [i])

#從線數(shù)組中刪除第i條線

others = list(np.delete(lineSets, i))

#生成除第i條線外的線集合

otherLines = MultiLineString(others)

#裁剪線

resultCollect = split(oneLine, otherLines)

#獲取線集合中線的數(shù)量

lineNum = len(resultCollect.geoms)

獲得的線裁剪結(jié)果(result collect)為多線集合，此集合內(nèi)幾何體排列順序與原始線坐標(biāo)位置一致，從而為判定或刪除(地層)線兩端外延部分奠定了條件。結(jié)合斷面圖特點(diǎn)，主要判定方法包括：① 判定特殊線型，如紅色的斷裂線和閉合曲線(地質(zhì)體、圖框)，此類線幾何體將不進(jìn)行任何刪除處理；斷裂線根據(jù)線屬性(顏色)判定，閉合曲線可根據(jù)首尾點(diǎn)的坐標(biāo)值確定。② 根據(jù)多線集合中幾何體的數(shù)量(lineNum≥3)判定線幾何體否已被裁剪，來決定是否刪除第一個(gè)和最后一個(gè)幾何體(線)。最后以MapGIS明碼格式保存線參數(shù)(節(jié)點(diǎn)坐標(biāo)值、顏色、寬度、線型等)。

2.2.3 折線成區(qū)和顏色填充

此過程相當(dāng)于MapGIS中的線轉(zhuǎn)弧段和拓?fù)涑蓞^(qū)。但MapGIS是隨機(jī)填充顏色，而本程序是根據(jù)地層、巖體的注釋點(diǎn)進(jìn)行顏色填充。

1) 折線轉(zhuǎn)多邊形

此處的折線轉(zhuǎn)多邊形相當(dāng)于MapGIS中的“線轉(zhuǎn)弧段”，Shapely提供了polygonize函數(shù)，代碼如下：

#折線的多邊形化，outLinSets為折線數(shù)組

polygonSet = polygonize(outLinSets)

# 提取多邊形的坐標(biāo)值

pgSet = list(polygonSet)

2) 讀取點(diǎn)文件

使用Pandas讀取wat文件更為方便，主要獲取注釋點(diǎn)(X，Y)坐標(biāo)值、字體顏色和標(biāo)注。當(dāng)然，也可以限制一些與顏色填充無關(guān)的注釋點(diǎn)(非強(qiáng)制)，如斷裂編號，為后續(xù)地層符號識別去除一些干擾因素。

#使用pandas讀取wat文件

df = pd.read_table(watFile, skiprows=2, sep=′,′,

header=None, usecols=[0, 1, 4, 13], encoding=′gbk′)

df.columns = [′x′, ′y′, ′text′, ′clr′]

#去除紅色斷裂注釋點(diǎn)(顏色為6)

df2 = df.loc[(df[′clr′] != 6), [′x′, ′y′, ′text′]]

3) 確定多邊形內(nèi)的注釋點(diǎn)

通過對多邊形和注釋點(diǎn)雙迭代運(yùn)算，利用Shapely中within函數(shù)，判定某個(gè)幾何體(注釋點(diǎn))是否存在另一個(gè)幾何體(多邊形)之內(nèi)。若為“True”，則保存在一個(gè)數(shù)組中，使注釋點(diǎn)與多邊形相互對應(yīng)，從而確定經(jīng)多邊形化(polygonize)后的區(qū)填充顏色值。主要代碼如下：

#迭代多邊形數(shù)組pgSet

for plg in pgSet:

pgInfoA = [ ]

#迭代注釋點(diǎn)數(shù)組

for i, pnt in enumerate(noteXY):

#判定注釋點(diǎn)是否在多邊形內(nèi)

isIn = Point(pnt).within(plg)

if isIn == True:

#保存在一個(gè)臨時(shí)數(shù)組中

pgInfoA.append (noteText[i])

#與多邊形對應(yīng)的注釋點(diǎn)數(shù)組

plgInfoSet.append (″.join(pgInfoA))

……

4) 生成區(qū)文件

在多邊形成區(qū)之前，預(yù)先編輯一個(gè)地層編碼和顏色對應(yīng)表(表1)，此表應(yīng)包含解釋圖中需著色所有地層、巖體或地質(zhì)體代碼(geoCode)；colorGIS為MapGIS專用顏色值，而colorHex為16進(jìn)制Hex顏色值，可供其他軟件調(diào)用。表1為內(nèi)蒙古克什克騰旗廣興地區(qū)音頻大地電磁測量項(xiàng)目的巖層(巖體)的顏色表。

表1 地層編碼和顏色對應(yīng)

確定區(qū)的填充顏色后，與表1進(jìn)行對比，依照MapGIS明碼文件格式[9]編輯區(qū)參數(shù)及區(qū)坐標(biāo)值，編寫wap文件。為保證識別率，在判斷時(shí)需注意：① 使用find函數(shù)，以避免同一個(gè)區(qū)可能會出現(xiàn)多個(gè)重復(fù)注釋點(diǎn)或不相關(guān)的注釋點(diǎn)；② 注意地層編碼順序，由簡單至復(fù)雜，比如P1d放置在P1d2和P1d1；③ NaN(空白區(qū))項(xiàng)，空白區(qū)或程序無法確定區(qū)，填充白色或者指定顏色。

3 實(shí)例

圖3顯示的是內(nèi)蒙古克什克騰旗廣興地區(qū)音頻大地電磁測量的地質(zhì)解釋素圖。為了使折線多邊形化，將曲線延長進(jìn)行“搭接”，并將上侏羅統(tǒng)滿克頭鄂博組和下二疊統(tǒng)大石寨組的不整合地質(zhì)界線設(shè)為圖中未出現(xiàn)的線顏色(如綠色)。在程序運(yùn)行之前，加載預(yù)先編輯的全區(qū)地層(巖體)地質(zhì)符號與顏色對應(yīng)表(表1)，作為區(qū)顏色填充依據(jù)。

圖3 手工繪制的地質(zhì)解釋斷面Fig.3 Section of geologic interpretation by manual drawing

圖4為經(jīng)程序處理后的繪制效果圖，圖面較為美觀，基本上符合最終圖件的要求。表示不整合地層界線的綠色折線變?yōu)轫樆牟ɡ司€，嚴(yán)格地隨波浪線位置填充顏色。除斷裂線外，地層(巖體)界線的外“搭接”部分均已刪除，此時(shí)，只要利用MapGIS程序的手工刪除斷裂線多余部分即可。

圖4 Python程序自動(dòng)化生成的地質(zhì)解釋斷面Fig.4 Section of geologic interpretation by automatic drawing (Python)

4 結(jié)論

Python語言在GIS方面同樣具有極為豐富的類庫。對于二維幾何體的分析與操作，Shapely庫功能強(qiáng)大且全面，再以Numpy為基礎(chǔ)使用Scipy解決相關(guān)的數(shù)學(xué)計(jì)算，以及Pandas的數(shù)據(jù)清洗功能，借助這些類庫可輕松地處理空間數(shù)據(jù)，并能縮短項(xiàng)目開發(fā)周期。相比常規(guī)方法，該繪圖程序主要有如下優(yōu)點(diǎn)：

1) 自動(dòng)生成波浪線(不整合界線)，可以直接進(jìn)行多邊形(區(qū))填充；

2) 通過線型和線形狀的判定，自動(dòng)刪除微短線，減少了繪圖工作環(huán)節(jié)；

3) 通過注釋點(diǎn)位置與多邊形的判定，識別區(qū)填充顏色。

總之，該程序脫離了商業(yè)軟件平臺，具有更好的擴(kuò)展性，依照項(xiàng)目制圖特點(diǎn)及規(guī)范要求，可高效、準(zhǔn)確地繪制出物探專題圖件，從而提高工作效率，使地質(zhì)物探技術(shù)人員將更多的精力投入到資料的分析與地質(zhì)解釋之中，從而得到更為豐富可靠的地質(zhì)成果。