harmonica

Harmonica - 重力・磁気

クイックリファレンス

import harmonica as hm
import numpy as np

# Forward model - prism gravity
prism = [-500, 500, -500, 500, -2000, -500]  # (west, east, south, north, bottom, top)
gravity = hm.prism_gravity(coordinates, prism, density=500, field='g_z')

# Terrain correction
layer = hm.prism_layer((easting, northing), surface=topo, reference=0,
                        properties={'density': 2670})
terrain_effect = layer.gravity(coordinates, field='g_z')

# Equivalent source gridding
eqs = hm.EquivalentSources(depth=10000, damping=10)
eqs.fit(coordinates, gravity_data)
grid = eqs.grid(spacing=5000, data_names=['gravity'])

# Upward continuation (requires gridded xarray)
upward = hm.upward_continuation(gravity_grid, height_displacement=1000)

主な関数

関数	目的
point_gravity	点質量による重力
prism_gravity	直方体による重力
tesseroid_gravity	球面プリズムによる重力（地域・全球）
prism_magnetic	プリズムによる磁気異常
prism_layer	地形からプリズム層を作成
EquivalentSources	等価ソースを使って散在データをグリッド化
upward_continuation	FFT ベースの上方接続
bouguer_correction	単純ブーゲー板補正

基本操作

順計算モデル - 直方体プリズム

# Define prism: (west, east, south, north, bottom, top) in meters
prism = [-500, 500, -500, 500, -2000, -500]
density = 500  # kg/m3 density contrast

# Observation grid
x_obs, y_obs = np.meshgrid(np.linspace(-5000, 5000, 100), np.linspace(-5000, 5000, 100))
z_obs = np.zeros_like(x_obs)

# Calculate gravity (mGal). Fields: 'g_z', 'g_north', 'g_east', 'potential'
gravity = hm.prism_gravity((x_obs.ravel(), y_obs.ravel(), z_obs.ravel()),
                           prism, density, field='g_z')

地形補正

import xarray as xr

topo = xr.open_dataarray('dem.nc')
layer = hm.prism_layer((topo.easting.values, topo.northing.values),
                       surface=topo.values, reference=0,
                       properties={'density': 2670})
terrain_effect = layer.gravity((obs_easting, obs_northing, obs_height), field='g_z')
bouguer_anomaly = free_air_anomaly - terrain_effect

等価ソースによるグリッド化

import verde as vd

# Project to Cartesian
projection = vd.get_projection(longitude, latitude)
easting, northing = projection(longitude, latitude)

eqs = hm.EquivalentSources(depth=10000, damping=10)
eqs.fit((easting, northing, altitude), gravity_mgal)
grid = eqs.grid(spacing=5000, data_names=['gravity'])

磁気順計算モデル

prism = [-500, 500, -500, 500, -2000, -500]
magnetization = hm.magnetic_vector(intensity=5.0, inclination=60, declination=10)
b_total = hm.prism_magnetic(coordinates, prism, magnetization, field='b_total')

微分フィルター

dx = hm.derivative_easting(gravity_grid)
dy = hm.derivative_northing(gravity_grid)
dz = hm.derivative_upward(gravity_grid)

thg = np.sqrt(dx**2 + dy**2)  # Total horizontal gradient
tilt = np.arctan2(dz, thg)     # Tilt angle

座標系

Harmonica は右手座標系を使用します。

• 東距（x）：東向きが正
• 北距（y）：北向きが正
• 上向き（z）：上向きが正（標高は正、深度は負）

単位は SI 単位系です。距離はメートル、密度は kg/m3、重力は mGal です。

代替手段との使い分け

ユースケース	ツール	理由
重力・磁気の順計算モデリング	Harmonica	目的に特化し、Fatiando エコシステムに属しています
ポテンシャル場の逆解析	SimPEG	正則化を備えた包括的な逆解析フレームワークです
商用重力データ処理	Oasis Montaj	業界標準の GUI と独自形式に対応しています
単純なブーゲー補正のみ	Custom numpy	単発の計算に必要な依存関係が少なくて済みます
等価ソースによるグリッド化	Harmonica	ポテンシャル場に最適なオープンソースの選択肢です
地域・全球スケール	Harmonica（テッセロイド）	球面幾何をネイティブに扱えます
磁気データの極化	Harmonica	グリッドデータ向けの FFT ベースのフィルターを備えています
教育・プロトタイピング	Harmonica	分かりやすい API と充実したドキュメントを備えています

Harmonica を選ぶ場合：順計算モデリング、地形補正、または等価ソースによるグリッド化を含む、オープンソースの重力・磁気データ処理が必要な場合です。座標変換とグリッド化には Verde と適切に連携します。Fatiando a Terra エコシステムの一部です。

SimPEG を選ぶ場合：地下の物性分布（密度または帯磁率モデル）を求めるために、ポテンシャル場データの逆解析が必要な場合です。

Oasis Montaj を選ぶ場合：独自形式、商用サポート、または GUI ベースの対話的な処理が必要な業務環境で作業する場合です。

一般的なワークフロー

地形補正とグリッド化を伴う重力調査の処理

• 生の重力観測データと観測点座標を読み込む
• 緯度補正、フリーエア補正、潮汐補正を適用する
• DEM を読み込み、hm.prism_layer() でプリズム層を構築する
• 観測点で地形効果を計算する
• フリーエア異常から地形効果を差し引き、ブーゲー異常を求める
• verde.get_projection() で座標を直交座標に投影する
• 観測点密度に偏りがある場合は、データをブロック化して縮約する
• hm.EquivalentSources() で等価ソースを適合させる
• ブーゲー異常を規則格子にグリッド化する
• vd.distance_mask() を適用し、データから離れた領域をマスクする
• 解釈用に微分フィルター（水平勾配、傾斜角）を適用する
• 上方接続を行い、広域的な特徴を強調する
• グリッドデータを NetCDF にエクスポートする

よくある問題

問題	解決策
重力の符号が誤っている	z 軸の規約（上向きが正）を確認する
等価ソースの適合が不良	depth と damping のパラメーターを調整する
地形補正が遅い	DEM の解像度を下げるか、より大きなプリズムを使用する
FFT フィルターで端点効果が生じる	upward_continuation の適用前にグリッドをパディングする
座標が一致しない	投影座標と地理座標のどちらを使うか一貫させる

ヒント

1. ローカルな調査では投影座標（メートル単位）を使用してください
2. 地域・全球スケールのモデリングではテッセロイドを使用してください
3. 等価ソースは不規則なデータ間隔を適切に処理します
4. 適切な密度を選択してください（上部地殻では2670 kg/m3が一般的です）
5. 符号規約を確認してください。深度は負のz値です

参考資料

• 重力・磁気補正 - 標準的な補正と異常
• 順方向モデリング - 順方向モデリングの詳細な手法

スクリプト

• scripts/gravity_processing.py - 重力調査データを処理します