In [1]:
# These two lines are necessary only if gempy is not installed
import sys, os
sys.path.append("../")
# Importing gempy
import gempy as gp
# Embedding matplotlib figures into the notebooks
%matplotlib inline
# Aux imports
import numpy as np
In this case instead loading a geo_data object directly, we will create one. The main atributes we need to pass are:
Additionaly we can pass the address to csv files (GeoModeller3D format) with the data.
In [2]:
# Importing the data from csv files and settign extent and resolution
geo_data = gp.create_data([696000,747000,6863000,6950000,-20000, 200],[50, 50, 50],
path_f = os.pardir+"/input_data/a_Foliations.csv",
path_i = os.pardir+"/input_data/a_Points.csv")
gp.get_raw_data(geo_data, 'interfaces').head()
Out[2]:
You can visualize the points in 3D (work in progress)
In [3]:
gp.visualize(geo_data)
Or a projection in 2D:
In [4]:
gp.plot_data(geo_data, direction='z')
Out[4]:
This model consist in 3 different depositional series. This mean that only data in the same depositional series affect the interpolation. To select with formations belong to witch series we will use the set_data_series function which takes a python dictionary as input.
We can see the unique formations with:
In [5]:
gp.get_series(geo_data)
Out[5]:
Setting the series we also give the specific order of the series. In python 3.6 and above the dictionaries conserve the key order so it is not necessary to give explicitly the order of the series.
Notice as well that the order of the formations within each series is not relevant for the result but in case of being wrong can lead to confusing color coding (work in progress).
In the representation given by get_series the elements get repeated but is only how Pandas print tables.
In [6]:
# Assigning series to formations as well as their order (timewise)
gp.set_data_series(geo_data, {"EarlyGranite_Series": 'EarlyGranite',
"BIF_Series":('SimpleMafic2', 'SimpleBIF'),
"SimpleMafic_Series":'SimpleMafic1'},
order_series = ["EarlyGranite_Series",
"BIF_Series",
"SimpleMafic_Series"], verbose=1)
Out[6]:
Now as in the previous chapter we just need to create the interpolator object and compute the model.
In [7]:
interp_data = gp.InterpolatorInput(geo_data)
In [8]:
sol = gp.compute_model(interp_data)
Now if we analyse the results we have a 3D array where the axis 0 represent the superposition of the series (potential fields). The color coding is working process yet.
In [9]:
import matplotlib.pyplot as plt
gp.plot_section(geo_data, sol[0,0,:], 11)
plt.show()
gp.plot_section(geo_data, sol[1,0,:], 11)
plt.show()
gp.plot_section(geo_data, sol[2,0,:], 11)
plt.show()
The axis 1 keeps the potential field:
In [10]:
gp.plot_potential_field(geo_data, sol[0,1,:], 11, cmap='inferno_r')
plt.show()
gp.plot_potential_field(geo_data, sol[1,1,:], 11, cmap='inferno_r')
plt.show()
gp.plot_potential_field(geo_data, sol[2,1,:], 11, cmap='inferno_r')
plt.show()
And the axis 2 keeps the faults network that in this model since there is not faults does not represent anything.
Additionally with can export the blocks to vtk in order to visualize them in Paraview. We are working in visualization in place as well.
In [11]:
gp.export_vtk_rectilinear(geo_data, sol[-1, 0, :], path=None)