Code

Gempy Front

This file is part of gempy.

gempy is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

Foobar is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with gempy. If not, see <http://www.gnu.org/licenses/>.

Module with classes and methods to perform implicit regional modelling based on the potential field method. Tested on Ubuntu 16

Created on 10/10 /2016

@author: Miguel de la Varga

gempy_front.compute_model(interp_data, output=’geology’, u_grade=None, get_potential_at_interfaces=False)[source]

Computes the geological model and any extra output given.

Parameters:
  • interp_data (gempy.interpolator.InterpolatorData) –
  • output ({'geology', 'gravity', 'gradients'}) – Only if theano functions has not been compiled yet
  • u_grade (array-like of {0, 1, 2}) – grade of the polynomial for the universal part of the Kriging interpolations. The value has to be either 0, 1 or 2 and the length has to be the number of series. By default the value depends on the number of points given as input to try to avoid singular matrix. NOTE: if during the computation of the model a singular matrix is returned try to reduce the u_grade of the series.
  • get_potential_at_interfaces (bool) – Get potential at interfaces
Returns:

depending on the chosen out returns different number of solutions:
if output is geology:
  1. Lithologies: block and scalar field
  2. Faults: block and scalar field for each faulting plane
if output is ‘gravity’:
  1. Weights: block and scalar field
  2. Faults: block and scalar field for each faulting plane
  3. Forward gravity
if output is gradients:
  1. Lithologies: block and scalar field
  2. Faults: block and scalar field for each faulting plane
  3. Gradients of scalar field x
  4. Gradients of scalar field y
  5. Gradients of scalar field z

In addition if get_potential_at_interfaces is True, the value of the potential field at each of the interfaces is given as well

Return type:

list of _np.array

gempy_front.compute_model_at(new_grid_array, interp_data, output=’geology’, u_grade=None, get_potential_at_interfaces=False)[source]

This function does the same as compute_model() plus the addion functionallity of passing a given array of point where evaluate the model instead of using the gempy.data_management.InputData grid.

Parameters:new_grid_array (_np.array) – 2D array with XYZ (columns) coorinates
Returns:
depending on the chosen out returns different number of solutions:
if output is geology:
  1. Lithologies: block and scalar field
  2. Faults: block and scalar field for each faulting plane
if output is ‘gravity’:
  1. Weights: block and scalar field
  2. Faults: block and scalar field for each faulting plane
  3. Forward gravity
if output is gradients:
  1. Lithologies: block and scalar field
  2. Faults: block and scalar field for each faulting plane
  3. Gradients of scalar field x
  4. Gradients of scalar field y
  5. Gradients of scalar field z

In addition if get_potential_at_interfaces is True, the value of the potential field at each of the interfaces is given as well

Return type:list of _np.array
gempy_front.create_data(extent, resolution=(50, 50, 50), **kwargs)[source]

Method to create a gempy.data_management.InputData object. It is analogous to gempy.InputData()

Parameters:
  • extent (list or array) – [x_min, x_max, y_min, y_max, z_min, z_max]. Extent for the visualization of data and default of for the grid class.
  • resolution (list or array) – [nx, ny, nz]. Resolution for the visualization of data and default of for the grid class.
Keyword Arguments:
 
  • Resolution ((Optional[list])) – [nx, ny, nz]. Defaults to 50
  • path_i – Path to the data bases of interfaces. Default os.getcwd(),
  • path_f – Path to the data bases of orientations. Default os.getcwd()
Returns:

gempy.data_management.InputData

gempy_front.create_from_geomodeller_xml(fp, resolution=(50, 50, 50), return_xml=False, **kwargs)[source]

Creates InputData object from a GeoModeller xml file. Automatically extracts and sets model extent, interface and orientation data as well as the stratigraphic pile.

Parameters:
  • fp (str) – Filepath for the GeoModeller xml file to be read.
  • resolution (tuple, optional) – Tuple containing the number of voxels in dimensions (x,y,z). Defaults to 50.
  • return_xml (bool, optional) – Toggles returning the ReadGeoModellerXML instance to leverage further info from the xml file (e.g. for stratigraphic pile ordering). Defaults to True.
  • **kwargs – Keyword arguments for create_data function.
Returns:

gp.data_management.InputData

gempy_front.data_to_pickle(geo_data, path=False)[source]

Save InputData object to a python pickle (serialization of python). Be aware that if the dependencies versions used to export and import the pickle differ it may give problems

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • path (str) – path where save the pickle (without .pickle)
Returns:

None

gempy_front.export_to_vtk(geo_data, path=None, name=None, lith_block=None, vertices=None, simplices=None)[source]

Export data to a vtk file for posterior visualizations

Parameters:
  • geo_data (gempy.InputData) – All values of a DataManagement object
  • block (numpy.array) – 3D array containing the lithology block
  • path (str) – path to the location of the vtk
Returns:

None

gempy_front.get_data(geo_data, itype=’all’, numeric=False, verbosity=0)[source]

Method to return the data stored in DataFrame within a gempy.interpolator.InterpolatorData object.

Parameters:
  • geo_data (gempy.interpolator.InterpolatorData) –
  • itype (str {'all', 'interfaces', 'orientaions', 'formations', 'series', 'faults', 'fautls_relations'}) – input data type to be retrieved.
  • numeric (bool) – if True it only returns numberical properties. This may be useful due to memory issues
  • verbosity (int) – Number of properties shown
