구부러진 측선에 점군 투영

소개

지구과학 데이터를 다루면 측선을 따라 점군을 투영하는 작업이 발생할 수 있습니다. 예를 들면 지진 활동의 깊이를 나타내는 단면도 작성이나 반사법 지진파 탐사 데이터를 중합하기위한 데이터 처리 등이다. 측선이 동서남북 방향의 완전한 직선이면 실장은 매우 간단하지만, 현실의 측선은 그렇게 되지 않는 경우가 많다.
그러나 Python으로 도형 처리를 실시하는 Shapely 모듈( 공식 문서 )을 이용하면, 구부러진 측선이라도 매우 간단하게 실장할 수 있다(아래 그림).

샘플 데이터 생성

간단히하기 위해, 여기서는 무작위로 생성 된 2 차원 정규 분포에 따른 점군을 고려한다.

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 로 변환해 draw 할 수 있다.

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로 그려진 측선에 대해 바로 단면도를 만들 수 있는 코드로 했다.