Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"""
Earth Gravity
=============
This is the magnitude of the gravity vector of the Earth (gravitational
+ centrifugal) at 10 km height. The data is on a regular grid with 0.5 degree
spacing at 10km ellipsoidal height. It was generated from the spherical
harmonic model EIGEN-6C4 [Forste_etal2014]_.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm
# Load the gravity grid
data = hm.datasets.fetch_gravity_earth()
print(data)
# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=150))
pc = data.gravity.plot.pcolormesh(
ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False
)
plt.colorbar(
pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
)
ax.set_title("Gravity of the Earth (EIGEN-6C4)")
ax.coastlines()
plt.tight_layout()
plt.show()
====================
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 :meth:`boule.Ellipsoid.normal_gravity` function from :mod:`boule` to
calculate the global gravity disturbance of the Earth using our sample gravity
data.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import boule as bl
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(
gravity disturbance for Earth using our sample gravity and topography data. One
thing to note is that the ETOPO1 topography is referenced to the geoid, not the
ellipsoid. Since we want to remove the masses between the surface of the Earth
and ellipsoid, we need to add the geoid height to the topography before Bouguer
correction.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr
import boule as bl
import harmonica as hm
# Load the global gravity, topography, and geoid grids
data = xr.merge(
[
hm.datasets.fetch_gravity_earth(),
hm.datasets.fetch_geoid_earth(),
hm.datasets.fetch_topography_earth(),
]
)
print(data)
# Calculate normal gravity and the disturbance
ellipsoid = bl.WGS84
gamma = ellipsoid.normal_gravity(data.latitude, data.height_over_ell)
disturbance = data.gravity - gamma
# Reference the topography to the ellipsoid
topography_ell = data.topography + data.geoid
# Calculate the Bouguer planar correction and the topography-free disturbance.
# Use the default densities for the crust and ocean water.