steadfastasart/geoscience-skills/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

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

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

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/gravity_processing.py](scripts/gravity_processing.py) - Process gravity survey data