Gravity Disturbances
Note
Click here to download the full example code
Gravity Disturbances¶
Gravity disturbances are the differences between the measured gravity and
a reference (normal) gravity produced by an ellipsoid. The disturbances are
what allows geoscientists to infer the internal structure of the Earth. We’ll
use the boule.Ellipsoid.normal_gravity function from boule to
calculate the global gravity disturbance of the Earth using our sample gravity
data.
Out:
<xarray.Dataset>
Dimensions: (latitude: 361, longitude: 721)
Coordinates:
* longitude (longitude) float64 -180.0 -179.5 -179.0 ... 179.5 180.0
* latitude (latitude) float64 -90.0 -89.5 -89.0 ... 89.0 89.5 90.0
Data variables:
gravity (latitude, longitude) float64 9.801e+05 ... 9.802e+05
height_over_ell (latitude, longitude) float64 1e+04 1e+04 ... 1e+04 1e+04
Attributes: (12/35)
generating_institute: gfz-potsdam
generating_date: 2018/11/07
product_type: gravity_field
body: earth
modelname: EIGEN-6C4
max_used_degree: 1277
... ...
maxvalue: 9.8018358E+05 mgal
minvalue: 9.7476403E+05 mgal
signal_wrms: 1.5467865E+03 mgal
grid_format: long_lat_value
attributes: longitude latitude gravity_ell
attributes_units: deg. deg. mgal
import boule as bl
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import harmonica as hm
# Load the global gravity grid
data = hm.datasets.fetch_gravity_earth()
print(data)
# Calculate normal gravity using the WGS84 ellipsoid
ellipsoid = bl.WGS84
gamma = ellipsoid.normal_gravity(data.latitude, data.height_over_ell)
# The disturbance is the observed minus normal gravity (calculated at the
# observation point)
disturbance = data.gravity - gamma
# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=160))
pc = disturbance.plot.pcolormesh(
ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap="seismic"
)
plt.colorbar(
pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.5
)
ax.set_title("Gravity of disturbance of the Earth")
ax.coastlines()
plt.show()
Total running time of the script: ( 0 minutes 1.772 seconds)