Point sources¶

We can compute the gravitational field of point masses through the harmonica.point_gravity function. It offers the possibility to define point sources and computation points either in Cartesian or in spherical coordinates.

Each point mass can be defined as a tuple containing its coordinates in the following order: easting, northing and upward (in Cartesian coordinates) or longitude, spherical latitude and radius (in spherical coordinates).

Cartesian coordinates¶

For example, lets define a single point mass and compute the gravitational potential it generates on a computation point located 100 meters above it.

Define the point source and its mass (in kg):

point = (400, 300, 200)
mass = 1e6

Define a computation point located 100m above the point source:

coordinates = (point[0], point[1], point[2] + 100)

And finally, compute the gravitational potential field it generates on the computation point:

potential = hm.point_gravity(coordinates, point, mass, field="potential")
print(potential, "J/kg")

6.672999999999999e-07 J/kg

Working with multiple sources¶

We can easily compute the gravitational field of a set of point sources by grouping its coordinates in arrays and with a single call of the harmonica.point_gravity.

Lets define a set of four point sources, along with their masses:

import numpy as np

easting = np.array([250, 750, 250, 750])
northing = np.array([250, 250, 750, 750])
upward = np.array([-100, -100, -100, -100])
points = (easting, northing, upward)

masses = np.array([1e6, -1e6, 2e6, -3e6])

We can define a set of computation points located on a regular grid at zero height:

import verde as vd

coordinates = vd.grid_coordinates(
    region=(-250, 1250, -250, 1250), shape=(40, 40), extra_coords=0
)

And finally calculate the vertical component of the gravitational acceleration generated by the whole set of point sources on every computation point:

g_z = hm.point_gravity(coordinates, points, masses, field="g_z")

Note

When passing multiple sources and coordinates to harmonica.point_gravity we calculate the field in parallel using multiple CPUs, speeding up the computation.

Lets plot this gravitational field:

import matplotlib.pyplot as plt

maxabs = vd.maxabs(g_z)
plt.pcolormesh(
    *coordinates[:2], g_z, vmin=-maxabs, vmax=maxabs, cmap="seismic"
)
plt.colorbar(label="mGal")
plt.gca().set_aspect("equal")
plt.xlabel("easting [m]")
plt.ylabel("northing [m]")
plt.show()

Spherical coordinates¶

Alternatively, we can compute the gravitational fields of point sources defined in a spherical coordinate system. To do so, we need to pass the coordinate_system argument as "spherical"`. The coordinates of the point source must now be passed as longitude, spherical latitude and radius, where the two former ones must be in decimal degrees and the latter in meters.

Lets define a single point source in the equator, at a longitude of 45 degrees and located on the surface of the WGS84 ellipsoid.

We can start by obtaining the WGS84 ellipsoid from boule:

import boule as bl

ellipsoid = bl.WGS84

Then we can define the point source in the equator along with its mass:

longitude, latitude = 45, 0
radius = ellipsoid.geocentric_radius(latitude, geodetic=False)
point = (longitude, latitude, radius)

mass = 1e6

We can also define a computation point located 1km above the point source:

coordinates = (point[0], point[1], point[2] + 1000)

And finally compute the gravitational field generated by the source on the computation point:

g_z = hm.point_gravity(
    coordinates, point, mass, field="g_z", coordinate_system="spherical"
)
print(g_z, "mGal")

6.6729999999999995e-06 mGal

Working with sources defined in geodetic coordinates¶

If our point sources and computation points are defined in geodetic coordinates, we can use the boule.Ellipsoid.geodetic_to_spherical method to convert them to spherical coordinates and use them to compute their gravitational field.

Lets define a set point sources in geodetic coordinates and their masses:

longitude = np.array([-71, -71, -69, -69])
latitude = np.array([-45, -43, -45, -43])
height = np.array([-10e3, -20e3, -30e3, -20e3])
points = (longitude, latitude, height)

masses = np.array([1e6, 2e6, -3e6, 5e6])

Then, build a regular grid of computation points defined in geodetic coordinates and located at 1km above the ellipsoid:

coordinates = vd.grid_coordinates(
    region=(-72, -68, -46, -42),
    shape=(101, 101),
    extra_coords=20e3,
)

Before we can start forward modelling these sources, we need to convert them to spherical coordinates. To do so, we can use the boule.Ellipsoid.geodetic_to_spherical method:

points_spherical = ellipsoid.geodetic_to_spherical(*points)
coordinates_spherical = ellipsoid.geodetic_to_spherical(*coordinates)

We can finally use these converted coordinates to compute the gravitational field the source generate on every computation point:

g_z = hm.point_gravity(
    coordinates_spherical,
    points_spherical,
    masses,
    field="g_z",
    coordinate_system="spherical",
)

Lets plot these results using cartopy:

import cartopy.crs as ccrs

plt.figure(figsize=(8, 6))
ax = plt.axes(projection=ccrs.Mercator())
maxabs = vd.maxabs(g_z)
tmp = ax.pcolormesh(
    *coordinates[:2],
    g_z,
    vmin=-maxabs,
    vmax=maxabs,
    cmap="seismic",
    transform=ccrs.PlateCarree(),
)
ax.set_title("Gravitational acceleration (downward)")
plt.colorbar(tmp, ax=ax, pad=0.04, shrink=0.73, label="mGal")
plt.show()

Harmonica v0.5.0

Point sources

Contents