"""
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 14
Created on 10/10 /2016
@author: Miguel de la Varga
"""
import os
from os import path
import sys
# This is for sphenix to find the packages
sys.path.append( path.dirname( path.dirname( path.abspath(__file__) ) ) )
import numpy as _np
# --DEP-- import pandas as _pn
import copy
from gempy.Visualization import PlotData2D, steano3D, vtkVisualization
from gempy.DataManagement import InputData, InterpolatorInput, GridClass
from gempy.sequential_pile import StratigraphicPile
from gempy.Topology import topology_analyze, topology_check_adjacency
import gempy.UncertaintyAnalysisPYMC2 # So far we use this type of import because the other one makes a copy and blows up some asserts
[docs]def data_to_pickle(geo_data, path=False):
"""
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
Args:
path (str): path where save the pickle
Returns:
None
"""
geo_data.data_to_pickle(path)
[docs]def read_pickle(path):
"""
Read InputData object from python pickle.
Args:
path (str): path where save the pickle
Returns:
gempy.DataManagement.InputData: Input Data object
"""
import pickle
with open(path, 'rb') as f:
# The protocol version used is detected automatically, so we do not
# have to specify it.
data = pickle.load(f)
return data
[docs]def get_series(geo_gata):
"""
Args:
geo_data (gempy.DataManagement.InputData object)
Returns:
Pandas.DataFrame: Return series and formations relations
"""
return geo_gata.series
[docs]def get_grid(geo_data):
"""
Args:
geo_data (gempy.DataManagement.InputData object)
Returns:
numpy.array: Return series and formations relations
"""
return geo_data.grid.grid
def get_resolution(geo_data):
return geo_data.resolution
def get_extent(geo_data):
return geo_data.extent
[docs]def get_data(geo_data, dtype='all', numeric=False, verbosity=0):
"""
Method that returns the interfaces and foliations pandas Dataframes. Can return both at the same time or only
one of the two
Args:
itype: input data type, either 'foliations', 'interfaces' or 'all' for both.
verbosity (int): Number of properties shown
Returns:
pandas.core.frame.DataFrame: Data frame with the raw data
"""
return geo_data.get_data(itype=dtype, numeric=numeric, verbosity=verbosity)
[docs]def create_data(extent, resolution=[50, 50, 50], **kwargs):
"""
Method to create a InputData object. It is analogous to gempy.InputData()
Args:
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 Args:
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 foliations. Default os.getcwd()
Returns:
GeMpy.DataManagement: Object that encapsulate all raw data of the project
dep: self.Plot(GeMpy_core.PlotData): Object to visualize data and results
"""
return InputData(extent, resolution, **kwargs)
[docs]def i_set_data(geo_data, dtype="foliations", action="Open"):
"""
Method to have interactive pandas tables in jupyter notebooks. The idea is to use this method to interact with
the table and i_close_set_data to recompute the parameters that depend on the changes made. I did not find a
easier solution than calling two different methods. After editing a dataframe is recommended to call
i_set_data(action= 'close') to recompute dependent variables
Args:
itype: input data type, either 'foliations' or 'interfaces'
Returns:
pandas.core.frame.DataFrame: Data frame with the changed data on real time
"""
if action == 'Close':
geo_data.i_close_set_data()
if action == 'Open':
geo_data.i_open_set_data(itype=dtype)
[docs]def select_series(geo_data, series):
"""
Return the formations of a given serie in string
:param series: list of int or list of str
:return: formations of a given serie in string separeted by |
"""
new_geo_data = copy.deepcopy(geo_data)
if type(series) == int or type(series[0]) == int:
new_geo_data.interfaces = geo_data.interfaces[geo_data.interfaces['order_series'].isin(series)]
new_geo_data.foliations = geo_data.foliations[geo_data.foliations['order_series'].isin(series)]
elif type(series[0]) == str:
new_geo_data.interfaces = geo_data.interfaces[geo_data.interfaces['series'].isin(series)]
new_geo_data.foliations = geo_data.foliations[geo_data.foliations['series'].isin(series)]
# Count faults
new_geo_data.set_faults(new_geo_data.count_faults())
# Change the dataframe with the series
new_geo_data.series = new_geo_data.series[new_geo_data.interfaces['series'].unique()]
new_geo_data.set_formation_number()
return new_geo_data
[docs]def set_series(geo_data, series_distribution=None, order_series=None, order_formations=None,
update_p_field=True, verbose=1):
"""
Method to define the different series of the project.
Args:
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:
self.series: A pandas DataFrame with the series and formations relations
self.interfaces: one extra column with the given series
self.foliations: one extra column with the given series
"""
geo_data.set_series(series_distribution=series_distribution, order=order_series)
geo_data.order_table()
if order_formations is not None:
geo_data.set_formation_number(order_formations)
# DEP
if verbose > 0:
return get_sequential_pile(geo_data)
else:
return None
def set_order_formations(geo_data, order_formations):
geo_data.set_formation_number(order_formations)
[docs]def set_interfaces(geo_data, interf_Dataframe, append=False):
"""
Method to change or append a Dataframe to interfaces in place.
Args:
interf_Dataframe: pandas.core.frame.DataFrame with the data
append: Bool: if you want to append the new data frame or substitute it
"""
geo_data.set_interfaces(interf_Dataframe, append=append)
# --DEP--
# # To update the interpolator parameters without calling a new object
# try:
# geo_data.interpolator._data = geo_data
# geo_data.interpolator._grid = geo_data.grid
#
# if update_p_field:
# geo_data.interpolator.compute_potential_fields()
# except AttributeError:
# pass
[docs]def set_foliations(geo_data, foliat_Dataframe, append=False, update_p_field=True):
"""
Method to change or append a Dataframe to foliations 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
"""
geo_data.set_foliations(foliat_Dataframe, append=append)
# To update the interpolator parameters without calling a new object
# try:
# geo_data.interpolator._data = geo_data
# geo_data.interpolator._grid = geo_data.grid
# # geo_data.interpolator._set_constant_parameteres(geo_data, geo_data.interpolator._grid)
# if update_p_field:
# geo_data.interpolator.compute_potential_fields()
# except AttributeError:
# pass
#DEP?
# def set_grid(geo_data, new_grid=None, extent=None, resolution=None, grid_type="regular_3D", **kwargs):
# """
# Method to initialize the class new_grid. So far is really simple and only has the regular new_grid type
#
# Args:
# grid_type (str): regular_3D or regular_2D (I am not even sure if regular 2D still working)
# **kwargs: Arbitrary keyword arguments.
#
# Returns:
# self.new_grid(GeMpy_core.new_grid): Object that contain different grids
# """
# if new_grid is not None:
# assert new_grid.shape[1] is 3 and len(new_grid.shape) is 2, 'The shape of new grid must be (n,3) where n is' \
# 'the number of points of the grid'
# geo_data.grid.grid = new_grid
# else:
# if not extent:
# extent = geo_data.extent
# if not resolution:
# resolution = geo_data.resolution
#
# geo_data.grid = geo_data.GridClass(extent, resolution, grid_type=grid_type, **kwargs)
[docs]def plot_data(geo_data, direction="y", data_type = 'all', series="all", legend_font_size=6, **kwargs):
"""
Plot the projection of the raw data (interfaces and foliations) in 2D following a
specific directions
Args:
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
"""
plot = PlotData2D(geo_data)
plot.plot_data(direction=direction, data_type=data_type, series=series, legend_font_size=legend_font_size, **kwargs)
# TODO saving options
[docs]def plot_section(geo_data, block, cell_number, direction="y", **kwargs):
"""
Plot a section of the block model
Args:
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'
,'nearest', 'bilinear', 'bicubic',
'spline16', 'spline36', 'hanning', 'hamming', 'hermite', 'kaiser',
'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell', 'sinc',
'lanczos'
**kwargs: imshow keywargs
Returns:
None
"""
plot = PlotData2D(geo_data)
plot.plot_block_section(cell_number, block=block, direction=direction, **kwargs)
# TODO saving options
[docs]def plot_potential_field(geo_data, potential_field, cell_number, N=20,
direction="y", plot_data=True, series="all", *args, **kwargs):
"""
Plot a potential field in a given direction.
Args:
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
serie: *Deprecated*
**kwargs: plt.contour kwargs
Returns:
None
"""
plot = PlotData2D(geo_data)
plot.plot_potential_field(potential_field, cell_number, N=N,
direction=direction, plot_data=plot_data, series=series,
*args, **kwargs)
[docs]def plot_data_3D(geo_data):
"""
Plot in vtk all the input data of a model
Args:
geo_data (gempy.DataManagement.InputData): Input data of the model
Returns:
None
"""
vv = vtkVisualization(geo_data)
vv.set_interfaces()
vv.set_foliations()
vv.render_model()
return None
[docs]def get_sequential_pile(geo_data):
"""
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 qt as interactive backend. (In notebook use:
%matplotlib qt5)
Args:
geo_data: gempy.DataManagement.InputData object
Returns:
interactive Matplotlib figure
"""
return StratigraphicPile(geo_data)
[docs]def set_interpolation_data(geo_data, **kwargs):
"""
InterpolatorInput 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.
Args:
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
be either 0, 3 or 9 (number of equations) 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.
rescaling_factor (float): rescaling factor of the input data to improve the stability when float32 is used. By
defaut the rescaling factor is calculated to obtein values between 0 and 1.
Keyword Args:
dtype ('str'): Choosing if using float32 or float64. This is important if is intended to use the GPU
See Also InterpolatorClass kwargs
Attributes:
geo_data: Original gempy.DataManagement.InputData object
geo_data_res: Rescaled data. It has the same structure has gempy.InputData
interpolator: Instance of the gempy.DataManagement.InterpolaorInput.InterpolatorClass. See Also
gempy.DataManagement.InterpolaorInput.InterpolatorClass docs
th_fn: Theano function which compute the interpolation
dtype: type of float
"""
in_data = InterpolatorInput(geo_data, **kwargs)
return in_data
[docs]def plot_surfaces_3D(geo_data, vertices_l, simplices_l,
#formations_names_l, formation_numbers_l,
alpha=1, plot_data=True,
size=(1920, 1080), fullscreen=False, bg_color=None):
"""
Plot in vtk the surfaces. For getting vertices and simplices See gempy.get_surfaces
Args:
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
"""
w = vtkVisualization(geo_data, bg_color=bg_color)
w.set_surfaces(vertices_l, simplices_l,
#formations_names_l, formation_numbers_l,
alpha)
if plot_data:
w.set_interfaces()
w.set_foliations()
w.render_model(size=size, fullscreen=fullscreen)
return w
[docs]def export_vtk_rectilinear(geo_data, block, path=None):
"""
Export data to a vtk file for posterior visualizations
Args:
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
"""
vtkVisualization.export_vtk_lith_block(geo_data, block, path)
# =====================================
# Functions for the InterpolatorData
# =====================================
# TODO check that is a interp_data object and if not try to create within the function one from the geo_data
[docs]def get_kriging_parameters(interp_data, verbose=0):
"""
Print the kringing parameters
Args:
interp_data (gempy.DataManagement.InterpolatorInput)
verbose (int): if > 0 print all the shape values as well.
Returns:
None
"""
return interp_data.interpolator.get_kriging_parameters(verbose=verbose)
[docs]def get_th_fn(interp_data, compute_all=True):
"""
Get theano function
Args:
interp_data (gempy.DataManagement.InterpolatorInput): Rescaled data.
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
Returns:
theano.function: 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 interp_data.compile_th_fn(compute_all=compute_all)
[docs]def compute_model(interp_data, output='geology', u_grade=None, get_potential_at_interfaces=False):
"""
Compute the geological model
Args:
interp_data (gempy.DataManagement.InterpolatorInput): Rescaled data.
u_grade (list): grade of the polynomial for the universal part of the Kriging interpolations. The value has to
be either 0, 3 or 9 (number of equations) 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:
numpy.array: if compute_all was chosen in gempy.DataManagement.InterpolatorInput, the first
row will be the lithology block model, the second the potential field and the third the fault
network block. if compute_all was False only the lithology block model will be computed. In
addition if get_potential_at_interfaces is True, the value of the potential field at each of
the interfaces is given as well
"""
if not getattr(interp_data, 'th_fn', None):
interp_data.compile_th_fn()
i = interp_data.get_input_data(u_grade=u_grade)
if output is 'geology':
lith_matrix, fault_matrix, potential_at_interfaces = interp_data.th_fn(*i)
# TODO check if this is necessary yet
if len(lith_matrix.shape) < 3:
_np.expand_dims(lith_matrix, 0)
_np.expand_dims(fault_matrix, 0)
interp_data.potential_at_interfaces = potential_at_interfaces
if get_potential_at_interfaces:
return lith_matrix, fault_matrix, interp_data.potential_at_interfaces
else:
return lith_matrix, fault_matrix
# TODO this should be a flag read from the compilation I guess
if output is 'gravity':
# TODO make asserts
lith_matrix, fault_matrix, potential_at_interfaces, grav = interp_data.th_fn(*i)
if len(lith_matrix.shape) < 3:
_np.expand_dims(lith_matrix, 0)
_np.expand_dims(fault_matrix, 0)
interp_data.potential_at_interfaces = potential_at_interfaces
if get_potential_at_interfaces:
return lith_matrix, fault_matrix, grav, interp_data.potential_at_interfaces
else:
return lith_matrix, fault_matrix, grav
[docs]def get_surfaces(interp_data, potential_lith=None, potential_fault=None, n_formation='all', step_size=1, original_scale=True):
"""
compute vertices and simplices of the interfaces for its vtk visualization or further
analysis
Args:
potential_lith (numpy.array): 1D numpy array with the solution of the computation of the model
containing the scalar field of potentials (second row of solution)
interp_data (gempy.DataManagement.InterpolatorInput): Interpolator object.
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
"""
try:
getattr(interp_data, 'potential_at_interfaces')
except:
raise AttributeError('You need to compute the model first')
def get_surface(potential_block, interp_data, pot_int, n_formation, step_size, original_scale):
assert n_formation > 0, 'Number of the formation has to be positive'
# In case the values are separated by series I put all in a vector
pot_int = interp_data.potential_at_interfaces.sum(axis=0)
#pot_int = np.trim_zeros(interp_data.potential_at_interfaces[n_formation]
from skimage import measure
if not potential_block.max() > pot_int[n_formation-1]:
pot_int[n_formation - 1] = potential_block.max()
print('Potential field of the surface is outside the block. Probably is due to float errors')
if not potential_block.min() < pot_int[n_formation - 1]:
pot_int[n_formation - 1] = potential_block.min()
print('Potential field of the surface is outside the block. Probably is due to float errors')
vertices, simplices, normals, values = measure.marching_cubes_lewiner(
potential_block.reshape(interp_data.geo_data_res.resolution[0],
interp_data.geo_data_res.resolution[1],
interp_data.geo_data_res.resolution[2]),
pot_int[n_formation-1],
step_size=step_size,
spacing=((interp_data.geo_data_res.extent[1] - interp_data.geo_data_res.extent[0]) / interp_data.geo_data_res.resolution[0],
(interp_data.geo_data_res.extent[3] - interp_data.geo_data_res.extent[2]) / interp_data.geo_data_res.resolution[1],
(interp_data.geo_data_res.extent[5] - interp_data.geo_data_res.extent[4]) / interp_data.geo_data_res.resolution[2]))
if original_scale:
vertices = interp_data.rescaling_factor * vertices + _np.array([interp_data._geo_data.extent[0],
interp_data._geo_data.extent[2],
interp_data._geo_data.extent[4]]).reshape(1, 3)
else:
vertices += _np.array([interp_data.geo_data_res.extent[0],
interp_data.geo_data_res.extent[2],
interp_data.geo_data_res.extent[4]]).reshape(1, 3)
return vertices, simplices
vertices = []
simplices = []
if potential_fault is not None:
assert len(_np.atleast_2d(potential_fault)) is interp_data.geo_data_res.n_faults, 'You need to pass a potential field per fault'
pot_int = interp_data.potential_at_interfaces[:interp_data.geo_data_res.n_faults + 1]
for n in interp_data.geo_data_res.interfaces['formation number'][
interp_data.geo_data_res.interfaces['isFault']].unique():
if n == 0:
continue
else:
v, s = get_surface(_np.atleast_2d(potential_fault)[n-1], interp_data, pot_int, n,
step_size=step_size, original_scale=original_scale)
vertices.append(v)
simplices.append(s)
if potential_lith is not None:
pot_int = interp_data.potential_at_interfaces[interp_data.geo_data_res.n_faults:]
# Compute the vertices of the lithologies
if n_formation == 'all':
for n in interp_data.geo_data_res.interfaces['formation number'][~interp_data.geo_data_res.interfaces['isFault']].unique():
if n == 0:
continue
else:
v, s = get_surface(potential_lith, interp_data, pot_int, n,
step_size=step_size, original_scale=original_scale)
vertices.append(v)
simplices.append(s)
else:
vertices, simplices = get_surface(potential_lith, interp_data, pot_int, n_formation,
step_size=step_size, original_scale=original_scale)
return vertices, simplices
def export_to_vtk(geo_data, path=None, lith_block=None, vertices=None, simplices=None):
if lith_block is not None:
vtkVisualization.export_vtk_lith_block(geo_data, lith_block, path=path)
if vertices is not None and simplices is not None:
vtkVisualization.export_vtk_surfaces(vertices, simplices)
[docs]def plot_surfaces_3D_real_time(interp_data, vertices_l, simplices_l,
#formations_names_l, formation_numbers_l,
alpha=1, plot_data=True, posterior=None, samples=None,
size=(1920, 1080), fullscreen=False):
"""
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).
Args:
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
"""
assert isinstance(interp_data, InterpolatorInput), 'The object has to be instance of the InterpolatorInput'
w = vtkVisualization(interp_data.geo_data_res, real_time=True)
w.set_surfaces(vertices_l, simplices_l,
#formations_names_l, formation_numbers_l,
alpha)
if posterior is not None:
assert isinstance(posterior, gempy.UncertaintyAnalysisPYMC2.Posterior), 'The object has to be instance of the Posterior class'
w.post = posterior
if samples is not None:
samp_i = samples[0]
samp_f = samples[1]
else:
samp_i = 0
samp_f = posterior.n_iter
w.create_slider_rep(samp_i, samp_f, samp_f)
w.interp_data = interp_data
if plot_data:
w.set_interfaces()
w.set_foliations()
w.render_model(size=size, fullscreen=fullscreen)
[docs]def topology_compute(geo_data, lith_block, fault_block,
cell_number=None, direction=None,
compute_areas=False, return_label_block=False):
"""
Computes model topology and returns graph, centroids and look-up-tables.
:param geo_data: geo_data object
:param lith_block: lithology block
:param fault_block: fault block
:param cell_number: (int) the slice position
:param direction: (str) "x", "y", or "z" - the slice direction
:return: (adjacency Graph object, centroid dict, labels-to-lith LOT dict, lith-to_labels LOT dict)
"""
fault_block = fault_block[::2].sum(axis=0)
if cell_number is None or direction is None: # topology of entire block
lb = lith_block.reshape(geo_data.resolution)
fb = fault_block.reshape(geo_data.resolution)
elif direction == "x":
lb = lith_block.reshape(geo_data.resolution)[cell_number, :, :]
fb = fault_block.reshape(geo_data.resolution)[cell_number, :, :]
elif direction == "y":
lb = lith_block.reshape(geo_data.resolution)[:, cell_number, :]
fb = fault_block.reshape(geo_data.resolution)[:, cell_number, :]
elif direction == "z":
lb = lith_block.reshape(geo_data.resolution)[:, :, cell_number]
fb = fault_block.reshape(geo_data.resolution)[:, :, cell_number]
return topology_analyze(lb, fb, geo_data.n_faults, areas_bool=compute_areas, return_block=return_label_block)
[docs]def topology_plot(geo_data, G, centroids, direction="y"):
"Plot topology graph."
PlotData2D.plot_topo_g(geo_data, G, centroids, direction=direction)
# =====================================
# Functions for Geophysics
# =====================================
def precomputations_gravity(interp_data, n_chunck=25, densities=None):
try:
getattr(interp_data, 'geophy')
except:
raise AttributeError('You need to set a geophysical object first. See set_geophysics_obj')
tz, select = interp_data.geophy.compute_gravity(n_chunck)
interp_data.interpolator.set_z_comp(tz, select)
if densities is not None:
set_densities(interp_data, densities)
return tz, select
def set_geophysics_obj(interp_data, ai_extent, ai_resolution, ai_z=None, range_max=None):
assert isinstance(interp_data, InterpolatorInput), 'The object has to be instance of the InterpolatorInput'
interp_data.set_geophysics_obj(ai_extent, ai_resolution, ai_z=ai_z, range_max=range_max)
return interp_data.geophy
def set_densities(interp_data, densities):
interp_data.interpolator.set_densities(densities)