在地质、物探工作中,自始至终都会应用不同形式的专题地图来直观地表现其研究成果^[1]。为获取最佳的地质效果,技术人员需不断地对地质情况深入研究,反复斟酌,经多次易稿后才能获取最终满意的成果,而伴随其中的专题图件也在不断的修改中完善。

在物探工作中,地质推断解释图是反映物探成果的必备专题图件之一。总体而言,地质推断解释图中的图元要素并非十分复杂;但由于修改次数较多,又受到专业制图软件功能的限制,工作效率难以满足快节奏项目需求。以目前国内应用范围最广的国产软件MapGIS为例,虽有强大编辑功能,但凭手工绘图方式,其效率依然难有明显提升。如电法勘探中的反演电阻率断面的地质推断解释图,不整合线是一种常见的地质现象,图中使用波浪线表示,而在MapGIS软件中波浪线则由一个特殊线型呈现,实际坐标位置处于曲线的中轴线上;若直接成区,线与区间的分界面无法满足制图要求,仍需徒手绘制或借助于第三方软件生成,操作繁琐且费时费力。笔者针对该类项目制图特点和需求,利用Python语言及其强大的第三方库,设计了自动制图程序,主要功能包括波浪线(不整合界线)、自动剪断线、微短线识别与删除以及多边形化成区等功能,目的是提高工作效率,减少制图过程中不必要的错误发生。

Python语言是一个纯净、简洁明了的编程语言,特别适合于初学者,并且还可帮助专业人员快速地解决复杂的问题,它可以轻松快速地对地理空间数据进行可视化处理^[2]。Python具有强大的科学计算、绘图功能,同时在地图绘制、地理空间数据的处理与转换方面具有丰富的类库^[3]。本文主要涉及了NumPy(N维数组和矩阵计算)、Shapely(二维图形处理库)、Pandas(数据清洗库)和Scipy(数学工具包和算法库)。另外,在程序调试中,还使用了Geopandas(地图绘制)和GeoJson(空间数据格式)库,用于预览绘制效果。

Shapely为Python中二维图形软件包,基于广为流传的GEOS库中的函数,进行几何体的操作和分析^[4]。本文主要涉及该库中的几何体裁剪、多边形化、节点抽稀和几何体空间位置判定等函数^[5]。

Scipy是构建在Python扩展库Numpy上的数学算法和便捷函数的集合。增加了强大的交互式Python任务,向用户提供了高级别的指令或类用于数据操作和可视化^[6,7]。本文主要利用该库提供的插值(interpolate)模块。

一般情况下,当物探工程师在绘制好地质解释图后,交付给制图人员进行线编辑、颜色填充、图面整饰等工作。利用MapGIS的制图流程大体包括自动裁剪线、线转弧段、拓扑重建、修改地层(岩体)颜色、手工去除微短线等。笔者基于Python语言及其Shapely、NumPy、Scipy和Pandas等扩展库,参照上述MapGIS绘制流程,同时也补充了不整合线(波浪线)的绘制,达到了自动制图工作的目的。程序实现的基本流程见图1。

MapGIS图形按点、线、区分类存储,与之对应的文件分别为点文件、线文件、区文件,明码格式文件名后缀分别为.wat、.wal和.wap^[8],详细格式可查阅MapGIS相关文献^[9]。为脱离MapGIS运行环境,扩展程序的适用性,此次调用了MapGIS的明码格式。

2.2.1 不整合界线的绘制

在上、下地层之间存有一个沉积间断面,叫不整合面,在地面的出露线称作不整合线^[10]。在断面图中,地层不整合界线一般用波浪线表示。由于MapGIS程序中提供波浪线线型无法满足区域填充的要求,首先通过程序生成一条适合区域填充的波浪线,从图形学上,可以一条曲线或折线为轴心,由多个连续的正弦曲线或余弦曲线绘制成波浪线。绘制过程主要包括曲线插值平滑曲线,等路径插值提取节点坐标,根据节点位置绘制正弦曲线或余弦曲线形成波浪线,以及节点抽稀减轻运算负荷等4部分(图2)。

