geostatspy

GeostatsPy - 地球統計解析

クイックリファレンス

import geostatspy.GSLIB as GSLIB
import geostatspy.geostats as geostats
import pandas as pd

df = pd.read_csv('data.csv')
df['npor'], tvpor, tnspor = geostats.nscore(df, 'porosity')  # Transform

lag, gamma, npairs = geostats.gamv(df, 'X', 'Y', 'npor',      # Variogram
    tmin=-9999, tmax=9999, xlag=50, xltol=25, nlag=15,
    azm=0, atol=22.5, bandwh=9999, bandwd=9999)

vario = GSLIB.make_variogram(nug=0.0, nst=1, it1=1, cc1=1.0,  # Model
                              azi1=0, hmaj1=300, hmin1=300)

est, var = geostats.kb2d(df, 'X', 'Y', 'npor', ..., vario=vario)  # Krige

主な関数

カテゴリ	関数
可視化	locmap, pixelplt, hist
バリオグラム	gamv, vmodel
クリギング	kb2d, kb3d
シミュレーション	sgsim, sisim
変換	nscore, backtr
デクラスタリング	declus

一般的な操作

1. 正規スコア変換

df['npor'], tvpor, tnspor = geostats.nscore(df, 'porosity')
original = geostats.backtr(nscore_data, tvpor, tnspor, zmin=0, zmax=0.3)

2. 実験バリオグラム

lag, gamma, npairs = geostats.gamv(
    df, 'X', 'Y', 'npor',
    tmin=-9999, tmax=9999,    # Trimming limits
    xlag=50, xltol=25,        # Lag distance, tolerance
    nlag=15, azm=0, atol=22.5, bandwh=9999, bandwd=9999)

3. バリオグラムモデル

# Types: 1=spherical, 2=exponential, 3=gaussian
vario = GSLIB.make_variogram(
    nug=0.0, nst=1,            # Nugget, number of structures
    it1=1, cc1=1.0,            # Type, sill contribution
    azi1=0, hmaj1=300, hmin1=300)  # Azimuth, major/minor range

4. クリギング（kb2d）

est, var = geostats.kb2d(
    df, 'X', 'Y', 'npor', tmin=-9999, tmax=9999,
    nx=50, xmn=25, xsiz=50,    # Grid X: ncells, origin, size
    ny=50, ymn=25, ysiz=50,    # Grid Y
    nxdis=1, nydis=1, ndmin=1, ndmax=10,
    radius=500, ktype=0, skmean=0.0, vario=vario)  # ktype: 0=simple, 1=ordinary

5. 逐次ガウスシミュレーション

sim = geostats.sgsim(
    df, 'X', 'Y', 'npor', wcol=-1, scol=-1,
    tmin=-9999, tmax=9999, itrans=0,
    ismooth=0, dession=0, dmession=0,
    zmin=-4, zmax=4, ltail=1, ltpar=0, utail=1, utpar=0,
    nsim=1, nx=50, xmn=25, xsiz=50, ny=50, ymn=25, ysiz=50,
    nz=1, zmn=0, zsiz=1, seed=73073,
    ndmin=1, ndmax=10, nodmax=10, radius=500, radius1=500,
    sang1=0, sang2=0, sang3=0, mxctx=10, mxcty=10, mxctz=1,
    ktype=0, vario=vario)

6. デクラスタリング

wts, cell_size, ncut = geostats.declus(
    df, 'X', 'Y', 'porosity', iminmax=1, noff=10, ncell=20, cmin=10, cmax=500)
declustered_mean = np.average(df['porosity'], weights=wts)

バリオグラムモデル

コード	モデル	用途
1	球状	最も一般的で、有限のレンジを持ちます
2	指数	漸近的にシルに達します
3	ガウス	非常に滑らかで、原点付近では放物線状です
4	パワー	非有界で、フラクタル状です

主なパラメーター

