The fatiando package has been deprecated. Please check out the new tools in the Fatiando a Terra website:

GravMag: Calculate the gravity disturbance and Bouguer anomaly for HawaiiΒΆ

Download source code:

GravMag: Calculate the gravity disturbance and Bouguer anomaly for Hawaii
from fatiando.gravmag import normal_gravity
from fatiando.vis import mpl
import numpy as np
import urllib

# Download the gravity and topography data
url = ''
urllib.urlretrieve(url + 'eigen-6c3stat-havai.gdf',
urllib.urlretrieve(url + 'etopo1-havai.gdf',
# Load them with numpy
lon, lat, height, gravity = np.loadtxt('eigen-6c3stat-havai.gdf', skiprows=34,
topo = np.loadtxt('etopo1-havai.gdf', skiprows=30, usecols=[-1], unpack=True)
shape = (151, 151)
area = (lon.min(), lon.max(), lat.min(), lat.max())

# First, lets calculate the gravity disturbance (e.g., the free-air anomaly)
# We'll do this using the closed form of the normal gravity for the WGS84
# ellipsoid
gamma = normal_gravity.gamma_closed_form(lat, height)
disturbance = gravity - gamma

# Now we can remove the effect of the Bouguer plate to obtain the Bouguer
# anomaly. We'll use the standard densities of 2.67^-3 for crust and 1.04
#^-3 for water.
bouguer = disturbance - normal_gravity.bouguer_plate(topo)

mpl.figure(figsize=(14, 3.5))
bm = mpl.basemap(area, projection='merc')
mpl.title('Gravity (mGal)')
mpl.contourf(lon, lat, gravity, shape, 60,, basemap=bm)
mpl.title('Gravity disturbance (mGal)')
amp = np.abs(disturbance).max()
mpl.contourf(lon, lat, disturbance, shape, 60,, basemap=bm,
             vmin=-amp, vmax=amp)
mpl.title('Bouguer anomaly (mGal)')
mpl.contourf(lon, lat, bouguer, shape, 60,, basemap=bm)