In [1]:
# Importing
import theano.tensor as T
import sys, os
sys.path.append("../GeMpy")
sys.path.append("../pygeomod")
import GeMpy_core
import Visualization
import importlib
importlib.reload(GeMpy_core)
importlib.reload(Visualization)
import numpy as np
import pandas as pn
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
os.environ['CUDA_LAUNCH_BLOCKING'] = '1'
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 [ ]:
# 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,
path_f = os.pardir+"/input_data/a_Foliations.csv",
path_i = os.pardir+"/input_data/a_Points.csv")
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"])
All input data is stored in pandas dataframes under, self.Data.Interances and self.Data.Foliations:
In [22]:
sandstone.Data.Foliations;
In [23]:
# Create a class Grid so far just regular grid
sandstone.create_grid()
sandstone.Grid.grid
Out[23]:
In [19]:
sandstone.Plot.plot_data(serie = sandstone.Data.series.columns.values[1])
This class will take the data from the class Data and calculate potential fields and block
In [24]:
sandstone.set_interpolator()
In [26]:
sandstone.Plot.plot_potential_field(10, n_pf=0)
In [31]:
sandstone.Plot.plot_potential_field(10, n_pf=1, cmap = "magma", plot_data = False)
In [29]:
sandstone.Plot.plot_potential_field(10, n_pf=2)
In [34]:
# 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 = 1)
In [35]:
sandstone.Plot.plot_block_section()
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})
In [32]:
%%timeit
sol = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0]
In [172]:
sandstone.block_export.profile.summary()
In [20]:
%%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 [22]:
sandstone.block_export.profile.summary()
In [ ]: