# 2.2: Centered Grid.ΒΆ

Geophysics Preprocessing builds on the centered grid (https://github.com/cgre-aachen/gempy/blob/master/notebooks/tutorials/ch1-3-Grids.ipynb) to precompute the constant part of forward physical computations as for example gravity:

$F_z = G_{\rho} ||| x \ln(y+r) + y \ln (x+r) - z \arctan (\frac{x y}{z r}) |^{x_2}_{x_1}|^{y_2}_{y_1}|^{ z_2}_{z_1}$

where we can compress the grid dependent terms as

$t_z = ||| x \ln (y+r) + y \ln (x+r)-z \arctan ( \frac{x y}{z r} ) |^{x_2}_{x_1}|^{y_2}_{y_1}|^{z_2}_{z_1}$

By doing this decomposition an keeping the grid constant we can compute the forward gravity by simply operate:

$F_z = G_{\rho} \cdot t_z$
# Importing gempy
from gempy.assets.geophysics import GravityPreprocessing

# Aux imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(1515)
pd.set_option('precision', 2)

g = GravityPreprocessing()

kernel_centers, kernel_dxyz_left, kernel_dxyz_right = g.create_irregular_grid_kernel(resolution=[10, 10, 20],


create_irregular_grid_kernel will create a constant kernel around the point 0,0,0. This kernel will be what we use for each device.

kernel_centers


Out:

array([[-100.        , -100.        ,   -6.        ],
[-100.        , -100.        ,   -7.2       ],
[-100.        , -100.        ,   -7.52912998],
...,
[ 100.        ,  100.        ,  -79.90178533],
[ 100.        ,  100.        , -100.17119644],
[ 100.        ,  100.        , -126.        ]])


$$t_z$$ is only dependent on distance and therefore we can use the kerenel created on the previous cell

tz = g.set_tz_kernel(resolution=[10, 10, 20], radius=100)
tz


Out:

array([-8.71768928e-05, -6.45647022e-05, -3.41579985e-05, ...,
-1.09610058e-02, -1.41543038e-02, -1.51096613e-02])


To compute tz we also need the edges of each voxel. The distance to the edges are stored on kernel_dxyz_left and kernel_dxyz_right. We can plot all the data as follows:

a, b, c = kernel_centers, kernel_dxyz_left, kernel_dxyz_right

fig = plt.figure(figsize=(13, 7))
plt.quiver(a[:, 0].reshape(11, 11, 21)[5, :, :].ravel(),
a[:, 2].reshape(11, 11, 21)[:, 5, :].ravel(),
np.zeros(231),
tz.reshape(11, 11, 21)[5, :, :].ravel(), label='$t_z$', alpha=.3
)

plt.plot(a[:, 0].reshape(11, 11, 21)[5, :, :].ravel(),
a[:, 2].reshape(11, 11, 21)[:, 5, :].ravel(), 'o', alpha=.3, label='Centers')

plt.plot(a[:, 0].reshape(11, 11, 21)[5, :, :].ravel() - b[:, 0].reshape(11, 11, 21)[5, :, :].ravel(),
a[:, 2].reshape(11, 11, 21)[:, 5, :].ravel(), '.', alpha=.3, label='Lefts')

plt.plot(a[:, 0].reshape(11, 11, 21)[5, :, :].ravel(),
a[:, 2].reshape(11, 11, 21)[:, 5, :].ravel() - b[:, 2].reshape(11, 11, 21)[:, 5, :].ravel(), '.', alpha=.6,
label='Ups')

plt.plot(a[:, 0].reshape(11, 11, 21)[5, :, :].ravel() + c[:, 0].reshape(11, 11, 21)[5, :, :].ravel(),
a[:, 2].reshape(11, 11, 21)[:, 5, :].ravel(), '.', alpha=.3, label='Rights')

plt.plot(a[:, 0].reshape(11, 11, 21)[5, :, :].ravel(),
a[:, 2].reshape(11, 11, 21)[:, 5, :].ravel() + c[:, 2].reshape(11, 11, 21)[5, :, :].ravel(), '.', alpha=.3,
label='Downs')

plt.xlim(-200, 200)
plt.ylim(-200, 0)
plt.legend()
plt.show()


Just the quiver:

fig = plt.figure(figsize=(13, 7))
plt.quiver(a[:, 0].reshape(11, 11, 21)[5, :, :].ravel(),
a[:, 2].reshape(11, 11, 21)[:, 5, :].ravel(),
np.zeros(231),
tz.reshape(11, 11, 21)[5, :, :].ravel()
)
plt.show()


Remember this is happening always in 3D:

fig = plt.figure()

ax.scatter(a[:, 0], a[:, 1], a[:, 2], c=tz)

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()


Total running time of the script: ( 0 minutes 0.369 seconds)

Gallery generated by Sphinx-Gallery