はじめに

地球科学データを扱っていると、測線沿いに点群を投影するような作業が生じることがある。例えば地震活動の深さを示す断面図作成であったり、反射法地震波探査データを重合するためのデータ処理などである。測線が東西南北方向の完全な直線であれば実装は極めて簡単だが、現実の測線はそうはいかないことが多い。
しかしPythonで図形処理を行うShapelyモジュール(公式ドキュメント)を利用すると、折れ曲がった測線であっても極めてシンプルに実装できる(下図)。

サンプルデータの生成

簡単のため、ここではランダムに生成した二次元正規分布に従った点群を考える。

XYZ_PolygonalInterpolation.py

import numpy as np
import matplotlib.pyplot as plt
## 0を中心とした正規分布に従う100点の点群を生成
X,Y = np.random.randn(100), np.random.randn(100) 

## 描画
fig, ax = plt.subplots()
ax.scatter(X,Y, s=10)
ax.axis('equal') #X軸とY軸を等しくする
plt.show()

測線の設定

ここに以下のような折れ曲がった測線を仮定し、Shapelyのbuffer()によって測線に一定の幅をもたせたポリゴンを生成する。(この辺の仕組みはSVGと同じである)
ポリゴンの内側にある点はcontains()により検出できる。
Shapelyオブジェクトはdescartesモジュール(デカルト? 同名のプロジェクトが沢山あるので注意)を用いることで簡単にmatplotlibのpatchへ変換し描画できる。

XYZ_PolygonalInterpolation.py

## ShapelyおよびShapelyオブジェクトをMatplotlibへ変換するdescartesモジュールを使う
from shapely.geometry import Point, Polygon, LineString
from descartes import PolygonPatch

## 折れ線の座標, 折れ線のShapelyオブジェクト
path = np.array([[-2,-1],[-0.5,-0.5],[0.5,0.5],[2,1]]) 
line = LineString(path)
## 折れ線に幅(0.5)をもたせて生成したポリゴンのShapelyオブジェクト。capとjoinで終端や変曲点の取り扱いを指定する
poly = line.buffer(0.5, cap_style=2, join_style=2)

## 点群の中でポリゴンに含まれるもののインデックス
idx = [i for i in range(len(X)) \
       if poly.contains(Point((X[i],Y[i])))]

## Matplotlibによる描画
fig, ax = plt.subplots()
ax.scatter(X,Y, c='black', s=2) #全ての点を黒で
ax.scatter(np.take(X,idx), np.take(Y,idx), c='red', s=5) #ポリゴン内の点を赤で
ax.add_patch(PolygonPatch(poly, facecolor=(0,0,0,0), edgecolor=(1,0,0,1))) #ポリゴンをdescartesモジュールを使って描画。色はRGBAで指定している。
ax.plot(path[:,0], path[:,1], color=(1,0,0,1)) #測線
ax.arrow(x=path[-1,0], y=path[-1,1], \
         dx=(path[-1,0]-path[-2,0])*0.2, \
         dy=(path[-1,1]-path[-2,1])*0.2, \
         width=0.1, color=(1,0,0,1)) #おまけの矢印。
ax.axis('equal'); plt.show()

図形への投影

Shapelyのproject()を用いることで、ポリゴン内の点を直線に投影し、始点からの距離を求めることができる。これを横軸に用いると折れ線沿いのX軸とすることができる。

XYZ_PolygonalInterpolation.py

## ポリゴン内の点を直線へ投影し、直線上の始点からの距離を得る
lengths = np.array([line.project(Point((X[i],X[i]))) for i in idx])
lmin, lmax = np.floor(lengths.min()), np.ceil(lengths.max()) #グラフ1および2描画用の最大最小値

## Matplotlibによる描画
fig, ax = plt.subplots(3) #3つのグラフを同時に描く
## グラフ0は前図の省略版
ax[0].scatter(X,Y, c='black', s=2)
ax[0].add_patch(PolygonPatch(poly, facecolor=(0,0,0,0), edgecolor=(1,0,0,1)))
ax[0].plot(path[:,0], path[:,1], color=(1,0,0,1))
ax[0].axis('equal')
## グラフ1は各点の直線からの距離
ax[1].scatter(x=lengths, y=distances) #代わりにyに震源深さ(Z)などを与えれば断面図になる
ax[1].set_ylabel('Distances from line')
ax[1].set_xlim(lmin, lmax) #グラフ2と合わせる
## グラフ2は投影された点群数のヒストグラム
ax[2].hist(lengths, bins=np.arange(lmin,lmax,0.5))
ax[2].set_xlim(lmin, lmax) #グラフ1と合わせる
ax[2].set_ylabel('Counts along line')
ax[2].set_xlabel('Length along line')
plt.show()

※画像の点線はわかりやすくするためにInkscapeで付加した

おわりに

色々な応用方法が考えられる。例えば自分の場合、shapefileを読めるようにしてQGISで描いたshapefileで描かれた測線に対してすぐ断面図を作ることができるようなコードにした。