Harmonica - Gravity and Magnetics
Quick Reference
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)Key Functions
| Function | Purpose |
|---|---|
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 |
Essential Operations
Forward Model - Rectangular Prism
# 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')Terrain Correction
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_effectEquivalent Source Gridding
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'])Magnetic Forward Model
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
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 angleCoordinate 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.
When to Use vs Alternatives
| Use Case | Tool | Why |
|---|---|---|
| Gravity/magnetic forward modelling | Harmonica | Purpose-built, Fatiando ecosystem |
| Potential field inversion | SimPEG | Full inversion framework with regularization |
| Commercial gravity processing | Oasis Montaj | Industry-standard GUI, proprietary formats |
| Simple Bouguer corrections only | Custom numpy | Fewer dependencies for one-off calculations |
| Equivalent source gridding | Harmonica | Best open-source option for potential fields |
| Regional/global scale | Harmonica (tesseroids) | Handles spherical geometry natively |
| Magnetic data reduction to pole | Harmonica | FFT-based filters for gridded data |
| Teaching/prototyping | Harmonica | Clean 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.
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
Common Issues
| Issue | Solution |
|---|---|
| Wrong gravity sign | Check z-axis convention (positive upward) |
| Poor equivalent source fit | Adjust depth and damping parameters |
| Slow terrain correction | Reduce DEM resolution or use larger prisms |
| Edge effects in FFT filters | Pad grid before applying upward_continuation |
| Coordinate mismatch | Ensure consistent use of projected vs geographic coords |
Tips
- Use projected coordinates (meters) for local surveys
- Use tesseroids for regional/global scale modelling
- Equivalent sources handle irregular data spacing well
- Choose appropriate density (2670 kg/m3 typical for upper crust)
- Check sign conventions - depths are negative z values
References
- Gravity/Magnetic Corrections - Standard corrections and anomalies
- Forward Modeling - Detailed forward modeling methods
Scripts
- scripts/gravity_processing.py - Process gravity survey data