1) 曲线插值

曲线插值的目的是光滑曲线,使图面更为美观,同时又要避免曲线过度变形。在Scipy库中,scipy.interpolate.splprep是适合于多维曲线的B样条插值法,经多次试验,s(光滑系数)取0,k(曲线拟合度)取2较为合适。Python实现代码如下:

from scipy.interpolate import splprep, splev

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

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

#获取XY坐标值

xySet = splev(u, tck)

2) 等路径提取节点坐标

在Matlab中提供了沿路径插值法(interpolate along the path),而Scipy并未提供类似的函数。在Shapely库中interpolate函数功能是沿线性几何体获得指定距离的坐标值,通过等距离迭代计算可获得类似效果。

#获得线的长度

myLine = LineString(xy0)

lineLen = myLine.length

# sizeUnit为单个正弦曲线的波长

segNum = int(lineLen/sizeUnit)+2

dataSet3=[]

for i in range(segNum):

#获取等路径的坐标值,存入数组中

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

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

3) 绘制波浪线

为了防止正弦(余弦)曲线变形,将预先生成的正弦曲线或余弦曲线,通过旋转、平移操作移至曲线的相应位置,旋转平移的代码如下:

#空间几何体的旋转、平移函数

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,交替绘制。此时,需要考虑节点的矢量方向,确保正弦(余弦)曲线衔接连续。代码如下:

if i % 2 == 0: #偶数点

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

#当后一点坐标值X大于前一点时,

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

else:

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

#翻转数组

xx1 = np.flipud(xx1)

yy1 = np.flipud(yy1)

else:#奇数点

(余弦曲线与正弦曲线平移类似,略)

#存储节点位置

allXPoints.append(xx1a)

allYPoints.append(yy1a)

4) 节点抽稀

此时波浪线的节点十分致密,将会加大计算机运行负荷,尤其是长剖面,线节点过多会导致MapGIS无法正常运行。线节点的抽稀方法采用道格拉斯—普克算法(douglas-peucker algorithm),在Shapely中提供了此功能(simplify函数)。经试验对比,抽稀系数为0.01较为适合。代码如下:

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

2.2.2 自动剪断线与微短线识别

从数学上,一个多边形(ploygon)由同一平面上闭合曲线构成,填充后可形成区。因此,绘制图件时,为形成一个闭合多边形,构成多边形轮廓的曲线必须有“搭接”。此环节可使用shapely.opt库中的split函数,通过迭代计算,进行不同几何体间相互分割,实现要求:① 被分割几何体为单个线几何体(linestring);② 分割几何体为被分割体外的其他所有线几何体的集合(multilinestring),确保分割几何体必须一次性完成剪断任务。关键代码如下:

#导入shapely库

from shapely.ops import polygonize, split

#遍历所有曲线数组

for i in range (len(lineSets)):

#生成第i条线(被裁剪)

oneLine = LineString(lineSets [i])

#从线数组中删除第i条线

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

#生成除第i条线外的线集合

otherLines = MultiLineString(others)

#裁剪线

resultCollect = split(oneLine, otherLines)

#获取线集合中线的数量

lineNum = len(resultCollect.geoms)

获得的线裁剪结果(result collect)为多线集合,此集合内几何体排列顺序与原始线坐标位置一致,从而为判定或删除(地层)线两端外延部分奠定了条件。结合断面图特点,主要判定方法包括:① 判定特殊线型,如红色的断裂线和闭合曲线(地质体、图框),此类线几何体将不进行任何删除处理;断裂线根据线属性(颜色)判定,闭合曲线可根据首尾点的坐标值确定。② 根据多线集合中几何体的数量(lineNum≥3)判定线几何体否已被裁剪,来决定是否删除第一个和最后一个几何体(线)。最后以MapGIS明码格式保存线参数(节点坐标值、颜色、宽度、线型等)。

2.2.3 折线成区和颜色填充

此过程相当于MapGIS中的线转弧段和拓扑成区。但MapGIS是随机填充颜色,而本程序是根据地层、岩体的注释点进行颜色填充。