Returns:

pandas.core.frame.DataFrame

gempy_front.get_grid(geo_data)[source]

Coordinates can be found in gempy.data_management.GridClass.values

Args:
geo_data (gempy.interpolator.InterpolatorData)
Returns:
gempy.data_management.GridClass
gempy_front.get_kriging_parameters(interp_data, verbose=0)[source]

Print the kringing parameters

Parameters:
  • interp_data (gempy.data_management.InputData) –
  • verbose (int) – if > 0 print all the shape values as well.
Returns:

None

gempy_front.get_sequential_pile(geo_data)[source]

Visualize an interactive stratigraphic pile to move around the formations and the series. IMPORTANT NOTE: To have the interactive properties it is necessary the use of an interactive backend. (In notebook use: %matplotlib qt5 or notebook)

Parameters:geo_data (gempy.interpolator.InterpolatorData) –
Returns:matplotlib.pyplot.Figure
gempy_front.get_series(geo_data)[source]
Parameters:geo_data (gempy.data_management.InputData) –
Returns:Return series and formations relations
Return type:DataFrame
gempy_front.get_surfaces(interp_data, potential_lith=None, potential_fault=None, n_formation=’all’, step_size=1, original_scale=True)[source]

Compute vertices and simplices of the interfaces for its vtk visualization and further analysis

Parameters:
  • interp_data (gempy.data_management.InputData) –
  • potential_lith (ndarray) – 1D numpy array with the solution of the computation of the model containing the scalar field of potentials (second row of lith solution)
  • potential_fault (ndarray) – 1D numpy array with the solution of the computation of the model containing the scalar field of the faults (every second row of fault solution)
  • n_formation (int or 'all') – Positive integer with the number of the formation of which the surface is returned. use method get_formation_number() to get a dictionary back with the values
  • step_size (int) – resolution of the method. This is every how many voxels the marching cube method is applied
  • original_scale (bool) – choosing if the coordinates of the vertices are given in the original or the rescaled coordinates
Returns:

vertices, simpleces

gempy_front.get_th_fn(interp_data)[source]

Get the compiled theano function

Parameters:interp_data (gempy.data_management.InputData) –
Returns:
Compiled function if C or CUDA which computes the interpolation given the input data
(XYZ of dips, dip, azimuth, polarity, XYZ ref interfaces, XYZ rest interfaces)
Return type:theano.compile.function_module.Function
gempy_front.interactive_df_change_df(geo_data, only_selected=False)[source]

Confirm and return the changes made to a dataframe using qgrid interactively. To update the gempy.data_management.InputData with the modify df use the correspondant set function.

Parameters:
  • geo_data (gempy.data_management.InputData) – the same gempy.data_management.InputData object used to call interactive_df_open()
  • only_selected (bool) –
Returns:

DataFrame

gempy_front.interactive_df_open(geo_data, itype)[source]