パラメーター	説明
nug	ナゲット効果（測定誤差＋マイクロスケールの変動）
sill	全分散（ナゲット＋構造による寄与）
range	相関が無視できる程度になる距離
azimuth	最大連続性の方向（北からの角度）
ktype	0＝単純クリギング（平均が既知）、1＝通常クリギング

代替手段との使い分け

用途	ツール	理由
GSLIB 形式のワークフロー	GeostatsPy	GSLIB プログラムを Python に直接移植しています
SGSIM / SISIM シミュレーション	GeostatsPy	GSLIB シミュレーションエンジンを完全に備えています
空間データのデクラスタリング	GeostatsPy	declus 関数を標準搭載しています
最新のバリオグラム API	scikit-gstat	API がすっきりしており、sklearn と統合できます
クリギングのみ（シミュレーションなし）	pykrige	用途を絞った API で、汎用クリギングをサポートしています
ランダム場の生成	gstools	柔軟な共分散モデルと場の生成機能を備えています
大規模な 3D 地質モデリング	SGeMS / Petrel	GUI ベースで、産業用途のワークフローに適しています
指示シミュレーション	GeostatsPy（sisim）	カテゴリカルな属性のシミュレーション

GeostatsPy を選ぶ場合：Python で GSLIB 互換のワークフローが必要な場合、特に逐次シミュレーション（SGSIM/SISIM）やデクラスタリングを行う場合、または GSLIB のパラメーター規約に慣れている場合に適しています。貯留層特性評価のワークフローに最適です。

scikit-gstat を選ぶ場合：バリオグラム解析とクリギングに、最新の scikit-learn 形式の API を使いたい場合に適しています。Python のデータサイエンスワークフローとの統合性にも優れています。

pykrige を選ぶ場合：クリギング補間のみ（シミュレーションなし）が必要で、外部ドリフトまたは回帰クリギング機能を備えた汎用クリギングを利用したい場合に適しています。

一般的なワークフロー

バリオグラム解析とクリギング補間

• 空間データを pandas の DataFrame に読み込む
• GSLIB.locmap() と GSLIB.hist() でデータを調べる
• クラスタリングを確認し、必要に応じて geostats.declus() でデクラスタリングする
• geostats.nscore() で正規スコア変換を適用する
• geostats.gamv() で実験バリオグラムを計算する（まずは等方性）
• 異方性を確認するため、方向別バリオグラムを調べる（方位角 0、45、90、135）
• GSLIB.make_variogram() でバリオグラムモデルを当てはめる
• 実験バリオグラムにモデルを重ね、適合度を確認する
• geostats.kb2d() でクリギングを実行する（通常クリギングでは ktype=1）
• 不確実性を定量化するため、SGSIM を実行する（50～100 回の実現）
• geostats.backtr() で結果を逆変換する
• 交差検証またはホールドアウトデータで検証する

よくある問題

問題	解決策
バリオグラムがシルに達しない	全レンジを捉えられるように nlag または xlag を増やします
クリギングで負の値が生成される	クリギングの前ではなく、後に逆変換します
SGSIM のアーティファクト	グリッドの定義（xmn、xsiz）がデータの範囲と一致しているか確認します
バリオグラムのペア数が少なすぎる	atol（角度許容差）または xltol（ラグ許容差）を増やします
バリオグラムのホール効果	周期性を示している可能性があります。ネスト構造を試してください

ヒント

1. 必ず正規スコアに変換します - ほとんどの手法はガウス分布を前提としています
2. 等方性から始めます - 方向別バリオグラムで妥当性が確認できた場合のみ、異方性を追加します
3. バリオグラムのペアを確認します - 各ラグで十分なペア数があることを確認します
4. 複数の実現を用います - 不確実性の定量化には50～100以上を使用します
5. 最後に逆変換します - 最終結果にのみ適用します

参考資料

• GSLIB Programs - GSLIBプログラムのリファレンス
• Simulation Methods - SGSIM、SISIM、その他の手法

スクリプト

• scripts/kriging_example.py - クリギングのワークフロー全体の例