1) 折线转多边形

此处的折线转多边形相当于MapGIS中的“线转弧段”,Shapely提供了polygonize函数,代码如下:

#折线的多边形化,outLinSets为折线数组

polygonSet = polygonize(outLinSets)

# 提取多边形的坐标值

pgSet = list(polygonSet)

2) 读取点文件

使用Pandas读取wat文件更为方便,主要获取注释点(X,Y)坐标值、字体颜色和标注。当然,也可以限制一些与颜色填充无关的注释点(非强制),如断裂编号,为后续地层符号识别去除一些干扰因素。

#使用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']

#去除红色断裂注释点(颜色为6)

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

3) 确定多边形内的注释点

通过对多边形和注释点双迭代运算,利用Shapely中within函数,判定某个几何体(注释点)是否存在另一个几何体(多边形)之内。若为“True”,则保存在一个数组中,使注释点与多边形相互对应,从而确定经多边形化(polygonize)后的区填充颜色值。主要代码如下:

#迭代多边形数组pgSet

for plg in pgSet:

pgInfoA = [ ]

#迭代注释点数组

for i, pnt in enumerate(noteXY):

#判定注释点是否在多边形内

isIn = Point(pnt).within(plg)

if isIn == True:

#保存在一个临时数组中

pgInfoA.append (noteText[i])

#与多边形对应的注释点数组

plgInfoSet.append (″.join(pgInfoA))

……

4) 生成区文件

在多边形成区之前,预先编辑一个地层编码和颜色对应表(表1),此表应包含解释图中需着色所有地层、岩体或地质体代码(geoCode);colorGIS为MapGIS专用颜色值,而colorHex为16进制Hex颜色值,可供其他软件调用。表1为内蒙古克什克腾旗广兴地区音频大地电磁测量项目的岩层(岩体)的颜色表。

确定区的填充颜色后,与表1进行对比,依照MapGIS明码文件格式^[9]编辑区参数及区坐标值,编写wap文件。为保证识别率,在判断时需注意:① 使用find函数,以避免同一个区可能会出现多个重复注释点或不相关的注释点;② 注意地层编码顺序,由简单至复杂,比如P₁d放置在P₁d²和P₁d¹;③ NaN(空白区)项,空白区或程序无法确定区,填充白色或者指定颜色。

图3显示的是内蒙古克什克腾旗广兴地区音频大地电磁测量的地质解释素图。为了使折线多边形化,将曲线延长进行“搭接”,并将上侏罗统满克头鄂博组和下二叠统大石寨组的不整合地质界线设为图中未出现的线颜色(如绿色)。在程序运行之前,加载预先编辑的全区地层(岩体)地质符号与颜色对应表(表1),作为区颜色填充依据。

图4为经程序处理后的绘制效果图,图面较为美观,基本上符合最终图件的要求。表示不整合地层界线的绿色折线变为顺滑的波浪线,严格地随波浪线位置填充颜色。除断裂线外,地层(岩体)界线的外“搭接”部分均已删除,此时,只要利用MapGIS程序的手工删除断裂线多余部分即可。

Python语言在GIS方面同样具有极为丰富的类库。对于二维几何体的分析与操作,Shapely库功能强大且全面,再以Numpy为基础使用Scipy解决相关的数学计算,以及Pandas的数据清洗功能,借助这些类库可轻松地处理空间数据,并能缩短项目开发周期。相比常规方法,该绘图程序主要有如下优点:

1) 自动生成波浪线(不整合界线),可以直接进行多边形(区)填充;

2) 通过线型和线形状的判定,自动删除微短线,减少了绘图工作环节;

3) 通过注释点位置与多边形的判定,识别区填充颜色。

总之,该程序脱离了商业软件平台,具有更好的扩展性,依照项目制图特点及规范要求,可高效、准确地绘制出物探专题图件,从而提高工作效率,使地质物探技术人员将更多的精力投入到资料的分析与地质解释之中,从而得到更为丰富可靠的地质成果。