In [1]:
# Importing
import theano.tensor as T
import theano
import sys, os
sys.path.append("../GeMpy")
# Importing GeMpy modules
import GeMpy
# Reloading (only for development purposes)
import importlib
importlib.reload(GeMpy)
# 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 inline
In [2]:
# Importing the data from csv files and settign extent and resolution
geo_data = GeMpy.import_data([696000,747000,6863000,6950000,-20000, 2000],[ 40, 40, 40],
path_f = os.pardir+"/input_data/a_Foliations.csv",
path_i = os.pardir+"/input_data/a_Points.csv")
# Assigning series to formations as well as their order (timewise)
GeMpy.set_data_series(geo_data, {"EarlyGranite_Series":geo_data.formations[-1],
"BIF_Series":(geo_data.formations[0], geo_data.formations[1]),
"SimpleMafic_Series":geo_data.formations[2]},
order_series = ["EarlyGranite_Series",
"BIF_Series",
"SimpleMafic_Series"], verbose=0)
In [4]:
# Preprocessing data to interpolate: This rescales the coordinates between 0 and 1 for stability issues. Here we can choose also the drift degree
# (in new updates I will change it to be possible to change the grade after compilation). From here we can set also the data type of the operations in case
# you want to use the GPU. Verbose is huge. There is a large list of strings that select what you want to print while the computation.
data_interp = GeMpy.set_interpolator(geo_data, u_grade = 2, dtype="float32", verbose=[])
In [5]:
# This cell will go to the backend
# Set all the theano shared parameters and return the symbolic variables (the input of the theano function)
input_data_T = data_interp.interpolator.tg.input_parameters_list()
# Prepare the input data (interfaces, foliations data) to call the theano function. Also set a few theano shared variables with the len
# of formations series and so on
input_data_P = data_interp.interpolator.data_prep()
# Compile the theano function.
debugging = theano.function(input_data_T, data_interp.interpolator.tg.whole_block_model(), on_unused_input='ignore',
allow_input_downcast=True, profile=True)
In [6]:
# Solve model calling the theano function
sol = debugging(input_data_P[0], input_data_P[1], input_data_P[2], input_data_P[3],input_data_P[4], input_data_P[5])
In [7]:
# Plot the block model.
GeMpy.plot_section(geo_data, 13, block = sol, direction='y', plot_data = False)
Out[7]:
In [11]:
p3.figure
In [12]:
# So far this is a simple 3D visualization. I have to adapt it into GeMpy
lith0 = sol == 0
lith1 = sol == 1
lith2 = sol == 2
lith3 = sol == 3
lith4 = sol == 4
np.unique(sol)
import ipyvolume.pylab as p3
p3.figure(width=800)
p3.scatter(geo_data.grid.grid[:,0][lith0],
geo_data.grid.grid[:,1][lith0],
geo_data.grid.grid[:,2][lith0], marker='box', color = 'blue', size = 0.1 )
p3.scatter(geo_data.grid.grid[:,0][lith1],
geo_data.grid.grid[:,1][lith1],
geo_data.grid.grid[:,2][lith1], marker='box', color = 'yellow', size = 1 )
p3.scatter(geo_data.grid.grid[:,0][lith2],
geo_data.grid.grid[:,1][lith2],
geo_data.grid.grid[:,2][lith2], marker='box', color = 'green', size = 1 )
p3.scatter(geo_data.grid.grid[:,0][lith3],
geo_data.grid.grid[:,1][lith3],
geo_data.grid.grid[:,2][lith3], marker='box', color = 'pink', size = 1 )
p3.scatter(geo_data.grid.grid[:,0][lith4],
geo_data.grid.grid[:,1][lith4],
geo_data.grid.grid[:,2][lith4], marker='box', color = 'red', size = 1 )
p3.xlim(np.min(geo_data.grid.grid[:,0]),np.min(geo_data.grid.grid[:,0])+2175.0*40)
p3.ylim(np.min(geo_data.grid.grid[:,1]),np.max(geo_data.grid.grid[:,1]))
p3.zlim(np.min(geo_data.grid.grid[:,2]),np.min(geo_data.grid.grid[:,2])+2175.0*40)#np.max(geo_data.grid.grid[:,2]))
p3.show()
In [11]:
# The profile at the moment sucks because all what is whithin a scan is not subdivided
debugging.profile.summary()
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
All input data is stored in pandas dataframes under, self.Data.Interances and self.Data.Foliations:
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
Now in the data frame we should have the series column too
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
#GeMpy.set_grid(geo_data)
#GeMpy.get_grid(geo_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]:
#GeMpy.plot_data(geo_data, 'y', geo_data.series.columns.values[1])
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 [10]:
%debug
In [9]:
geo_data.interpolator.results
Out[9]:
In [10]:
geo_data.interpolator.tg.c_o_T.get_value(), geo_data.interpolator.tg.a_T.get_value()
Out[10]:
In [9]:
geo_data.interpolator.compile_potential_field_function()
Out[9]:
In [31]:
geo_data.interpolator.compute_potential_fields('BIF_Series',verbose = 3)
Out[31]:
In [18]:
geo_data.interpolator.potential_fields
Out[18]:
In [11]:
geo_data.interpolator.results
Out[11]:
In [41]:
geo_data.interpolator.tg.c_resc.get_value()
Out[41]:
Now we could visualize the individual potential fields as follow:
In [ ]:
GeMpy.plot_potential_field(geo_data,10, n_pf=0)
In [12]:
GeMpy.plot_potential_field(geo_data,13, n_pf=1, cmap = "magma", plot_data = True,
verbose = 5)
In [13]:
GeMpy.plot_potential_field(geo_data, 10, n_pf=2)
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 [7]:
GeMpy.compute_block_model(geo_data)
Out[7]:
In [ ]:
#GeMpy.set_interpolator(geo_data, u_grade = 0, compute_potential_field=True)
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 [8]:
GeMpy.plot_section(geo_data, 13, direction='y')
Out[8]:
In [14]:
"""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})
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)
The following profile is with a 2 core laptop. Nothing spectacular.
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 [16]:
%%timeit
# Compute the block
GeMpy.compute_block_model(geo_data, [0,1,2], verbose = 0)
In [17]:
geo_data.interpolator._interpolate.profile.summary()
In [ ]: