Rectangular prisms

Rectangular prisms

One of the geometric bodies that Harmonica is able to forward model are rectangular prisms. We can compute their gravitational fields through the harmonica.prism_gravity function.

Rectangular prisms can only be defined in Cartesian coordinates, therefore they are intended to be used to forward model bodies in a small geographic region, in which we can neglect the curvature of the Earth.

Each prism can be defined as a tuple containing its boundaries in the following order: west, east, south, north, bottom, top in Cartesian coordinates and in meters.

Lets define a single prism and compute the gravitational potential it generates on a computation point located at 500 meters above its uppermost surface.

Define the prism and its density (in kg per cubic meters):

prism = (-100, 100, -100, 100, -1000, 500)
density = 3000

Define a computation point located 500 meters above the top surface of the prism:

coordinates = (0, 0, 1000)

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

potential = hm.prism_gravity(coordinates, prism, density, field="potential")
print(potential, "J/kg")
0.011051569953086739 J/kg

Passing multiple prisms

We can compute the gravitational field of a set of prisms by passing a list of them, where each prism is defined as mentioned before, and then making a single call of the harmonica.prism_gravity function.

Lets define a set of four prisms, along with their densities:

prisms = [
    [2e3, 3e3, 2e3, 3e3, -10e3, -1e3],
    [3e3, 4e3, 7e3, 8e3, -9e3, -1e3],
    [7e3, 8e3, 1e3, 2e3, -7e3, -1e3],
    [8e3, 9e3, 6e3, 7e3, -8e3, -1e3],
densities = [2670, 3300, 2900, 2980]

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=(0, 10e3, 0, 10e3), shape=(40, 40), extra_coords=0

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

g_z = hm.prism_gravity(coordinates, prisms, densities, field="g_z")


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

Lets plot this gravitational field:

import pygmt
grid = vd.make_xarray_grid(
   coordinates, g_z, data_names="g_z", extra_coords_names="extra")
fig = pygmt.Figure()
   region=(0, 10e3, 0, 10e3),
   frame=["WSne", "x+leasting (m)", "y+lnorthing (m)"],
fig.colorbar(cmap=True, position="JMR", frame=["a2", "x+lmGal"])

Prism layer

One of the most common usage of prisms is to model geologic structures. Harmonica offers the possibility to define a layer of prisms through the harmonica.prism_layer function: a regular grid of prisms of equal size along the horizontal dimensions and with variable top and bottom boundaries. It returns a xarray.Dataset with the coordinates of the centers of the prisms and their corresponding physical properties.

The harmonica.DatasetAccessorPrismsLayer Dataset accessor can be used to obtain some properties of the layer like its shape and size or the boundaries of any prism in the layer. Moreover, we can use the harmonica.DatasetAcessorPrismsLayer.gravity method to compute the gravitational field of the prism layer on any set of computation points.

Lets create a simple prism layer, whose lowermost boundaries will be set on zero and their uppermost boundary will approximate a sinusoidal function. We can start by setting the region of the layer and the horizontal dimensions of the prisms:

region = (0, 100e3, -40e3, 40e3)
spacing = 2000

Then we can define a regular grid where the centers of the prisms will fall:

easting, northing = vd.grid_coordinates(region=region, spacing=spacing)

We need to define a 2D array for the uppermost surface of the layer. We will sample a trigonometric function for this simple example:

import numpy as np

wavelength = 24 * spacing
surface = np.abs(np.sin(easting * 2 * np.pi / wavelength))

Lets assign the same density to each prism through a 2d array with the same value: 2700 kg per cubic meter.

density = np.full_like(surface, 2700)

Now we can define the prism layer specifying the reference level to zero:

prisms = hm.prism_layer(
    coordinates=(easting, northing),
    properties={"density": density},

Lets define a grid of observation points at 1 km above the zeroth height:

region_pad = vd.pad_region(region, 10e3)
coordinates = vd.grid_coordinates(
    region_pad, spacing=spacing, extra_coords=1e3

And compute the gravitational field generated by the prism layer on them:

gravity = prisms.prism_layer.gravity(coordinates, field="g_z")

Finally, lets plot the gravitational field:

grid = vd.make_xarray_grid(
   coordinates, gravity, data_names="gravity", extra_coords_names="extra")

fig = pygmt.Figure()
title = "Gravitational acceleration of a layer of prisms"
   frame=[f"WSne+t{title}", "x+leasting (m)", "y+lnorthing (m)"],
fig.colorbar(cmap=True, position="JMR", frame=["a.02", "x+lmGal"])