Gempy: critical points

Reference: https://docs.gempy.org/getting_started/get_started.html#sphx-glr-download-getting-started-get-started-py

*GemPy core code is written in Python. However for efficiency (and other reasons) most of heavy computations happend in optimize compile code, either C or CUDA for GPU. To do so, GemPy rely on the library theano. To guarantee maximum optimization theano requires to compile the code for every Python kernel.

*GemPy is a surface based interpolator. This means that all the input data we add has to be refered to a surface. The surfaces always mark the bottom of a unit. By default GemPy surfaces are empty:

*It is decisive to remember that, in GemPy, interface position points mark the bottom of a layer. 

*GemPy input data consist on surface points and orientations (perpendicular to the layers)

 

*As we generate our Data from CSV-files, we also have to define our model's real extent in 𝑥, 𝑦 and 𝑧, as well as declare a desired resolution for each axis. This resolution will in turn determine the number of voxels used during modeling. Here, we rely on a medium resolution of 50x50x50, amounting to 125,000 voxels. The model extent should be chosen in a way that it contains all relevant data in a representative space. As our model voxels are not cubes, but prisms, the resolution can take a different shape than the extent. We don't recommend going much higher than 100 cells in every direction (1,000,000 voxels), as higher resolutions will become increasingly expensive to compute.

 

*Every fault is treated as an independent series and have to be at set at the top of the pile. The relative order between the distinct faults defines the tectonic relation between them (first entry is the youngest).

 

Command lines

*gp.set_interpolator(geo_model, theano_optimizer='fast_compile', verbose=[])

: Optimize compile code

 

*gp.get_data() : list the input data

*gp.set_interpolator() : This function rescales the extent and coordinates of the original data (and store it in the attribute geo_data_res which behaves as a usual InputData object) and adds mathematical parameters that are needed for conducting the interpolation.

*gp.plot_2d(geo, model, cell_number = ?, direction='z')

e.g.) gp.plot_2d(geo_model, show_data=False, show_scalar=True, show_lith=False, cell_number=25, direction='y')
plt.show()

With cell_number=25 and remembering that we defined our resolution to be 50 cells in each direction, we have chosen a section going through the middle of our block. We have moved 25 cells in direction='y', the plot thus depicts a plane parallel to the 𝑥- and 𝑦-axes. Setting plot_data=True, we could plot original data together with the results. Changing the values for cell_number and direction, we can move through our 3D block model and explore it by looking at different 2D planes.

Pole_vector (np.ndarray) – [s4] (numpy.ndarray[float, 3]): 2D numpy array where axis 1 is the gradient values G_x, G_y, G_z of the pole while axis 0 is n number of orientations

Orientation (np.ndarray) – [s5] (numpy.ndarray[float, 3]): 2D numpy array where axis 1 is are orientation values [dip, azimuth, polarity] of the pole while axis 0 is n number of orientations. —  

Dip is the inclination angle of 0 to 90 degrees measured from the horizontal plane downwards.  

Azimuth is the dip direction defined by a 360 degrees clockwise rotation, i.e. 0 = North, 90 = East, 180 = South, and 270 = West.

Polarity defines where the upper (geologically younger) side of the orientation plane is and can be declared to be either normal (1) or reversed (-1).

The orientation plane is perpendicular to the gradient