Example 1: Sandstone Model


In [1]:
# Importing
import theano.tensor as T
import sys, os
sys.path.append("../GeMpy")

# Importing GeMpy modules
import GeMpy_core
import Visualization

# Reloading (only for development purposes)
import importlib
importlib.reload(GeMpy_core)
importlib.reload(Visualization)

# Usuful packages
import numpy as np
import pandas as pn

import matplotlib.pyplot as plt

# This was to choose the gpu
os.environ['CUDA_LAUNCH_BLOCKING'] = '1'

# Default options of printin
np.set_printoptions(precision = 6, linewidth= 130, suppress =  True)

%matplotlib inline
#%matplotlib notebook

First we make a GeMpy instance with most of the parameters default (except range that is given by the project). Then we also fix the extension and the resolution of the domain we want to interpolate. Finally we compile the function, only needed once every time we open the project (the guys of theano they are working on letting loading compiled files, even though in our case it is not a big deal).

General note. So far the reescaling factor is calculated for all series at the same time. GeoModeller does it individually for every potential field. I have to look better what this parameter exactly means


In [2]:
# Setting extend, grid and compile
# Setting the extent
sandstone = GeMpy_core.GeMpy()

# Create Data class with raw data
sandstone.import_data( [696000,747000,6863000,6950000,-20000, 2000],[ 40, 40, 80],
                         path_f = os.pardir+"/input_data/a_Foliations.csv",
                         path_i = os.pardir+"/input_data/a_Points.csv")

All input data is stored in pandas dataframes under, self.Data.Interances and self.Data.Foliations:


In [3]:
sandstone.Data.Foliations.head()


Out[3]:
X Y Z azimuth dip polarity formation series G_x G_y G_z
0 739426.627684 6.891935e+06 75.422691 220.000 70.0 1 SimpleBIF Default serie -0.604023 -7.198463e-01 0.342020
1 717311.112372 6.891941e+06 -1497.488309 90.000 60.0 1 SimpleBIF Default serie 0.866025 5.302876e-17 0.500000
2 723415.321182 6.891939e+06 -5298.154309 10.000 20.0 1 SimpleMafic1 Default serie 0.059391 3.368241e-01 0.939693
3 742907.686575 6.891935e+06 -2786.869309 250.000 60.0 1 SimpleMafic1 Default serie -0.813798 -2.961981e-01 0.500000
4 712584.536312 6.891942e+06 -582.769334 90.014 60.0 1 SimpleMafic1 Default serie 0.866025 -2.116099e-04 0.500000

In case of disconformities, we can define which formation belong to which series using a dictionary. Until python 3.6 is important to specify the order of the series otherwise is random


In [4]:
sandstone.Data.set_series({"EarlyGranite_Series":sandstone.Data.formations[-1], 
                      "BIF_Series":(sandstone.Data.formations[0], sandstone.Data.formations[1]),
                      "SimpleMafic_Series":sandstone.Data.formations[2]}, 
                       order = ["EarlyGranite_Series",
                              "BIF_Series",
                              "SimpleMafic_Series"])


Out[4]:
EarlyGranite_Series BIF_Series SimpleMafic_Series
0 EarlyGranite SimpleMafic2 SimpleMafic1
1 EarlyGranite SimpleBIF SimpleMafic1

Now in the data frame we should have the series column too


In [5]:
sandstone.Data.Foliations.head()


Out[5]:
X Y Z azimuth dip polarity formation series G_x G_y G_z
0 739426.627684 6.891935e+06 75.422691 220.000 70.0 1 SimpleBIF BIF_Series -0.604023 -7.198463e-01 0.342020
1 717311.112372 6.891941e+06 -1497.488309 90.000 60.0 1 SimpleBIF BIF_Series 0.866025 5.302876e-17 0.500000
2 723415.321182 6.891939e+06 -5298.154309 10.000 20.0 1 SimpleMafic1 SimpleMafic_Series 0.059391 3.368241e-01 0.939693
3 742907.686575 6.891935e+06 -2786.869309 250.000 60.0 1 SimpleMafic1 SimpleMafic_Series -0.813798 -2.961981e-01 0.500000
4 712584.536312 6.891942e+06 -582.769334 90.014 60.0 1 SimpleMafic1 SimpleMafic_Series 0.866025 -2.116099e-04 0.500000

Next step is the creating of a grid. So far only regular. By default it takes the extent and the resolution given in the import_data method.


In [6]:
# Create a class Grid so far just regular grid
sandstone.create_grid()
sandstone.Grid.grid


Out[6]:
array([[  696000.      ,  6863000.      ,   -20000.      ],
       [  696000.      ,  6863000.      ,   -19721.519531],
       [  696000.      ,  6863000.      ,   -19443.037109],
       ..., 
       [  747000.      ,  6950000.      ,     1443.037964],
       [  747000.      ,  6950000.      ,     1721.519043],
       [  747000.      ,  6950000.      ,     2000.      ]], dtype=float32)

