Back to Details

harmonica

Compare original and translation side by side

🇺🇸

Original

English
🇨🇳

Translation

Chinese

Harmonica - Gravity and Magnetics

Harmonica - 重力与磁数据处理

Quick Reference

快速参考

python
import harmonica as hm
import numpy as np
python
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')
prism = [-500, 500, -500, 500, -2000, -500] # (西, 东, 南, 北, 底部深度, 顶部深度) 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')
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'])
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)

向上延拓（需要网格化的xarray数据）

upward = hm.upward_continuation(gravity_grid, height_displacement=1000)
undefined
upward = hm.upward_continuation(gravity_grid, height_displacement=1000)
undefined

Key Functions

核心函数

FunctionPurpose
point_gravity
Gravity from point masses
prism_gravity
Gravity from rectangular prisms
tesseroid_gravity
Gravity from spherical prisms (regional/global)
prism_magnetic
Magnetic anomaly from prisms
prism_layer
Create layer of prisms from topography
EquivalentSources
Grid scattered data with equivalent sources
upward_continuation
FFT-based upward continuation
bouguer_correction
Simple Bouguer plate correction
函数用途
point_gravity
点质量的重力计算
prism_gravity
矩形棱柱体的重力计算
tesseroid_gravity
球形棱柱体的重力计算（区域/全球尺度）
prism_magnetic
棱柱体的磁异常计算
prism_layer
基于地形创建棱柱体层
EquivalentSources
利用等效源对离散数据进行网格化
upward_continuation
基于FFT的向上延拓
bouguer_correction
简易布格平板校正

Essential Operations

关键操作

Forward Model - Rectangular Prism

正演模型 - 矩形棱柱体

python
undefined
python
undefined

Define prism: (west, east, south, north, bottom, top) in meters

定义棱柱体：(西, 东, 南, 北, 底部深度, 顶部深度)，单位为米

prism = [-500, 500, -500, 500, -2000, -500] density = 500 # kg/m3 density contrast
prism = [-500, 500, -500, 500, -2000, -500] density = 500 # 密度差，单位为千克/立方米

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)
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'

计算重力（单位为毫伽）。可选场分量：'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')
undefined
gravity = hm.prism_gravity((x_obs.ravel(), y_obs.ravel(), z_obs.ravel()), prism, density, field='g_z')
undefined

Terrain Correction

地形校正

python
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
python
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

Equivalent Source Gridding

等效源网格化

python
import verde as vd
python
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'])
undefined
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'])
undefined

Magnetic Forward Model

磁学正演模型

python
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')
python
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')

Derivative Filters

导数滤波器

python
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
python
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)  # 总水平梯度
tilt = np.arctan2(dz, thg)     # 倾斜角

Coordinate System

坐标系

Harmonica uses a right-handed coordinate system:
  • Easting (x): positive east
  • Northing (y): positive north
  • Upward (z): positive up (heights positive, depths negative)
Units are SI: meters for distance, kg/m3 for density, mGal for gravity.
Harmonica采用右手坐标系：
  • 东向(x): 东为正方向
  • 北向(y): 北为正方向
  • 向上(z): 向上为正方向（高度为正，深度为负）
单位采用国际单位制：距离单位为米，密度单位为千克/立方米，重力单位为毫伽（mGal）。

When to Use vs Alternatives

适用场景与替代工具对比

Use CaseToolWhy
Gravity/magnetic forward modellingHarmonicaPurpose-built, Fatiando ecosystem
Potential field inversionSimPEGFull inversion framework with regularization
Commercial gravity processingOasis MontajIndustry-standard GUI, proprietary formats
Simple Bouguer corrections onlyCustom numpyFewer dependencies for one-off calculations
Equivalent source griddingHarmonicaBest open-source option for potential fields
Regional/global scaleHarmonica (tesseroids)Handles spherical geometry natively
Magnetic data reduction to poleHarmonicaFFT-based filters for gridded data
Teaching/prototypingHarmonicaClean API, good documentation
Choose Harmonica when: You need open-source gravity/magnetic processing with forward modelling, terrain corrections, or equivalent source gridding. It integrates well with Verde for projections and gridding. Part of the Fatiando a Terra ecosystem.
Choose SimPEG when: You need to invert potential field data for subsurface property distributions (density or susceptibility models).
Choose Oasis Montaj when: You work in an industry setting that requires proprietary formats, commercial support, or GUI-based interactive processing.
使用场景工具原因
重力/磁学正演模拟Harmonica专属工具，属于Fatiando生态系统
位场反演SimPEG带正则化的完整反演框架
商用重力数据处理Oasis Montaj行业标准GUI，支持专有格式
仅需简易布格校正自定义numpy代码依赖更少，适合一次性计算
等效源网格化Harmonica位场处理领域最佳开源方案
区域/全球尺度建模Harmonica（tesseroids）原生支持球面几何
磁数据化极处理Harmonica基于FFT的网格化数据滤波器
教学/原型开发HarmonicaAPI简洁，文档完善
选择Harmonica的场景：需要开源的重力/磁数据处理工具，支持正演模拟、地形校正或等效源网格化。它与Verde的投影和网格化工具集成良好，属于Fatiando a Terra生态系统。
选择SimPEG的场景：需要对位场数据进行反演，以得到地下属性分布（密度或磁化率模型）。
选择Oasis Montaj的场景：工作环境要求使用专有格式、商用支持或基于GUI的交互式处理。

Common Workflows

常见工作流程

Process Gravity Survey with Terrain Correction and Gridding

带地形校正与网格化的重力测量数据处理

  • Load raw gravity observations and station coordinates
  • Apply latitude, free-air, and tidal corrections
  • Load DEM and build prism layer with 
    hm.prism_layer()
  • Compute terrain effect at observation points
  • Subtract terrain effect from free-air anomaly to get Bouguer anomaly
  • Project coordinates to Cartesian with 
    verde.get_projection()
  • Block-reduce data if station density is uneven
  • Fit equivalent sources with 
    hm.EquivalentSources()
  • Grid the Bouguer anomaly onto a regular grid
  • Apply 
    vd.distance_mask()
    to mask areas far from data
  • Apply derivative filters (horizontal gradient, tilt angle) for interpretation
  • Perform upward continuation to enhance regional features
  • Export gridded data to NetCDF
  • 加载原始重力观测数据和测站坐标
  • 应用纬度校正、自由空气校正和潮汐校正
  • 加载DEM并使用
    hm.prism_layer()
    构建棱柱体层
  • 计算观测点处的地形影响
  • 从自由空气异常中减去地形影响得到布格异常
  • 使用
    verde.get_projection()
    将坐标投影到笛卡尔坐标系
  • 若测站密度不均，对数据进行分块降采样
  • 使用
    hm.EquivalentSources()
    拟合等效源
  • 将布格异常网格化到规则网格
  • 应用
    vd.distance_mask()
    对远离数据的区域进行掩膜
  • 应用导数滤波器（水平梯度、倾斜角）用于解释
  • 进行向上延拓以增强区域特征
  • 将网格化数据导出为NetCDF格式

Common Issues

常见问题

IssueSolution
Wrong gravity signCheck z-axis convention (positive upward)
Poor equivalent source fitAdjust 
depth
and 
damping
parameters
Slow terrain correctionReduce DEM resolution or use larger prisms
Edge effects in FFT filtersPad grid before applying 
upward_continuation
Coordinate mismatchEnsure consistent use of projected vs geographic coords
问题解决方案
重力符号错误检查z轴约定（向上为正）
等效源拟合效果差调整
depth
和
damping
参数
地形校正速度慢降低DEM分辨率或使用更大的棱柱体
FFT滤波器边缘效应在应用
upward_continuation
前对网格进行填充
坐标不匹配确保统一使用投影坐标或地理坐标

Tips

提示

  1. Use projected coordinates (meters) for local surveys
  2. Use tesseroids for regional/global scale modelling
  3. Equivalent sources handle irregular data spacing well
  4. Choose appropriate density (2670 kg/m3 typical for upper crust)
  5. Check sign conventions - depths are negative z values
  1. 使用投影坐标（米）处理局部测量数据
  2. 使用tesseroids处理区域/全球尺度建模
  3. 等效源能很好地处理不规则数据间距
  4. 选择合适的密度（上地壳典型值为2670千克/立方米）
  5. 检查符号约定 - 深度为z轴负值

References

参考资料

  • Gravity/Magnetic Corrections - Standard corrections and anomalies
  • Forward Modeling - Detailed forward modeling methods
  • 重力/磁校正 - 标准校正方法与异常类型
  • 正演模拟 - 详细的正演模拟方法

Scripts

脚本示例

  • scripts/gravity_processing.py - Process gravity survey data
  • scripts/gravity_processing.py - 重力测量数据处理脚本