gempy

GemPy - 3D 地質モデリング

クイックリファレンス

import gempy as gp

# Create model
geo_model = gp.create_geomodel(
    project_name='Model',
    extent=[0, 1000, 0, 1000, 0, 500],  # [xmin, xmax, ymin, ymax, zmin, zmax]
    resolution=[50, 50, 25]
)

# Add surface points and orientations
gp.add_surface_points(geo_model, x=[100, 500, 900], y=[500, 500, 500],
                      z=[400, 350, 400], surface='TopFormation')
gp.add_orientations(geo_model, x=[500], y=[500], z=[375],
                    pole_vector=[0, 0, 1], surface='TopFormation')

# Compute and visualize
gp.set_interpolator(geo_model)
sol = gp.compute_model(geo_model)
gp.plot_2d(geo_model, cell_number=[25], direction='y')

主なクラス

クラス	目的
GeoModel	メインコンテナ。すべてのモデルデータを保持します
StructuralFrame	地質学的な関係を管理します
Grid	補間に使用する計算メッシュです
Solutions	モデルの結果（岩相、スカラー場）です

基本的な操作

地質学的な関係を定義する

gp.map_stack_to_surfaces(geo_model, mapping={
    'Strata1': ['TopFormation', 'BaseFormation'],
    'Basement': ['Basement']
})
geo_model.structural_frame.structural_groups[0].structural_relation = \
    gp.data.StackRelationType.ERODE

断層を追加する

gp.add_surface_points(geo_model, x=[300, 300, 300], y=[200, 500, 800],
                      z=[100, 250, 400], surface='Fault1')
gp.add_orientations(geo_model, x=[300], y=[500], z=[250],
                    pole_vector=[1, 0, 0.5], surface='Fault1')
gp.map_stack_to_surfaces(geo_model, mapping={
    'Fault_Series': ['Fault1'],
    'Strata1': ['Layer1', 'Layer2'],
})
geo_model.structural_frame.structural_groups[0].structural_relation = \
    gp.data.StackRelationType.FAULT

ファイルからデータを読み込む

import pandas as pd

points_df = pd.read_csv('surface_points.csv')  # Columns: X, Y, Z, surface
gp.add_surface_points(geo_model, x=points_df['X'], y=points_df['Y'],
                      z=points_df['Z'], surface=points_df['surface'])

ori_df = pd.read_csv('orientations.csv')  # Columns: X, Y, Z, dip, azimuth, surface
gp.add_orientations(geo_model, x=ori_df['X'], y=ori_df['Y'], z=ori_df['Z'],
                    dip=ori_df['dip'], azimuth=ori_df['azimuth'],
                    surface=ori_df['surface'])

可視化

gp.plot_2d(geo_model, cell_number=[25], direction='y', show_data=True)
gp.plot_3d(geo_model, show_data=True, show_surfaces=True)  # Requires PyVista

結果にアクセスする

sol = gp.compute_model(geo_model)
lithology = sol.raw_arrays.lith_block
lith_3d = lithology.reshape(geo_model.grid.regular_grid.resolution)

構造関係の種類

種類	説明
ERODE	新しい地層が古い地層を侵食します（不整合）
ONLAP	新しい地層が古い地層の上にオンラップします
FAULT	地層面が断層面です
INTRUSION	貫入岩体です

代替手段との使い分け

シナリオ	推奨
地表データからの3D暗黙的地質モデリング	GemPy — 専用設計で、Python ネイティブです
構造フレームを使った複雑な褶曲モデリング	LoopStructural — 褶曲への対応が優れています
商用レベルの地下モデリング	SKUA-GOCAD — 業界標準のプロプライエタリなツールです
簡単な2D断面図のみが必要	GemPy でも可能ですが、過剰な場合があります。手動補間も検討してください
地質モデルから重力・磁気のフォワードモデルが必要	GemPy — ポテンシャル場に標準対応しています

GemPy を選ぶ場合：地表接触点と走向・傾斜から、断層や不整合の関係を含む暗黙的な3D地質モデリングが必要な場合です。特に Python のワークフローと統合する場合に適しています。標準的なケースでは、LoopStructural よりも習得しやすくなっています。

GemPy を避ける場合：地質構造が複雑な褶曲を中心とする場合（LoopStructural を使用してください）、または本番環境向けの高品質な貯留層モデリングが必要な場合（商用ツールを使用してください）です。

一般的なワークフロー

地表データから3D地質モデルを構築する

• 地表接触点（各地層につき最低2点）と走向・傾斜（各地層につき最低1つ）を読み込むか定義します
• 適切な範囲と解像度で GeoModel を作成します
• gp.add_surface_points() で地表点を追加します
• gp.add_orientations() で走向・傾斜を追加します
• gp.map_stack_to_surfaces() で地質スタックをマッピングします
• 構造関係（ERODE、ONLAP、FAULT）を設定します
• gp.set_interpolator() で補間器を設定します
• gp.compute_model() でモデルを計算します
• gp.plot_2d() で2D断面図を作成し、検証します
• gp.plot_3d() で3Dの結果を可視化するか、VTK にエクスポートします

よくある問題

問題	解決策
モデルを計算できない	各地層につき最低2点と1つの走向・傾斜があることを確認します
端部にアーティファクトが発生する	モデルの範囲をデータより広く設定します
断層のずれが正しくない	pole_vector の方向を確認します
メモリエラー	グリッドの解像度を下げます

参考資料

• 構造関係 — ERODE、ONLAP、FAULT の詳細ガイド
• データ要件 — 入力データの形式と制約
• トラブルシューティング — よくある問題と解決策

スクリプト

• scripts/validate_input.py — 入力データファイルを検証します
• scripts/export_model.py — モデルを VTK、numpy、または CSV にエクスポートします