Open the qgrid interactive DataFrame (http://qgrid.readthedocs.io/en/latest/). To seve the changes see: interactive_df_change_df()

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • itype (str {'all', 'interfaces', 'orientaions', 'formations', 'series', 'faults', 'fautls_relations'}) – input data type to be retrieved.
Returns:

Interactive DF

Return type:

DataFrame

gempy_front.plot_data(geo_data, direction=’y’, data_type=’all’, series=’all’, legend_font_size=6, **kwargs)[source]

Plot the projection of the raw data (interfaces and orientations) in 2D following a specific directions

Parameters:
  • direction (str) – xyz. Caartesian direction to be plotted
  • series (str) – series to plot
  • **kwargs – seaborn lmplot key arguments. (TODO: adding the link to them)
Returns:

None

gempy_front.plot_data_3D(geo_data, **kwargs)[source]

Plot in vtk all the input data of a model :param geo_data: Input data of the model :type geo_data: gempy.DataManagement.InputData

Returns:None
gempy_front.plot_gradient(geo_data, scalar_field, gx, gy, gz, cell_number, q_stepsize=5, direction=’y’, plot_scalar=True, **kwargs)[source]

Plot the gradient of the scalar field in a given direction.

Parameters:
  • geo_data (gempy.DataManagement.InputData) – Input data of the model
  • scalar_field (numpy.array) – scalar field to plot with the gradient
  • gx (numpy.array) – gradient in x-direction
  • gy (numpy.array) – gradient in y-direction
  • gz (numpy.array) – gradient in z-direction
  • cell_number (int) – position of the array to plot
  • q_stepsize (int) – step size between arrows to indicate gradient
  • direction (str) – xyz. Caartesian direction to be plotted
  • plot_scalar (bool) – boolean to plot scalar field
  • **kwargs – plt.contour kwargs
Returns:

None

gempy_front.plot_scalar_field(geo_data, potential_field, cell_number, N=20, direction=’y’, plot_data=True, series=’all’, *args, **kwargs)[source]

Plot a potential field in a given direction.

Parameters:
  • cell_number (int) – position of the array to plot
  • potential_field (str) – name of the potential field (or series) to plot
  • n_pf (int) – number of the potential field (or series) to plot
  • direction (str) – xyz. Caartesian direction to be plotted
  • serieDeprecated
  • **kwargs – plt.contour kwargs
Returns:

None

gempy_front.plot_section(geo_data, block, cell_number, direction=’y’, **kwargs)[source]

Plot a section of the block model

Parameters:
  • cell_number (int) – position of the array to plot
  • direction (str) – xyz. Caartesian direction to be plotted
  • interpolation (str) – Type of interpolation of plt.imshow. Default ‘none’. Acceptable values are ‘none’
  • 'bilinear', 'bicubic', (,'nearest',) –
  • 'spline36', 'hanning', 'hamming', 'hermite', 'kaiser', ('spline16',) –
  • 'catrom', 'gaussian', 'bessel', 'mitchell', 'sinc', ('quadric',) –
  • 'lanczos'**kwargs: imshow keywargs
Returns:

None

gempy_front.plot_surfaces_3D(geo_data, vertices_l, simplices_l, alpha=1, plot_data=True, size=(1920, 1080), fullscreen=False, bg_color=None)[source]

Plot in vtk the surfaces. For getting vertices and simplices See gempy.get_surfaces

Parameters:
  • vertices_l (numpy.array) – 2D array (XYZ) with the coordinates of the points
  • simplices_l (numpy.array) – 2D array with the value of the vertices that form every single triangle
  • formations_names_l (list) – Name of the formation of the surfaces
  • formation_numbers_l (list) – formation_numbers (int)
  • alpha (float) – Opacity
  • plot_data (bool) – Default True
  • size (tuple) – Resolution of the window
  • fullscreen (bool) – Launch window in full screen or not
Returns:

None

gempy_front.plot_surfaces_3D_real_time(interp_data, vertices_l, simplices_l, alpha=1, plot_data=True, posterior=None, samples=None, size=(1920, 1080), fullscreen=False)[source]

Plot in vtk the surfaces in real time. Moving the input data will affect the surfaces. IMPORTANT NOTE it is highly recommended to have the flag fast_run in the theano optimization. Also note that the time needed to compute each model increases linearly with every potential field (i.e. fault or discontinuity). It may be better to just modify each potential field individually to increase the speed (See gempy.select_series).

Parameters:
  • vertices_l (numpy.array) – 2D array (XYZ) with the coordinates of the points
  • simplices_l (numpy.array) – 2D array with the value of the vertices that form every single triangle
  • formations_names_l (list) – Name of the formation of the surfaces
  • formation_numbers_l (list) – formation_numbers (int)
  • alpha (float) – Opacity
  • plot_data (bool) – Default True
  • size (tuple) – Resolution of the window
  • fullscreen (bool) – Launch window in full screen or not
Returns:

None

gempy_front.plot_topology(geo_data, G, centroids, direction=’y’)[source]

Plot the topology adjacency graph in 2-D.

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • G (skimage.future.graph.rag.RAG) –
  • centroids (dict) – Centroid node coordinates as a dictionary with node id’s (int) as keys and (x,y,z) coordinates as values.
Keyword Args
direction (str): “x”, “y” or “z” specifying the slice direction for 2-D topology analysis. Default None.
Returns:Nothing, it just plots.
gempy_front.read_pickle(path)[source]

Read InputData object from python pickle.

Parameters:path (str) – path where save the pickle
Returns:gempy.data_management.InputData
gempy_front.rescale_data(geo_data, rescaling_factor=None)[source]

Rescale the data of a gempy.data_management.InputData object between 0 and 1 due to stability problem of the float32.

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • rescaling_factor (float) – factor of the rescaling. Default to maximum distance in one the axis
Returns:

Rescaled data

Return type:

gempy.data_management.InputData

gempy_front.rescale_factor_default(geo_data)[source]

Returns the default rescaling factor for a given geo_data

Parameters:geo_data (gempy.data_management.InputData) –
Returns:rescaling factor
Return type:float
gempy_front.select_series(geo_data, series)[source]

Return the formations of a given serie in string

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • series (list of int or list of str) – Subset of series to be selected
Returns:

New object only containing the selected series

Return type:

gempy.data_management.InputData

gempy_front.set_formations(geo_data, formations=None, formations_order=None, formations_values=None)[source]

Function to order and change the value of the model formations. The values of the formations will be the final numerical value that each formation will take in the interpolated geological model (lithology block) :param geo_data: :type geo_data: gempy.data_management.InputData :param formations_order: List with a given order of the formations. Due to the interpolation algorithm

this order is only relevant to keep consistent the colors of layers and input data. The order ultimately is determined by the geometric sedimentary order
Parameters:
  • formations (list of str) – same as formations order. you can use any
  • formations_values (list of floats or int) – values of the formations will be the final

numerical value that each formation will take in the interpolated geological model (lithology block)

Returns:formations dataframe already updated in place
Return type:DataFrame
gempy_front.set_grid(geo_data, grid)[source]

Set a new gempy.data_management.GridClass object into a gempy.data_management.InputData object.

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • grid (gempy.data_management.GridClass) –
gempy_front.set_interfaces(geo_data, interf_dataframe, append=False)[source]

Method to change or append a Dataframe to interfaces in place.

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • interf_dataframe (DataFrame) –
  • append (Bool) – if you want to append the new data frame or substitute it
gempy_front.set_interpolation_data(geo_data, **kwargs)[source]
Create a gempy.interpolator.InterpolatorData. InterpolatorData is a class that contains all the
preprocessing operations to prepare the data to compute the model. Also is the object that has to be manipulated to vary the data without recompile the modeling function.
Parameters:
  • geo_data (gempy.DataManagement.InputData) – All values of a DataManagement object
  • compile_theano (bool) – select if the theano function is compiled during the initialization. Default: True
  • compute_all (bool) –

    If true the solution gives back the block model of lithologies, the potential field and the block model of faults. If False only return the block model of lithologies. This may be important to speed

    up the computation. Default True
  • u_grade (list) – grade of the polynomial for the universal part of the Kriging interpolations. The value has to
  • either 0, 3 or 9 (be) –
  • on the number of points given as input to try to avoid singular matrix. NOTE (depends) – if during the computation
  • the model a singular matrix is returned try to reduce the u_grade of the series. (of) –
  • rescaling_factor (float) – rescaling factor of the input data to improve the stability when float32 is used. By
  • the rescaling factor is calculated to obtein values between 0 and 1. (defaut) –
Keyword Arguments:
 
  • dtype ('str') – Choosing if using float32 or float64. This is important if is intended to use the GPU
  • Also InterpolatorClass kwargs (See) –
gempy_front.geo_data

Original gempy.DataManagement.InputData object

gempy_front.geo_data_res

Rescaled data. It has the same structure has gempy.InputData

gempy_front.interpolator

Instance of the gempy.DataManagement.InterpolaorInput.InterpolatorClass. See Also gempy.DataManagement.InterpolaorInput.InterpolatorClass docs th_fn: Theano function which compute the interpolation

gempy_front.dtype

type of float

gempy_front.set_orientation_from_interfaces(geo_data, indices_array)[source]
Create and set orientations from at least 3 points of the gempy.data_management.InputData.interfaces
Dataframe
Parameters:
  • geo_data (gempy.data_management.InputData) –
  • indices_array (array-like) – 1D or 2D array with the pandas indices of the gempy.data_management.InputData.interfaces. If 2D every row of the 2D matrix will be used to create an orientation
  • verbose
Returns:

Already updated inplace

Return type:

gempy.data_management.InputData.orientations

gempy_front.set_orientations(geo_data, orient_dataframe, append=False)[source]

Method to change or append a dataframe to orientations in place. A equivalent Pandas Dataframe with [‘X’, ‘Y’, ‘Z’, ‘dip’, ‘azimuth’, ‘polarity’, ‘formation’] has to be passed.

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • interf_dataframe (DataFrame) –
  • append (Bool) – if you want to append the new data frame or substitute it
gempy_front.set_series(geo_data, series_distribution=None, order_series=None, order_formations=None, verbose=0)[source]

Function to set in place the different series of the project with their correspondent formations

Parameters:
  • geo_data (gempy.data_management.InputData) –
  • series_distribution (dict or DataFrame) – with the name of the series as key and the name of the formations as values.
  • order_series (Optional[list]) – only necessary if passed a dict (python < 3.6)order of the series by default takes the dictionary keys which until python 3.6 are random. This is important to set the erosion relations between the different series
  • order_formations (Optional[list]) – only necessary if passed a dict (python < 3.6)order of the series by default takes the dictionary keys which until python 3.6 are random. This is important to set the erosion relations between the different series
  • verbose (int) – if verbose is True plot hte sequential pile

Notes

The following dataframes will be modified in place
  1. geo_data.series: A pandas DataFrame with the series and formations relations
  2. geo_data.interfaces: one extra column with the given series
  3. geo_data.orientations: one extra column with the given series
gempy_front.topology_compute(geo_data, lith_block, fault_block, cell_number=None, direction=None, compute_areas=False, return_label_block=False)[source]

Computes model topology and returns graph, centroids and look-up-tables.

Parameters:
  • geo_data (gempy.data_management.InputData) – GemPy’s data object for the model.
  • lith_block (ndarray) – Lithology block model.
  • fault_block (ndarray) – Fault block model.
  • cell_number (int) – Cell number for 2-D slice topology analysis. Default None.
  • direction (str) – “x”, “y” or “z” specifying the slice direction for 2-D topology analysis. Default None.
  • compute_areas (bool) – If True computes adjacency areas for connected nodes in voxel number. Default False.
  • return_label_block (bool) – If True additionally returns the uniquely labeled block model as np.ndarray. Default False.
Returns:

G: Region adjacency graph object (skimage.future.graph.rag.RAG) containing the adjacency topology graph

(G.adj).

centroids (dict): Centroid node coordinates as a dictionary with node id’s (int) as keys and (x,y,z) coordinates

as values. {node id (int): tuple(x,y,z)}

labels_unique (np.array): List of all labels used. lith_to_labels_lot (dict): Dictionary look-up-table to go from lithology id to node id. labels_to_lith_lot (dict): Dictionary look-up-table to go from node id to lithology id.

Return type:

tuple

Data Management

This file is part of gempy.

gempy is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

gempy is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with gempy. If not, see <http://www.gnu.org/licenses/>.

class data_management.GridClass[source]

Class to generate grids to pass later on to a InputData class.

create_custom_grid(custom_grid)[source]

Give the coordinates of an external generated grid

Parameters:custom_grid (numpy.ndarray like) – XYZ (in columns) of the desired coordinates
Returns:Unraveled 3D numpy array where every row correspond to the xyz coordinates of a regular grid
Return type:numpy.ndarray
create_regular_grid_3d(extent, resolution)[source]

Method to create a 3D regular grid where is interpolated

Parameters:
  • extent (list) – [x_min, x_max, y_min, y_max, z_min, z_max]
  • resolution (list) – [nx, ny, nz].
Returns:

Unraveled 3D numpy array where every row correspond to the xyz coordinates of a regular grid

Return type:

numpy.ndarray

class data_management.InputData(extent, resolution=[50, 50, 50], path_i=None, path_o=None, path_f=None, **kwargs)[source]

Class that contains all raw data of our models

Parameters:
  • extent (list) – [x_min, x_max, y_min, y_max, z_min, z_max]
  • Resolution (Optional[list]) – [nx, ny, nz]. Defaults to 50
  • path_i (Optional[str]) – Path to the data bases of interfaces. Default os.getcwd(),
  • path_o (Optional[str]) – Path to the data bases of orientations. Default os.getcwd()
  • kwargs – key words passed to gempy.data_management.GridClass
interfaces

pn.core.frame.DataFrames – Pandas data frame containing the necessary information respect the interface points of the model

orientations

pn.core.frame.DataFrames – Pandas data frame containing the necessary information respect the orientations of the model

formations

pn.core.frame.DataFrames – Pandas data frame containing the formations names and the value used for each voxel in the final model and the lithological order

series

pn.core.frame.DataFrames – Pandas data frame containing the series and the formations contained on them

faults

pn.core.frame.DataFrames – Pandas data frame containing the series and if they are faults or not (otherwise they are lithologies) and in case of being fault if is finite

faults_relations

pn.core.frame.DataFrames – Pandas data frame containing the offsetting relations between each fault and the rest of the series (either other faults or lithologies)

grid

gempy.data_management.GridClass – grid object containing mainly the coordinates to interpolate the model

extent

list – [x_min, x_max, y_min, y_max, z_min, z_max]

resolution

Optional[list] – [nx, ny, nz]

add_interface(**kwargs)[source]

Adds interface to dataframe.

Parameters:**kwargs – X, Y, Z, formation, labels, order_series, series
Returns:None
add_orientation(**kwargs)[source]

Adds orientation to dataframe.

Parameters:**kwargs – G_x, G_y, G_z, X, Y, Z, azimuth, dip, formation, labels, order_series, polarity, series

Returns: Nothing

calculate_gradient()[source]

Calculate the gradient vector of module 1 given dip and azimuth to be able to plot the orientations

orientations

extra columns with components xyz of the unity vector.

calculate_orientations()[source]

Calculate and update the orientation data (azimuth and dip) from gradients in the data frame.

count_faults()[source]

Read the string names of the formations to detect automatically the number of faults.

data_to_pickle(path=False)[source]

Save InputData object to a python pickle (serialization of python). Be aware that if the dependencies versions used to export and import the pickle differ it may give problems

Parameters:path (str) – path where save the pickle
Returns:None
drop_interface(index)[source]

Drops interface from dataframe identified by index

Parameters:index – dataframe index
Returns:None
drop_orientations(index)[source]

Drops orientation from dataframe identified by index

Parameters:index – dataframe index
Returns:None
get_data(itype=’all’, numeric=False, verbosity=0)[source]

Method that returns the interfaces and orientations pandas Dataframes. Can return both at the same time or only one of the two

Parameters:
  • itype – input_data data type, either ‘orientations’, ‘interfaces’ or ‘all’ for both.
  • numeric (bool) – Return only the numerical values of the dataframe. This is much lighter database for storing traces
  • verbosity (int) – Number of properties shown
Returns:

Data frame with the raw data

Return type:

pandas.core.frame.DataFrame

get_formations()[source]
Returns:Returns a list of formations
Return type:pandas.core.frame.DataFrame
import_data_csv(path_i, path_o, **kwargs)[source]

Method to import interfaces and orientations from csv. The format is the same as the export 3D model data of GeoModeller (check in the input_data data folder for an example).

Parameters:
  • path_i (str) – path to the csv table
  • path_o (str) – path to the csv table
  • **kwargs – kwargs of :func: ~pn.read_csv
orientations

pandas.core.frame.DataFrame – Pandas data frame with the orientations data

Interfaces

pandas.core.frame.DataFrame – Pandas data frame with the interfaces data

static load_data_csv(data_type, path=’/home/miguel/PycharmProjects/gempy/docs’, **kwargs)[source]

Method to load either interface or orientations data csv files. Normally this is in which GeoModeller exports it

Parameters:
  • data_type (str) – ‘interfaces’ or ‘orientations’
  • path (str) – path to the files. Default os.getcwd()
  • **kwargs – Arbitrary keyword arguments.
Returns:

Data frame with the raw data

Return type:

pandas.core.frame.DataFrame

modify_interface(index, **kwargs)[source]

Allows modification of the x,y and/or z-coordinates of an interface at specified dataframe index.

Parameters:
  • index – dataframe index of the orientation point
  • **kwargs – X, Y, Z (int or float)
Returns:

None

modify_orientation(index, recalculate_gradient=False, recalculate_orientations=False, **kwargs)[source]

Allows modification of orientation data at specified dataframe index.

Parameters:
  • index – dataframe index of the orientation point
  • **kwargs – G_x, G_y, G_z, X, Y, Z, azimuth, dip, formation, labels, order_series, polarity
Returns:

None

order_table()[source]

First we sort the dataframes by the series age. Then we set a unique number for every formation and resort the formations. All inplace

static plane_fit(point_list)[source]

Fit plane to points in PointSet Fit an d-dimensional plane to the points in a point set. adjusted from: http://stackoverflow.com/questions/12299540/plane-fitting-to-4-or-more-xyz-points

Parameters:point_list (array_like) – array of points XYZ
Returns:Return a point, p, on the plane (the point-cloud centroid), and the normal, n.
reset_indices()[source]

Resets dataframe indices for orientations and interfaces.

Returns:None
set_annotations()[source]

Add a column in the Dataframes with latex names for each input_data paramenter.

Returns:None
set_fault_relation(rel_matrix=None)[source]

Method to set the faults that offset a given sequence and therefore also another fault

Parameters:rel_matrix (numpy.array) – 2D Boolean array with the logic. Rows affect (offset) columns
set_faults(series_name=None)[source]

Set a flag to the series that are faults.

Parameters:series_name (list or array_like) – Name of the series which are faults
set_grid(custom_grid=None, extent=None, resolution=None, grid_type=None, **kwargs)[source]

Method to initialize the class GridClass. You can pass either a custom set of points or create a regular grid

Parameters:
  • grid_type (str) – regular_3D or None
  • custom_grid (array_like) – 2D array with XYZ columns. To exploit gempy functionality the indexing has to be ij (See Also numpy.meshgrid documentation)
  • **kwargs – Arbitrary keyword arguments.
Returns:

Object that contain different grids

Return type:

self.grid(gempy.GridClass)

set_interfaces(interf_Dataframe, append=False, order_table=True)[source]

Method to change or append a Dataframe to interfaces in place. A equivalent Pandas Dataframe with [‘X’, ‘Y’, ‘Z’, ‘formation’] has to be passed.

Parameters:
  • interf_Dataframe – pandas.core.frame.DataFrame with the data
  • append – Bool: if you want to append the new data frame or substitute it
set_orientations(foliat_Dataframe, append=False, order_table=True)[source]
Method to change or append a Dataframe to orientations in place. A equivalent Pandas Dataframe with

[‘X’, ‘Y’, ‘Z’, ‘dip’, ‘azimuth’, ‘polarity’, ‘formation’] has to be passed.

Args:
interf_Dataframe: pandas.core.frame.DataFrame with the data append: Bool: if you want to append the new data frame or substitute it
set_series(series_distribution=None, order=None)[source]

Method to define the different series of the project.

Parameters:
  • series_distribution (dict) – with the name of the serie as key and the name of the formations as values.
  • order (Optional[list]) – order of the series by default takes the dictionary keys which until python 3.6 are random. This is important to set the erosion relations between the different series
Returns:

A pandas DataFrame with the series and formations relations self.interfaces: one extra column with the given series self.orientations: one extra column with the given series

Return type:

self.series

data_management.get_orientation(normal)[source]

Get orientation (dip, azimuth, polarity ) for points in all point set

Theano Graph

DEP– I need to update this string Function that generates the symbolic code to perform the interpolation. Calling this function creates

both the theano functions for the potential field and the block.
returns:theano function for the potential field theano function for the block
class theano_graph.TheanoGraph(output=’geology’, optimizer=’fast_compile’, verbose=[0], dtype=’float32’, is_fault=None, is_lith=None)[source]

This class is used to help to divide the construction of the graph into sensical parts. All its methods buildDEP2 a part of the graph. Every method can be seen as a branch and collection of branches until the last method that will be the whole tree. Every part of the graph could be compiled separately but as we increase the complexity the input of each of these methods is more and more difficult to provide (if you are in a branch close to the trunk you need all the results of the branches above)

b_vector()[source]

Creation of the independent vector b to solve the kriging system

Parameters:verbose – -deprecated-
Returns:independent vector
Return type:theano.tensor.vector
block_fault(slope=50)[source]

Compute the part of the block model of a given series (dictated by the bool array yet to be computed)

Returns:Value of lithology at every interpolated point
Return type:theano.tensor.vector
block_series(slope=50, weights=None)[source]

Compute the part of the block model of a given series (dictated by the bool array yet to be computed)

Returns:Value of lithology at every interpolated point
Return type:theano.tensor.vector
compare(a, b, slice_init, Z_x, l, n_formation, drift)[source]

Treshold of the points to interpolate given 2 potential field values. TODO: This function is the one we need to change for a sigmoid function

Parameters:
  • a (scalar) – Upper limit of the potential field
  • b (scalar) – Lower limit of the potential field
  • n_formation (scalar) – Value given to the segmentation, i.e. lithology number
  • Zx (vector) – Potential field values at all the interpolated points
Returns:

segmented values

Return type:

theano.tensor.vector

compute_a_fault(len_i_0, len_i_1, len_f_0, len_f_1, n_form_per_serie_0, n_form_per_serie_1, u_grade_iter, fault_matrix, final_block)[source]

Function that loops each fault, generating a potential field for each on them with the respective block model

Parameters:
  • len_i_0 – Lenght of rest of previous series
  • len_i_1 – Lenght of rest for the computed series
  • len_f_0 – Lenght of dips of previous series
  • len_f_1 – Length of dips of the computed series
  • n_form_per_serie_0 – Number of formations of previous series
  • n_form_per_serie_1 – Number of formations of the computed series
Returns:

block model derived from the faults that afterwards is used as a drift for the “real” data

Return type:

theano.tensor.matrix

compute_a_series(len_i_0, len_i_1, len_f_0, len_f_1, n_form_per_serie_0, n_form_per_serie_1, u_grade_iter, final_block, scalar_field_at_form)[source]

Function that loops each series, generating a potential field for each on them with the respective block model

Parameters:
  • len_i_0 – Lenght of rest of previous series
  • len_i_1 – Lenght of rest for the computed series
  • len_f_0 – Lenght of dips of previous series
  • len_f_1 – Length of dips of the computed series
  • n_form_per_serie_0 – Number of formations of previous series
  • n_form_per_serie_1 – Number of formations of the computed series
Returns:

final block model

Return type:

theano.tensor.matrix

contribution_gradient_interface(grid_val=None, weights=None)[source]

Computation of the contribution of the foliations at every point to interpolate

Returns:Contribution of all foliations (input) at every point to interpolate
Return type:theano.tensor.vector
contribution_interface(grid_val=None, weights=None)[source]

Computation of the contribution of the interfaces at every point to interpolate

Returns:Contribution of all interfaces (input) at every point to interpolate
Return type:theano.tensor.vector
contribution_interface_gradient(direction=’x’, grid_val=None, weights=None)[source]

Computation of the contribution of the foliations at every point to interpolate

Returns:Contribution of all foliations (input) at every point to interpolate
Return type:theano.tensor.vector
contribution_universal_drift(grid_val=None, weights=None, a=0, b=100000000)[source]

Computation of the contribution of the universal drift at every point to interpolate

Returns:Contribution of the universal drift (input) at every point to interpolate
Return type:theano.tensor.vector
cov_gradients(verbose=0)[source]

Create covariance function for the gradients

Returns:covariance of the gradients. Shape number of points in dip_pos x number of points in dip_pos
Return type:theano.tensor.matrix
cov_interface_gradients()[source]

Create covariance function for the gradiens :returns:

covariance of the gradients. Shape number of points in rest x number of
points in dip_pos
Return type:theano.tensor.matrix
cov_interfaces()[source]

Create covariance function for the interfaces

Returns:covariance of the interfaces. Shape number of points in rest x number of points in rest
Return type:theano.tensor.matrix
covariance_matrix()[source]

Set all the previous covariances together in the universal cokriging matrix

Returns:Multivariate covariance
Return type:theano.tensor.matrix
extend_dual_kriging()[source]

Tile the dual kriging vector to cover all the points to interpolate.So far I just make a matrix with the dimensions len(DK)x(grid) but in the future maybe I have to try to loop all this part so consume less memory

Returns:Matrix with the Dk parameters repeated for all the points to interpolate
Return type:theano.tensor.matrix
faults_contribution(weights=None, a=0, b=100000000)[source]

Computation of the contribution of the faults drift at every point to interpolate. To get these we need to compute a whole block model with the faults data

Returns:Contribution of the faults drift (input) at every point to interpolate
Return type:theano.tensor.vector
faults_matrix()[source]

This function creates the part of the graph that generates the faults function creating a “block model” at the references and the rest of the points. Then this vector has to be appended to the covariance function

Returns:
  • theano.tensor.matrix: Drift matrix for the interfaces. Shape number of points in rest x n faults. This drif is a simple addition of an arbitrary number
  • theano.tensor.matrix: Drift matrix for the gradients. Shape number of points in dips x n faults. For discrete values this matrix will be null since the derivative of a constant is 0
Return type:list
input_parameters_list()[source]

Create a list with the symbolic variables to use when we compile the theano function

Returns:
[self.dips_position_all, self.dip_angles_all, self.azimuth_all, self.polarity_all,
self.ref_layer_points_all, self.rest_layer_points_all]
Return type:list
matrices_shapes()[source]

Get all the lengths of the matrices that form the covariance matrix

Returns:length_of_CG, length_of_CGI, length_of_U_I, length_of_faults, length_of_C
scalar_field_at_all(weights=None)[source]

Compute the potential field at all the interpolation points, i.e. grid plus rest plus ref :returns: Potential fields at all points :rtype: theano.tensor.vector

solve_kriging()[source]

Solve the kriging system. This has to get substituted by a more efficient and stable method QR decomposition in all likelihood

Returns:Dual kriging parameters
Return type:theano.tensor.vector
static squared_euclidean_distances(x_1, x_2)[source]

Compute the euclidian distances in 3D between all the points in x_1 and x_2

Parameters:
  • x_1 (theano.tensor.matrix) – shape n_points x number dimension
  • x_2 (theano.tensor.matrix) – shape n_points x number dimension
Returns:

Distancse matrix. shape n_points x n_points

Return type:

theano.tensor.matrix

universal_matrix()[source]

Create the drift matrices for the potential field and its gradient

Returns:Drift matrix for the interfaces. Shape number of points in rest x 3**degree drift (except degree 0 that is 0)

theano.tensor.matrix: Drift matrix for the gradients. Shape number of points in dips x 3**degree drift (except degree 0 that is 0)

Return type:theano.tensor.matrix
x_to_interpolate(verbose=0)[source]

here I add to the grid points also the references points(to check the value of the potential field at the interfaces). Also here I will check what parts of the grid have been already computed in a previous series to avoid to recompute.

Returns:The 3D points of the given grid plus the reference and rest points
Return type:theano.tensor.matrix

Geophysics

This file is part of gempy.

gempy is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

gempy is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with gempy. If not, see <http://www.gnu.org/licenses/>.

Topology

@author: Alexander Schaaf

topology.check_adjacency(G, n1, n2)[source]

Check if n2 is adjacent/shares edge with n1.

topology.classify_edges(G, centroids, lith_block, fault_block)[source]

Classifies edges by type into stratigraphic or fault in “G.adj”. Accessible via G.adj[node1][node2][“edge_type”]

Parameters:
  • G (skimage.future.graph.rag.RAG) – Topology graph object.
  • centroids (dict) – Centroid dictionary {node id (int): tuple(x,y,z)}
  • lith_block (np.ndarray) – Shaped lithology block model.
  • fault_block (np.ndarray) – Shaped fault block model.

Returns:

topology.compare_graphs(G1, G2)[source]

Compares two graph objects using the Jaccard-Index.

Parameters:
  • G1 (skimage.future.graph.rag.RAG) – Topology graph object.
  • G2 (skimage.future.graph.rag.RAG) – Another topology graph object.
Returns:

Jaccard-Index

Return type:

(float)

topology.compute_adj_shape(n1, n2, labels_block)[source]

Compute shape of adjacency area: number of voxels in X, Y and Z (column height). (Fabian)

topology.compute_areas(G, labels_block)[source]

Computes adjacency areas and stores them in G.adj[n1][n2][“area”].

Parameters:
  • G (skimage.future.graph.rag.RAG) – Topology graph object.
  • labels_block (np.ndarray) – Uniquely labeled block model.
topology.get_centroids(label_block)[source]

Get node centroids in 2d and 3d as {node id (int): tuple(x,y,z)}.

topology.labels_lithology_lot(labels_unique, labels, block_original, verbose=0)[source]

Create LOT from label to lithology id.

topology.lithology_labels_lot(lithologies, labels, block_original, labels_unique, verbose=0)[source]

Create LOT from lithology id to label.

topology.topology_analyze(lith_block, fault_block, n_faults, areas_bool=False, return_block=False)[source]

Analyses the block models adjacency topology. Every lithological entity is described by a uniquely labeled node (centroid) and its connections to other entities by edges.

Parameters:
  • lith_block (np.ndarray) – Lithology block model
  • fault_block (np.ndarray) – Fault block model
  • n_faults (int) – Number of faults.
Keyword Arguments:
 
  • areas_bool (bool) – If True computes adjacency areas for connected nodes in voxel number. Default False.
  • return_block (bool) – If True additionally returns the uniquely labeled block model as np.ndarray.
Returns:

G: Region adjacency graph object (skimage.future.graph.rag.RAG) containing the adjacency topology graph

(G.adj).

centroids (dict): Centroid node coordinates as a dictionary with node id’s (int) as keys and (x,y,z) coordinates

as values.

labels_unique (np.array): List of all labels used. lith_to_labels_lot (dict): Dictionary look-up-table to go from lithology id to node id. labels_to_lith_lot (dict): Dictionary look-up-table to go from node id to lithology id.

Return type:

tuple

Posterior Analysis

@author: Alexander Schaaf, Miguel de la Varga

posterior_analysis.calculate_gradient(dip, az, pol)[source]

Calculates the gradient from dip, azimuth and polarity values.

posterior_analysis.calculate_ie(lbs)[source]

Computes the per-voxel and total information entropy of given block models.

posterior_analysis.calculate_information_entropy(lith_prob)[source]

Calculates information entropy for the given probability array.

posterior_analysis.calculate_information_entropy_total(ie, absolute=False)[source]

Calculate total information entropy (float) from an information entropy array.

posterior_analysis.calculate_probability_lithology(lith_blocks)[source]

Blocks must be just the lith blocks!

posterior_analysis.change_input_data(db, interp_data, i)[source]

Changes input data in interp_data to posterior input data at iteration i.

Parameters:
  • interp_data (gempy.data_management.InterpolationData) – An interp_data object with the structure we want to
  • compute.
  • i (int) – Iteration we want to recompute
Returns:

interp_data with the data of the given iteration

Return type:

gempy.data_management.InterpolationData

posterior_analysis.change_input_data_avrg(db, interp_data)[source]

Calculates average posterior parameters and updates interp_data.geo_data_res dataframes.

posterior_analysis.compute_posterior_models_all(db, interp_data, indices, u_grade=None, get_potential_at_interfaces=False)[source]

Computes block models from stored input parameters for all iterations.

Parameters:
  • () (db) – loaded pymc database (e.g. hdf5)
  • interp_data (gp.data_management.InterpolatorData) – GemPy interpolator object
  • indices (list or np.array) – Trace indices specifying which models from the database will be calculated.
  • u_grade (list, optional) –
  • get_potential_at_interfaces

Returns:

posterior_analysis.modify_plane_dip(dip, group_id, data_obj)[source]

Modify a dip angle of a plane identified by a group_id, recalculate the gradient and move the points vertically. Currently only supports the modification of dip angle - azimuth and polarity will stay the same.

Parameters:
  • dip (float) – Desired dip angle of the plane.
  • group_id (str) – Group id identifying the data points belonging to the plane.
  • ( (data_obj) – obj:): Data object to be modified (geo_data or interp_data.geo_data_res_no_basement)
Returns:

Directly modifies the given data object.

posterior_analysis.move_plane_points(normal, centroid, data_obj, interf_f)[source]

Moves interface points to fit plane of given normal and centroid in data object.

Visualization