harmonica

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_effect

Equivalent 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 angle

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.

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

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

References

• Gravity/Magnetic Corrections - Standard corrections and anomalies
• Forward Modeling - Detailed forward modeling methods

Scripts

• scripts/gravity_processing.py - Process gravity survey data