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](references/corrections.md) - Standard corrections and anomalies
• [Forward Modeling](references/forward_modeling.md) - Detailed forward modeling methods

Scripts

• [scripts/gravityprocessing.py](scripts/gravityprocessing.py) - Process gravity survey data

Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

• Author: SteadfastAsArt
• Source: SteadfastAsArt/geoscience-skills
• License: MIT

Install and usage instructions live in the source repository linked above.