Plotting raw data

The object Plot is created automatically as we call the methods above. This object contains some methods to plot the data and the results.

It is possible to plot a 2D projection of the data in a specific direction using the following method. Also is possible to choose the series you want to plot. Additionally all the key arguments of seaborn lmplot can be used.


In [7]:
sandstone.Plot.plot_data(series = sandstone.Data.series.columns.values[1])


Class Interpolator

This class will take the data from the class Data and calculate potential fields and block. We can pass as key arguments all the variables of the interpolation. I recommend not to touch them if you do not know what are you doing. The default values should be good enough. Also the first time we execute the method, we will compile the theano function so it can take a bit of time.


In [8]:
sandstone.set_interpolator()

Now we could visualize the individual potential fields as follow:

Early granite


In [9]:
sandstone.Plot.plot_potential_field(10, n_pf=0)


BIF Series


In [10]:
sandstone.Plot.plot_potential_field(13, n_pf=1, cmap = "magma",  plot_data = True,
                                        verbose = 5 )


SImple mafic


In [11]:
sandstone.Plot.plot_potential_field(10, n_pf=2)


Optimizing the export of lithologies

But usually the final result we want to get is the final block. The method compute_block_model will compute the block model, updating the attribute block. This attribute is a theano shared function that can return a 3D array (raveled) using the method get_value().


In [12]:
# Reset the block
sandstone.Interpolator.block.set_value(np.zeros_like(sandstone.Grid.grid[:,0]))

# Compute the block
sandstone.Interpolator.compute_block_model([0,1,2], verbose = 0)

In [13]:
sandstone.Interpolator.block.get_value(), np.unique(sandstone.Interpolator.block.get_value())


Out[13]:
(array([ 0.,  0.,  0., ...,  0.,  0.,  0.], dtype=float32),
 array([ 0.,  1.,  2.,  3.,  4.], dtype=float32))

And again after computing the model in the Plot object we can use the method plot_block_section to see a 2D section of the model


In [14]:
sandstone.Plot.plot_block_section(13, interpolation = 'nearest',  direction='y')
plt.savefig("sandstone_example.png")


Export to vtk. (Under development)


In [16]:
"""Export model to VTK

Export the geology blocks to VTK for visualisation of the entire 3-D model in an
external VTK viewer, e.g. Paraview.

..Note:: Requires pyevtk, available for free on: https://github.com/firedrakeproject/firedrake/tree/master/python/evtk

**Optional keywords**:
    - *vtk_filename* = string : filename of VTK file (default: output_name)
    - *data* = np.array : data array to export to VKT (default: entire block model)
"""
vtk_filename = "noddyFunct2"

extent_x = 10
extent_y = 10
extent_z = 10

delx = 0.2
dely = 0.2
delz = 0.2
from pyevtk.hl import gridToVTK
# Coordinates
x = np.arange(0, extent_x + 0.1*delx, delx, dtype='float64')
y = np.arange(0, extent_y + 0.1*dely, dely, dtype='float64')
z = np.arange(0, extent_z + 0.1*delz, delz, dtype='float64')

# self.block = np.swapaxes(self.block, 0, 2)
gridToVTK(vtk_filename, x, y, z, cellData = {"geology" : sol})


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-16-ff637538da86> in <module>()
     26 
     27 # self.block = np.swapaxes(self.block, 0, 2)
---> 28 gridToVTK(vtk_filename, x, y, z, cellData = {"geology" : sol})

NameError: name 'sol' is not defined

Performance Analysis

One of the advantages of theano is the posibility to create a full profile of the function. This has to be included in at the time of the creation of the function. At the moment it should be active (the downside is larger compilation time and I think also a bit in the computation so be careful if you need a fast call)

CPU

The following profile is with a 2 core laptop. Nothing spectacular.


In [ ]:
%%timeit
# Reset the block
sandstone.Interpolator.block.set_value(np.zeros_like(sandstone.Grid.grid[:,0]))

# Compute the block
sandstone.Interpolator.compute_block_model([0,1,2], verbose = 0)

Looking at the profile we can see that most of time is in pow operation (exponential). This probably is that the extent is huge and we are doing it with too much precision. I am working on it


In [ ]:
esandstone.Interpolator._interpolate.profile.summary()

GPU


In [ ]:
%%timeit
# Reset the block
sandstone.block.set_value(np.zeros_like(sandstone.grid[:,0]))

# Compute the block
sandstone.compute_block_model([0,1,2], verbose = 0)

In [ ]:
sandstone.block_export.profile.summary()

In [ ]: