harmonica

Original：🇺🇸 English

Translated

1 scriptsChecked / no sensitive code detected

Gravity and magnetic data processing and forward modelling using Fatiando a Terra. Use when Claude needs to: (1) Compute gravity forward models (point masses, prisms, tesseroids), (2) Apply terrain/Bouguer corrections, (3) Grid scattered potential field data with equivalent sources, (4) Perform upward/downward continuation, (5) Calculate magnetic anomalies from magnetized bodies, (6) Apply derivative filters (gradients, tilt angle), (7) Process regional or local gravity surveys.

2installs

Sourcesteadfastasart/geoscience-skills

Added on2026-05-16

NPX Install

npx skill4agent add steadfastasart/geoscience-skills harmonica

SKILL.md Content

View Translation Comparison →

Harmonica - Gravity and Magnetics

Quick Reference

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

# 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

python

# 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

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

# 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

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

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

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

harmonica

NPX Install

Tags

SKILL.md Content

Harmonica - Gravity and Magnetics

Quick Reference

Key Functions

Essential Operations

Forward Model - Rectangular Prism

Terrain Correction

Equivalent Source Gridding

Magnetic Forward Model

Derivative Filters

Coordinate System

When to Use vs Alternatives

Common Workflows

Process Gravity Survey with Terrain Correction and Gridding

Common Issues

Tips

References

Scripts