2.1 Forward Gravity: Simple example

Importing gempy

import gempy as gp
from gempy.assets.geophysics import GravityPreprocessing

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

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

cwd = os.getcwd()
if 'examples' not in cwd:
    data_path = os.getcwd() + '/examples/'
else:
    data_path = cwd + '/../../'

geo_model = gp.load_model('Greenstone', path=data_path + 'data/gempy_models/Greenstone')

Out:

Active grids: ['regular']

geo_model.stack

	order_series	BottomRelation	isActive	isFault	isFinite
EarlyGranite_Series	1	Erosion	True	False	False
BIF_Series	2	Erosion	True	False	False
SimpleMafic_Series	3	Erosion	True	False	False
Basement	4	Erosion	False	False	False

geo_model.surfaces

	surface	series	order_surfaces	color	id	value_0
3	EarlyGranite	EarlyGranite_Series	1	#728f02	1	2.61
0	SimpleMafic2	BIF_Series	1	#015482	2	2.92
1	SimpleBIF	BIF_Series	2	#9f0052	3	3.10
2	SimpleMafic1	SimpleMafic_Series	1	#ffbe00	4	2.92
4	basement	Basement	1	#443988	5	2.61

gp.plot_2d(geo_model)

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7fcc4670beb0>

Creating grid

First we need to define the location of the devices. For this example we can make a map:

grav_res = 20
X = np.linspace(7.050000e+05, 747000, grav_res)
Y = np.linspace(6863000, 6925000, grav_res)
Z = 300
xyz = np.meshgrid(X, Y, Z)
xy_ravel = np.vstack(list(map(np.ravel, xyz))).T
xy_ravel

Out:

array([[7.05000000e+05, 6.86300000e+06, 3.00000000e+02],
       [7.07210526e+05, 6.86300000e+06, 3.00000000e+02],
       [7.09421053e+05, 6.86300000e+06, 3.00000000e+02],
       ...,
       [7.42578947e+05, 6.92500000e+06, 3.00000000e+02],
       [7.44789474e+05, 6.92500000e+06, 3.00000000e+02],
       [7.47000000e+05, 6.92500000e+06, 3.00000000e+02]])

We can see the location of the devices relative to the model data:

gp.plot_2d(geo_model, direction='z', show=False)
plt.scatter(xy_ravel[:, 0], xy_ravel[:, 1], s=1)
plt.show()

Now we need to create the grid centered on the devices (see: https://github.com/cgre-aachen/gempy/blob/master/notebooks/tutorials/ch1-3-Grids.ipynb)

geo_model.set_centered_grid(xy_ravel, resolution=[10, 10, 15], radius=5000)

Out:

Active grids: ['regular' 'centered']

Grid Object. Values:
array([[ 6.96510000e+05,  6.86367000e+06, -1.97980000e+04],
       [ 6.96510000e+05,  6.86367000e+06, -1.93940000e+04],
       [ 6.96510000e+05,  6.86367000e+06, -1.89900000e+04],
       ...,
       [ 7.52000000e+05,  6.93000000e+06, -3.10768481e+03],
       [ 7.52000000e+05,  6.93000000e+06, -4.31811404e+03],
       [ 7.52000000e+05,  6.93000000e+06, -6.00000000e+03]])

geo_model.grid.centered_grid.kernel_centers

Out:

array([[-5000.        , -5000.        ,  -300.        ],
       [-5000.        , -5000.        ,  -360.        ],
       [-5000.        , -5000.        ,  -383.36972966],
       ...,
       [ 5000.        ,  5000.        , -3407.68480754],
       [ 5000.        ,  5000.        , -4618.11403801],
       [ 5000.        ,  5000.        , -6300.        ]])

Now we need to compute the component tz (see https://github.com/cgre-achen/gempy/blob/master/notebooks/tutorials/ch2-2-Cell_selection.ipynb)

g = GravityPreprocessing(geo_model.grid.centered_grid)

tz = g.set_tz_kernel()

tz

Out:

array([-0.00435884, -0.0035374 , -0.00260207, ..., -0.60455378,
       -0.888396  , -0.98280245])

Compiling the gravity graph

If geo_model has already a centered grid, the calculation of tz happens automatically. This theano graph will return gravity as well as the lithologies. In addition we need either to pass the density block (see below). Or the position of density on the surface(in the future the name) to compute the density block at running time.

geo_model.surfaces

	surface	series	order_surfaces	color	id	value_0
3	EarlyGranite	EarlyGranite_Series	1	#728f02	1	2.61
0	SimpleMafic2	BIF_Series	1	#015482	2	2.92
1	SimpleBIF	BIF_Series	2	#9f0052	3	3.10
2	SimpleMafic1	SimpleMafic_Series	1	#ffbe00	4	2.92
4	basement	Basement	1	#443988	5	2.61

In this case the densities of each layer are at the loc 1 (0 is the id)

# New way
gp.set_interpolator(geo_model, output=['gravity'], pos_density=1, gradient=False,
                    theano_optimizer='fast_run')

Out:

Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization:  fast_run
Device:  cpu
Precision:  float64
Number of faults:  0
Compilation Done!
Kriging values:
                        values
range                86591.22
$C_o$             178524761.9
drift equations  [3, 3, 3, 3]
nugget grad              0.01
nugget scalar             0.0

<gempy.core.interpolator.InterpolatorModel object at 0x7fcc46f045e0>

Once we have created a gravity interpolator we can call it from compute model as follows:

sol = gp.compute_model(geo_model)
grav = sol.fw_gravity

gp.plot_2d(geo_model, direction=['z'], height=7, show_results=False, show_data=True,
           show=False)
plt.scatter(xy_ravel[:, 0], xy_ravel[:, 1], s=1)
plt.imshow(sol.fw_gravity.reshape(grav_res, grav_res),
           extent=(xy_ravel[:, 0].min() + (xy_ravel[0, 0] - xy_ravel[1, 0]) / 2,
                   xy_ravel[:, 0].max() - (xy_ravel[0, 0] - xy_ravel[1, 0]) / 2,
                   xy_ravel[:, 1].min() + (xy_ravel[0, 1] - xy_ravel[30, 1]) / 2,
                   xy_ravel[:, 1].max() - (xy_ravel[0, 1] - xy_ravel[30, 1]) / 2),
           cmap='viridis_r', origin='lower')
plt.show()

Plotting lithologies

If we want to compute the lithologies we will need to create a normal interpolator object as seen in the Chapter 1 of the tutorials

Now we can plot all together (change the alpha parameter to see the gravity overlying):

sphinx_gallery_thumbnail_number = 4

gp.plot_2d(geo_model, cell_number=[-1], direction=['z'], show=False,
           kwargs_regular_grid={'alpha': .5})

plt.scatter(xy_ravel[:, 0], xy_ravel[:, 1], s=1)
plt.imshow(grav.reshape(grav_res, grav_res),
           extent=(xy_ravel[:, 0].min() + (xy_ravel[0, 0] - xy_ravel[1, 0]) / 2,
                   xy_ravel[:, 0].max() - (xy_ravel[0, 0] - xy_ravel[1, 0]) / 2,
                   xy_ravel[:, 1].min() + (xy_ravel[0, 1] - xy_ravel[30, 1]) / 2,
                   xy_ravel[:, 1].max() - (xy_ravel[0, 1] - xy_ravel[30, 1]) / 2),
           cmap='viridis_r', origin='lower', alpha=.8)
cbar = plt.colorbar()
cbar.set_label(r'$\mu$gal')
plt.show()