Source code for Visualization

"""
    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/>.


Module with classes and methods to visualized structural geology data and potential fields of the regional modelling based on
the potential field method.
Tested on Ubuntu 14

Created on 23/09/2016

@author: Miguel de la Varga
"""

import warnings
try:
    import vtk
except ImportError:
    warnings.warn('Vtk package is not installed. No vtk visualization available.')
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import numpy as np
from os import path
import sys
# This is for sphenix to find the packages
sys.path.append( path.dirname( path.dirname( path.abspath(__file__) ) ) )
from IPython.core.debugger import Pdb
from gempy.colors import color_lot, cmap, norm
#from gempy import compute_model, get_surfaces
import gempy as gp
# TODO: inherit pygeomod classes
# import sys, os
#sns.set_context('talk')
plt.style.use(['seaborn-white', 'seaborn-talk'])

[docs]class PlotData2D(object): """ Class to make the different plot related with gempy Args: geo_data(gempy.InputData): All values of a DataManagement object block(numpy.array): 3D array containing the lithology block **kwargs: Arbitrary keyword arguments. Keyword Args: potential_field(numpy.ndarray): 3D array containing a individual potential field verbose(int): Level of verbosity during the execution of the functions (up to 5). Default 0 """ def __init__(self, geo_data, color_lot=color_lot, cmap=cmap, norm=norm, **kwargs): self._data = geo_data self._color_lot = color_lot self._cmap = cmap self._norm = norm self.formation_names = self._data.interfaces['formation'].unique() self.formation_numbers = self._data.interfaces['formation number'].unique() if 'potential_field' in kwargs: self._potential_field_p = kwargs['potential_field'] # TODO planning the whole visualization scheme. Only data, potential field # and block. 2D 3D? Improving the iteration # with pandas framework # self._set_style() # # def _set_style(self): # """ # Private function to set some plotting options # # """ # # plt.style.use(['seaborn-white', 'seaborn-talk']) # # sns.set_context("paper") # # matplotlib.rc("font", family="Helvetica")
[docs] def plot_data(self, direction="y", data_type='all', series="all", legend_font_size=8, **kwargs): """ Plot the projecton of the raw data (interfaces and foliations) in 2D following a specific directions Args: direction(str): xyz. Caartesian direction to be plotted data_type (str): type of data to plot. 'all', 'interfaces' or 'foliations' series(str): series to plot **kwargs: seaborn lmplot key arguments. (TODO: adding the link to them) Returns: Data plot """ if 'scatter_kws' not in kwargs: kwargs['scatter_kws'] = {"marker": "D", "s": 100, "edgecolors": "black", "linewidths": 1} x, y, Gx, Gy = self._slice(direction)[4:] extent = self._slice(direction)[3] aspect = (extent[1] - extent[0])/(extent[3] - extent[2]) if series == "all": series_to_plot_i = self._data.interfaces[self._data.interfaces["series"]. isin(self._data.series.columns.values)] series_to_plot_f = self._data.foliations[self._data.foliations["series"]. isin(self._data.series.columns.values)] else: series_to_plot_i = self._data.interfaces[self._data.interfaces["series"] == series] series_to_plot_f = self._data.foliations[self._data.foliations["series"] == series] # Change dictionary keys numbers for formation names for i in zip(self.formation_names, self.formation_numbers): self._color_lot[i[0]] = self._color_lot[i[1]] if data_type == 'all': p = sns.lmplot(x, y, data=series_to_plot_i, fit_reg=False, aspect=aspect, hue="formation", #scatter_kws=scatter_kws, legend=True, legend_out=True, palette=self._color_lot, **kwargs) # Plotting orientations plt.quiver(series_to_plot_f[x], series_to_plot_f[y], series_to_plot_f[Gx], series_to_plot_f[Gy], pivot="tail") if data_type == 'interfaces': p = sns.lmplot(x, y, data=series_to_plot_i, fit_reg=False, aspect=aspect, hue="formation", #scatter_kws=scatter_kws, legend=True, legend_out=True, palette=self._color_lot, **kwargs) if data_type == 'foliations': plt.quiver(series_to_plot_f[x], series_to_plot_f[y], series_to_plot_f[Gx], series_to_plot_f[Gy], pivot="tail") # # # code for moving legend outside of plot # box = p.ax.get_position() # get figure position # # reduce width of box to make room for outside legend # p.ax.set_position([box.x0, box.y0, box.width * 0.93, box.height]) # # put legend outside # p.ax.legend(loc="center right", title="Formation", # bbox_to_anchor=(1.25, 0.5), ncol=1, # prop={'size': legend_font_size}) plt.xlabel(x) plt.ylabel(y)
def _slice(self, direction, cell_number=25): """ Slice the 3D array (blocks or potential field) in the specific direction selected in the plotting functions """ _a, _b, _c = (slice(0, self._data.resolution[0]), slice(0, self._data.resolution[1]), slice(0, self._data.resolution[2])) if direction == "x": _a = cell_number x = "Y" y = "Z" Gx = "G_y" Gy = "G_z" extent_val = self._data.extent[2], self._data.extent[3], self._data.extent[4], self._data.extent[5] elif direction == "y": _b = cell_number x = "X" y = "Z" Gx = "G_x" Gy = "G_z" extent_val = self._data.extent[0], self._data.extent[1], self._data.extent[4], self._data.extent[5] elif direction == "z": _c = cell_number x = "X" y = "Y" Gx = "G_x" Gy = "G_y" extent_val = self._data.extent[0], self._data.extent[1], self._data.extent[2], self._data.extent[3] else: raise AttributeError(str(direction) + "must be a cartesian direction, i.e. xyz") return _a, _b, _c, extent_val, x, y, Gx, Gy
[docs] def plot_block_section(self, cell_number=13, block=None, direction="y", interpolation='none', plot_data=False, **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: Block plot """ if block is not None: import theano import numpy assert (type(block) is theano.tensor.sharedvar.TensorSharedVariable or type(block) is numpy.ndarray), \ 'Block has to be a theano shared object or numpy array.' if type(block) is numpy.ndarray: _block = block else: _block = block.get_value() else: try: _block = self._data.interpolator.tg.final_block.get_value() except AttributeError: raise AttributeError('There is no block to plot') plot_block = _block.reshape(self._data.resolution[0], self._data.resolution[1], self._data.resolution[2]) _a, _b, _c, extent_val, x, y = self._slice(direction, cell_number)[:-2] if plot_data: self.plot_data(direction, 'all') # TODO: plot_topo option - need fault_block for that # DEP? selecting_colors = np.unique(plot_block) if 'cmap' not in kwargs: kwargs['cmap'] = self._cmap if 'norm' not in kwargs: kwargs['norm'] = self._norm im = plt.imshow(plot_block[_a, _b, _c].T, origin="bottom",# cmap=self._cmap, norm=self._norm, extent=extent_val, interpolation=interpolation, **kwargs) import matplotlib.patches as mpatches colors = [im.cmap(im.norm(value)) for value in self.formation_numbers] patches = [ mpatches.Patch(color=colors[i], label=self.formation_names[i]) for i in range(len(self.formation_names))] plt.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) plt.xlabel(x) plt.ylabel(y)
[docs] def plot_potential_field(self, 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: Potential field plot """ if plot_data: self.plot_data(direction, 'all') _a, _b, _c, extent_val, x, y = self._slice(direction, cell_number)[:-2] plt.contour(potential_field.reshape( self._data.resolution[0], self._data.resolution[1], self._data.resolution[2])[_a, _b, _c].T, N, extent=extent_val, *args, **kwargs) if 'colorbar' in kwargs: plt.colorbar() # plt.title(self._data.interfaces['series'].unique()[n_pf]) plt.xlabel(x) plt.ylabel(y)
def plot_topo_g(geo_data, G, centroids, direction="y"): if direction == "y": c1, c2 = (0, 2) e1 = geo_data.extent[1] - geo_data.extent[0] e2 = geo_data.extent[5] - geo_data.extent[4] d1 = geo_data.extent[0] d2 = geo_data.extent[4] if len(list(centroids.items())[0][1]) == 2: c1, c2 = (0, 1) r1 = geo_data.resolution[0] r2 = geo_data.resolution[2] elif direction == "x": c1, c2 = (1, 2) e1 = geo_data.extent[3] - geo_data.extent[2] e2 = geo_data.extent[5] - geo_data.extent[4] d1 = geo_data.extent[2] d2 = geo_data.extent[4] if len(list(centroids.items())[0][1]) == 2: c1, c2 = (0, 1) r1 = geo_data.resolution[1] r2 = geo_data.resolution[2] elif direction == "z": c1, c2 = (0, 1) e1 = geo_data.extent[1] - geo_data.extent[0] e2 = geo_data.extent[3] - geo_data.extent[2] d1 = geo_data.extent[0] d2 = geo_data.extent[2] if len(list(centroids.items())[0][1]) == 2: c1, c2 = (0, 1) r1 = geo_data.resolution[0] r2 = geo_data.resolution[1] for edge in G.edges_iter(): a, b = edge # if G.adj[a][b]["edge_type"] == "stratigraphic": # plt.plot(np.array([centroids[a][c1], centroids[b][c1]]) * e1 / r1 + d1, # np.array([centroids[a][c2], centroids[b][c2]]) * e2 / r2 + d2, "gray", linewidth=2) # elif G.adj[a][b]["edge_type"] == "fault": # plt.plot(np.array([centroids[a][c1], centroids[b][c1]]) * e1 / r1 + d1, # np.array([centroids[a][c2], centroids[b][c2]]) * e2 / r2 + d2, "black", linewidth=2) plt.plot(np.array([centroids[a][c1], centroids[b][c1]]) * e1 / r1 + d1, np.array([centroids[a][c2], centroids[b][c2]]) * e2 / r2 + d2, "black", linewidth=0.75) for node in G.nodes_iter(): plt.plot(centroids[node][c1] * e1 / r1 + d1, centroids[node][c2] * e2 / r2 +d2, marker="o", color="black", markersize=10, alpha=0.75) plt.text(centroids[node][c1] * e1 / r1 + d1, centroids[node][c2] * e2 / r2 + d2, str(node), color="white", size=6, ha="center", va="center", weight="ultralight", family="monospace")
[docs] @staticmethod def annotate_plot(frame, label_col, x, y, **kwargs): """ Annotate the plot of a given DataFrame using one of its columns Should be called right after a DataFrame or series plot method, before telling matplotlib to show the plot. Parameters ---------- frame : pandas.DataFrame plot_col : str The string identifying the column of frame that was plotted label_col : str The string identifying the column of frame to be used as label kwargs: Other key-word args that should be passed to plt.annotate Returns ------- None Notes ----- After calling this function you should call plt.show() to get the results. This function only adds the annotations, it doesn't show them. """ import matplotlib.pyplot as plt # Make sure we have pyplot as plt for label, x, y in zip(frame[label_col], frame[x], frame[y]): plt.annotate(label, xy=(x + 0.2, y + 0.15), **kwargs)
class steano3D(): def __init__(self, geo_data): self._data = geo_data def plot3D_steno(self, block, project, plot=True, **kwargs): import steno3d import numpy as np steno3d.login() description = kwargs.get('description', 'Nothing') proj = steno3d.Project( title=project, description=description, public=True, ) mesh = steno3d.Mesh3DGrid(h1=np.ones(self._data.resolution[0]) * (self._data.extent[0] - self._data.extent[1]) / (self._data.resolution[0] - 1), h2=np.ones(self._data.resolution[1]) * (self._data.extent[2] - self._data.extent[3]) / (self._data.resolution[1] - 1), h3=np.ones(self._data.resolution[2]) * (self._data.extent[4] - self._data.extent[5]) / (self._data.resolution[2] - 1)) data = steno3d.DataArray( title='Lithologies', array=block) vol = steno3d.Volume(project=proj, mesh=mesh, data=[dict(location='CC', data=data)]) vol.upload() if plot: return vol.plot()
[docs]class vtkVisualization: """ Class to visualize data and results in 3D. Init will create all the render properties while the method render model will lunch the window. Using set_interfaces, set_foliations and set_surfaces in between can be chosen what will be displayed. Args: geo_data(gempy.InputData): All values of a DataManagement object ren_name (str): Name of the renderer window verbose (int): Verbosity for certain functions Attributes: renWin(vtk.vtkRenderWindow()) camera_list (list): list of cameras for the distinct renderers ren_list (list): list containing the vtk renderers """ def __init__(self, geo_data, ren_name='GemPy 3D-Editor', verbose=0, color_lot=color_lot, real_time=False, bg_color=None): self.real_time = real_time # self.C_LOT = self.color_lot_create(geo_data) self.geo_data = geo_data self.interp_data = None self.C_LOT = color_lot # Number of renders self.n_ren = 4 self.formation_number = geo_data.interfaces['formation number'].unique() self.formation_name = geo_data.interfaces['formation'].unique() # Extents self.extent = geo_data.extent _e = geo_data.extent self._e_dx = _e[1] - _e[0] self._e_dy = _e[3] - _e[2] self._e_dz = _e[5] - _e[4] self._e_d_avrg = (self._e_dx + self._e_dy + self._e_dz) / 3 # Resolution self.res = geo_data.resolution # create render window, settings self.renwin = vtk.vtkRenderWindow() self.renwin.SetWindowName(ren_name) # Set 4 renderers. ie 3D, X,Y,Z projections self.ren_list = self.create_ren_list() # create interactor and set interactor style, assign render window self.interactor = vtk.vtkRenderWindowInteractor() self.interactor.SetRenderWindow(self.renwin) # 3d model camera for the 4 renders self.camera_list = self._create_cameras(self.extent, verbose=verbose) # Setting the camera and the background color to the renders self.set_camera_backcolor(color=bg_color) # Creating the axis for e, r in enumerate(self.ren_list): # add axes actor to all renderers axe = self._create_axes(self.camera_list[e]) r.AddActor(axe) r.ResetCamera()
[docs] def render_model(self, size=(1920, 1080), fullscreen=False): """ Method to launch the window Args: size (tuple): Resolution of the window fullscreen (bool): Launch window in full screen or not Returns: """ # initialize and start the app if fullscreen: self.renwin.FullScreenOn() self.renwin.SetSize(size) self.interactor.Initialize() self.interactor.Start() # close_window(interactor) del self.renwin, self.interactor
[docs] def create_surface_points(self, vertices): """ Method to create the points that form the surfaces Args: vertices (numpy.array): 2D array (XYZ) with the coordinates of the points Returns: vtk.vtkPoints: with the coordinates of the points """ Points = vtk.vtkPoints() for v in vertices: Points.InsertNextPoint(v) return Points
[docs] def create_surface_triangles(self, simplices): """ Method to create the Triangles that form the surfaces Args: simplices (numpy.array): 2D array with the value of the vertices that form every single triangle Returns: vtk.vtkTriangle """ Triangles = vtk.vtkCellArray() Triangle = vtk.vtkTriangle() for s in simplices: Triangle.GetPointIds().SetId(0, s[0]) Triangle.GetPointIds().SetId(1, s[1]) Triangle.GetPointIds().SetId(2, s[2]) Triangles.InsertNextCell(Triangle) return Triangles
[docs] def create_surface(self, vertices, simplices, fn, alpha=1): """ Method to create the polydata that define the surfaces Args: vertices (numpy.array): 2D array (XYZ) with the coordinates of the points simplices (numpy.array): 2D array with the value of the vertices that form every single triangle fn (int): Formation number alpha (float): Opacity Returns: vtk.vtkActor, vtk.vtkPolyDataMapper, vtk.vtkPolyData """ surf_polydata = vtk.vtkPolyData() surf_polydata.SetPoints(self.create_surface_points(vertices)) surf_polydata.SetPolys(self.create_surface_triangles(simplices)) surf_polydata.Modified() surf_mapper = vtk.vtkPolyDataMapper() surf_mapper.SetInputData(surf_polydata) surf_mapper.Update() surf_actor = vtk.vtkActor() surf_actor.SetMapper(surf_mapper) surf_actor.GetProperty().SetColor(self.C_LOT[fn]) surf_actor.GetProperty().SetOpacity(alpha) return surf_actor, surf_mapper, surf_polydata
[docs] def create_sphere(self, X, Y, Z, fn, n_sphere=0, n_render=0, n_index=0, r=0.03): """ Method to create the sphere that represent the interfaces points Args: X: X coord Y: Y coord Z: Z corrd fn (int): Formation number n_sphere (int): Number of the sphere n_render (int): Number of the render where the sphere belongs n_index (int): index value in the PandasDataframe of InupData.interfaces r (float): radio of the sphere Returns: vtk.vtkSphereWidget """ s = vtk.vtkSphereWidget() s.SetInteractor(self.interactor) s.SetRepresentationToSurface() s.r_f = self._e_d_avrg * r s.PlaceWidget(X - s.r_f, X + s.r_f, Y - s.r_f, Y + s.r_f, Z - s.r_f, Z + s.r_f) s.GetSphereProperty().SetColor(self.C_LOT[fn]) s.SetCurrentRenderer(self.ren_list[n_render]) s.n_sphere = n_sphere s.n_render = n_render s.index = n_index s.AddObserver("EndInteractionEvent", self.sphereCallback) # EndInteractionEvent s.On() return s
[docs] def create_foliation(self, X, Y, Z, fn, Gx, Gy, Gz, n_plane=0, n_render=0, n_index=0, alpha=0.5): """ Method to create a plane given a foliation Args: X : X coord Y: Y coord Z: Z corrd fn (int): Formation number Gx (str): Component of the gradient x Gy (str): Component of the gradient y Gz (str): Component of the gradient z n_plane (int): Number of the plane n_render (int): Number of the render where the plane belongs n_index (int): index value in the PandasDataframe of InupData.interfaces alpha: Opacity of the plane Returns: vtk.vtkPlaneWidget """ d = vtk.vtkPlaneWidget() d.SetInteractor(self.interactor) d.SetRepresentationToSurface() # Position source = vtk.vtkPlaneSource() source.SetCenter(X, Y, Z) source.SetNormal(Gx, Gy, Gz) source.Update() d.SetInputData(source.GetOutput()) d.SetHandleSize(0.05) min_extent = np.min([self._e_dx, self._e_dy, self._e_dz]) d.SetPlaceFactor(min_extent/10) d.PlaceWidget() d.SetNormal(Gx, Gy, Gz) d.GetPlaneProperty().SetColor(self.C_LOT[fn]) d.GetHandleProperty().SetColor(self.C_LOT[fn]) d.GetHandleProperty().SetOpacity(alpha) d.SetCurrentRenderer(self.ren_list[n_render]) d.n_plane = n_plane d.n_render = n_render d.index = n_index d.AddObserver("EndInteractionEvent", self.planesCallback) d.On() return d
[docs] def set_surfaces(self, vertices, simplices, #formations, fns, alpha=1): """ Create all the surfaces and set them to the corresponding renders for their posterior visualization with render_model Args: vertices (list): list of 3D numpy arrays containing the points that form each plane simplices (list): list of 3D numpy arrays containing the verticies that form every triangle formations (list): ordered list of strings containing the name of the formations to represent fns (list): ordered list of formation numbers (int) alpha: Opacity of the plane Returns: None """ self.surf_rend_1 = [] # self.s_rend_2 = [] # self.s_rend_3 = [] # self.s_rend_4 = [] # self.s_mapper = [] # self.s_polydata = [] formations = self.formation_name fns = self.formation_number assert type( vertices) is list, 'vertices and simpleces have to be a list of arrays even when only one formation' \ 'is passed' assert 'DefaultBasement' not in formations, 'Remove DefaultBasement from the list of formations' for v, s, fn in zip(vertices, simplices, fns): act, map, pol = self.create_surface(v, s, fn, alpha) self.surf_rend_1.append(act) # self.s_mapper.append(map) # self.s_polydata.append(pol) #print(self.s_rend_1) self.ren_list[0].AddActor(act) self.ren_list[1].AddActor(act) self.ren_list[2].AddActor(act) self.ren_list[3].AddActor(act)
[docs] def set_interfaces(self): """ Create all the interfaces points and set them to the corresponding renders for their posterior visualization with render_model Returns: None """ self.s_rend_1 = [] self.s_rend_2 = [] self.s_rend_3 = [] self.s_rend_4 = [] for e, val in enumerate(self.geo_data.interfaces.iterrows()): index = val[0] row = val[1] self.s_rend_1.append(self.create_sphere(row['X'], row['Y'], row['Z'], row['formation number'], n_sphere=e, n_render=0, n_index=index)) self.s_rend_2.append(self.create_sphere(row['X'], row['Y'], row['Z'], row['formation number'], n_sphere=e, n_render=1, n_index=index)) self.s_rend_3.append(self.create_sphere(row['X'], row['Y'], row['Z'], row['formation number'], n_sphere=e, n_render=2, n_index=index)) self.s_rend_4.append(self.create_sphere(row['X'], row['Y'], row['Z'], row['formation number'], n_sphere=e, n_render=3, n_index=index))
[docs] def set_foliations(self): """ Create all the foliations and set them to the corresponding renders for their posterior visualization with render_model Returns: None """ self.f_rend_1 = [] self.f_rend_2 = [] self.f_rend_3 = [] self.f_rend_4 = [] for e, val in enumerate(self.geo_data.foliations.iterrows()): index = val[0] row = val[1] self.f_rend_1.append(self.create_foliation(row['X'], row['Y'], row['Z'], row['formation number'], row['G_x'], row['G_y'], row['G_z'], n_plane=e, n_render=0, n_index=index)) self.f_rend_2.append(self.create_foliation(row['X'], row['Y'], row['Z'], row['formation number'], row['G_x'], row['G_y'], row['G_z'], n_plane=e, n_render=1, n_index=index)) self.f_rend_3.append(self.create_foliation(row['X'], row['Y'], row['Z'], row['formation number'], row['G_x'], row['G_y'], row['G_z'], n_plane=e, n_render=2, n_index=index)) self.f_rend_4.append(self.create_foliation(row['X'], row['Y'], row['Z'], row['formation number'], row['G_x'], row['G_y'], row['G_z'], n_plane=e, n_render=3, n_index=index))
def create_slider_rep(self, min=0, max=10, val=0): slider_rep = vtk.vtkSliderRepresentation2D() slider_rep.SetMinimumValue(min) slider_rep.SetMaximumValue(max) slider_rep.SetValue(val) slider_rep.GetPoint1Coordinate().SetCoordinateSystemToDisplay() slider_rep.GetPoint1Coordinate().SetValue(0, 40) slider_rep.GetPoint2Coordinate().SetCoordinateSystemToDisplay() slider_rep.GetPoint2Coordinate().SetValue(300, 40) self.slider_w = vtk.vtkSliderWidget() self.slider_w.SetInteractor(self.interactor) self.slider_w.SetRepresentation(slider_rep) self.slider_w.SetCurrentRenderer(self.ren_list[0]) self.slider_w.SetAnimationModeToAnimate() self.slider_w.On() self.slider_w.AddObserver('EndInteractionEvent', self.slideCallback) def slideCallback(self, obj, event): self.post.change_input_data(self.interp_data, obj.GetRepresentation().GetValue()) # lith_block, fault_block = gp.compute_model(self.interp_data) try: for surf in self.surf_rend_1: self.ren_list[0].RemoveActor(surf) self.ren_list[1].RemoveActor(surf) self.ren_list[2].RemoveActor(surf) self.ren_list[3].RemoveActor(surf) ver, sim = self.update_surfaces_real_time(self.interp_data) self.set_surfaces(ver, sim) except AttributeError: print('no surf') pass try: for sph in zip(self.s_rend_1, self.s_rend_2, self.s_rend_3, self.s_rend_4, self.geo_data.interfaces.iterrows()): row_i = sph[4][1] sph[0].PlaceWidget(row_i['X'] - sph[0].r_f, row_i['X'] + sph[0].r_f, row_i['Y'] - sph[0].r_f, row_i['Y'] + sph[0].r_f, row_i['Z'] - sph[0].r_f, row_i['Z'] + sph[0].r_f) sph[1].PlaceWidget(row_i['X'] - sph[1].r_f, row_i['X'] + sph[1].r_f, row_i['Y'] - sph[1].r_f, row_i['Y'] + sph[1].r_f, row_i['Z'] - sph[1].r_f, row_i['Z'] + sph[1].r_f) sph[2].PlaceWidget(row_i['X'] - sph[2].r_f, row_i['X'] + sph[2].r_f, row_i['Y'] - sph[2].r_f, row_i['Y'] + sph[2].r_f, row_i['Z'] - sph[2].r_f, row_i['Z'] + sph[2].r_f) sph[3].PlaceWidget(row_i['X'] - sph[3].r_f, row_i['X'] + sph[3].r_f, row_i['Y'] - sph[3].r_f, row_i['Y'] + sph[3].r_f, row_i['Z'] - sph[3].r_f, row_i['Z'] + sph[3].r_f) except AttributeError: pass try: for fol in (zip(self.f_rend_1, self.f_rend_2, self.f_rend_3, self.f_rend_4, self.geo_data.foliations.iterrows())): row_f = fol[4][1] fol[0].SetCenter(row_f['X'], row_f['Y'], row_f['Z']) fol[0].SetNormal(row_f['G_x'], row_f['G_y'], row_f['G_z']) # self.ren_list[0].RemoveActor(fol[0]) # self.ren_list[1].RemoveActor(fol[1]) # self.ren_list[2].RemoveActor(fol[2]) # self.ren_list[3].RemoveActor(fol[3]) # self.set_foliations() except AttributeError: pass
[docs] def sphereCallback(self, obj, event): """ Function that rules what happens when we move a sphere. At the moment we update the other 3 renderers and update the pandas data frame """ # Resetting the xy camera when a sphere is moving to be able to change only 2D fp = self.ren_list[1].GetActiveCamera().GetFocalPoint() p = self.ren_list[1].GetActiveCamera().GetPosition() dist = np.sqrt((p[0] - fp[0]) ** 2 + (p[1] - fp[1]) ** 2 + (p[2] - fp[2]) ** 2) self.ren_list[1].GetActiveCamera().SetPosition(fp[0], fp[1], fp[2] + dist) self.ren_list[1].GetActiveCamera().SetViewUp(0.0, 1.0, 0.0) # Resetting the yz camera when a sphere is moving to be able to change only 2D fp = self.ren_list[2].GetActiveCamera().GetFocalPoint() p = self.ren_list[2].GetActiveCamera().GetPosition() dist = np.sqrt((p[0] - fp[0]) ** 2 + (p[1] - fp[1]) ** 2 + (p[2] - fp[2]) ** 2) self.ren_list[2].GetActiveCamera().SetPosition(fp[0] + dist, fp[1], fp[2]) self.ren_list[2].GetActiveCamera().SetViewUp(0.0, -1.0, 0.0) self.ren_list[2].GetActiveCamera().Roll(90) # Resetting the xz camera when a sphere is moving to be able to change only 2D fp = self.ren_list[3].GetActiveCamera().GetFocalPoint() p = self.ren_list[3].GetActiveCamera().GetPosition() dist = np.sqrt((p[0] - fp[0]) ** 2 + (p[1] - fp[1]) ** 2 + (p[2] - fp[2]) ** 2) self.ren_list[3].GetActiveCamera().SetPosition(fp[0], fp[1] - dist, fp[2]) self.ren_list[3].GetActiveCamera().SetViewUp(-1.0, 0.0, 0.0) self.ren_list[3].GetActiveCamera().Roll(90) # Get new position of the sphere new_center = obj.GetCenter() # Get which sphere we are moving index = obj.index # Modify Pandas DataFrame self.geo_data.interface_modify(index, X=new_center[0], Y=new_center[1], Z=new_center[2]) # Check what render we are working with render = obj.n_render r_f = obj.r_f # Update the other renderers if render != 0: self.s_rend_1[obj.n_sphere].PlaceWidget(new_center[0] - r_f, new_center[0] + r_f, new_center[1] - r_f, new_center[1] + r_f, new_center[2] - r_f, new_center[2] + r_f) if render != 1: self.s_rend_2[obj.n_sphere].PlaceWidget(new_center[0] - r_f, new_center[0] + r_f, new_center[1] - r_f, new_center[1] + r_f, new_center[2] - r_f, new_center[2] + r_f) if render != 2: self.s_rend_3[obj.n_sphere].PlaceWidget(new_center[0] - r_f, new_center[0] + r_f, new_center[1] - r_f, new_center[1] + r_f, new_center[2] - r_f, new_center[2] + r_f) if render != 3: self.s_rend_4[obj.n_sphere].PlaceWidget(new_center[0] - r_f, new_center[0] + r_f, new_center[1] - r_f, new_center[1] + r_f, new_center[2] - r_f, new_center[2] + r_f) if self.real_time: for surf in self.surf_rend_1: self.ren_list[0].RemoveActor(surf) self.ren_list[1].RemoveActor(surf) self.ren_list[2].RemoveActor(surf) self.ren_list[3].RemoveActor(surf) vertices, simpleces = self.update_surfaces_real_time(self.interp_data) # print(vertices[0][60]) self.set_surfaces(vertices, simpleces)
# self.renwin.Render()
[docs] def planesCallback(self, obj, event): """ Function that rules what happens when we move a plane. At the moment we update the other 3 renderers and update the pandas data frame """ # Resetting the xy camera when a sphere is moving to be able to change only 2D fp = self.ren_list[1].GetActiveCamera().GetFocalPoint() p = self.ren_list[1].GetActiveCamera().GetPosition() dist = np.sqrt((p[0] - fp[0]) ** 2 + (p[1] - fp[1]) ** 2 + (p[2] - fp[2]) ** 2) self.ren_list[1].GetActiveCamera().SetPosition(fp[0], fp[1], fp[2] + dist) self.ren_list[1].GetActiveCamera().SetViewUp(0.0, 1.0, 0.0) # Resetting the yz camera when a sphere is moving to be able to change only 2D fp = self.ren_list[2].GetActiveCamera().GetFocalPoint() p = self.ren_list[2].GetActiveCamera().GetPosition() dist = np.sqrt((p[0] - fp[0]) ** 2 + (p[1] - fp[1]) ** 2 + (p[2] - fp[2]) ** 2) self.ren_list[2].GetActiveCamera().SetPosition(fp[0] + dist, fp[1], fp[2]) self.ren_list[2].GetActiveCamera().SetViewUp(0.0, -1.0, 0.0) self.ren_list[2].GetActiveCamera().Roll(90) # Resetting the xz camera when a sphere is moving to be able to change only 2D fp = self.ren_list[3].GetActiveCamera().GetFocalPoint() p = self.ren_list[3].GetActiveCamera().GetPosition() dist = np.sqrt((p[0] - fp[0]) ** 2 + (p[1] - fp[1]) ** 2 + (p[2] - fp[2]) ** 2) self.ren_list[3].GetActiveCamera().SetPosition(fp[0], fp[1] - dist, fp[2]) self.ren_list[3].GetActiveCamera().SetViewUp(-1.0, 0.0, 0.0) self.ren_list[3].GetActiveCamera().Roll(90) # Get new position of the plane and GxGyGz new_center = obj.GetCenter() new_normal = obj.GetNormal() # Get which plane we are moving index = obj.index new_source = vtk.vtkPlaneSource() new_source.SetCenter(new_center) new_source.SetNormal(new_normal) new_source.Update() # Modify Pandas DataFrame # update the gradient vector components and its location self.geo_data.foliation_modify(index, X=new_center[0], Y=new_center[1], Z=new_center[2], G_x=new_normal[0], G_y=new_normal[1], G_z=new_normal[2]) # update the dip and azimuth values according to the new gradient self.geo_data.calculate_orientations() # Update the other renderers render = obj.n_render # TODO: get the plane source, change the stuff there, and then update the widget accordingly? if render != 0: self.f_rend_1[obj.n_plane].SetInputData(new_source.GetOutput()) self.f_rend_1[obj.n_plane].SetNormal(new_normal) self.f_rend_1[obj.n_plane].SetCenter(new_center[0], new_center[1], new_center[2]) if render != 1: self.f_rend_2[obj.n_plane].SetInputData(new_source.GetOutput()) self.f_rend_2[obj.n_plane].SetNormal(new_normal) self.f_rend_2[obj.n_plane].SetCenter(new_center[0], new_center[1], new_center[2]) if render != 2: self.f_rend_3[obj.n_plane].SetInputData(new_source.GetOutput()) self.f_rend_3[obj.n_plane].SetNormal(new_normal) self.f_rend_3[obj.n_plane].SetCenter(new_center[0], new_center[1], new_center[2]) if render != 3: self.f_rend_4[obj.n_plane].SetInputData(new_source.GetOutput()) self.f_rend_4[obj.n_plane].SetNormal(new_normal) self.f_rend_4[obj.n_plane].SetCenter(new_center[0], new_center[1], new_center[2]) if self.real_time: for surf in self.surf_rend_1: self.ren_list[0].RemoveActor(surf) self.ren_list[1].RemoveActor(surf) self.ren_list[2].RemoveActor(surf) self.ren_list[3].RemoveActor(surf) vertices, simpleces = self.update_surfaces_real_time(self.interp_data) # print(vertices[0][60]) self.set_surfaces(vertices, simpleces)
[docs] def create_ren_list(self): """ Create a list of the 4 renderers we use. One general view and 3 cartersian projections Returns: list: list of renderers """ # viewport dimensions setup xmins = [0, 0.6, 0.6, 0.6] xmaxs = [0.6, 1, 1, 1] ymins = [0, 0, 0.33, 0.66] ymaxs = [1, 0.33, 0.66, 1] # create list of renderers, set vieport values ren_list = [] for i in range(self.n_ren): # append each renderer to list of renderers ren_list.append(vtk.vtkRenderer()) # add each renderer to window self.renwin.AddRenderer(ren_list[-1]) # set viewport for each renderer ren_list[-1].SetViewport(xmins[i], ymins[i], xmaxs[i], ymaxs[i]) return ren_list
def _create_cameras(self, extent, verbose=0): """ Create the 4 cameras for each renderer """ _e = extent _e_dx = _e[1] - _e[0] _e_dy = _e[3] - _e[2] _e_dz = _e[5] - _e[4] _e_d_avrg = (_e_dx + _e_dy + _e_dz) / 3 _e_max = np.argmax(_e) # General camera model_cam = vtk.vtkCamera() model_cam.SetPosition(_e[_e_max] * 5, _e[_e_max] * 5, _e[_e_max] * 5) model_cam.SetFocalPoint(np.min(_e[0:2]) + _e_dx / 2, np.min(_e[2:4]) + _e_dy / 2, np.min(_e[4:]) + _e_dz / 2) model_cam.SetViewUp(-0.239, 0.155, 0.958) # XY camera RED xy_cam = vtk.vtkCamera() xy_cam.SetPosition(np.min(_e[0:2]) + _e_dx / 2, np.min(_e[2:4]) + _e_dy / 2, _e[_e_max] * 4) xy_cam.SetFocalPoint(np.min(_e[0:2]) + _e_dx / 2, np.min(_e[2:4]) + _e_dy / 2, np.min(_e[4:]) + _e_dz / 2) # YZ camera GREEN yz_cam = vtk.vtkCamera() yz_cam.SetPosition(_e[_e_max] * 4, np.min(_e[2:4]) + _e_dy / 2, np.min(_e[4:]) + _e_dz / 2) yz_cam.SetFocalPoint(np.min(_e[0:2]) + _e_dx / 2, np.min(_e[2:4]) + _e_dy / 2, np.min(_e[4:]) + _e_dz / 2) yz_cam.SetViewUp(0, -1, 0) yz_cam.Roll(90) # XZ camera BLUE xz_cam = vtk.vtkCamera() xz_cam.SetPosition(np.min(_e[0:2]) + _e_dx / 2, - _e[_e_max] * 4, np.min(_e[4:]) + _e_dz / 2) xz_cam.SetFocalPoint(np.min(_e[0:2]) + _e_dx / 2, np.min(_e[2:4]) + _e_dy / 2, np.min(_e[4:]) + _e_dz / 2) xz_cam.SetViewUp(0, 1, 0) xz_cam.Roll(0) # camera position debugging if verbose == 1: print("RED XY:", xy_cam.GetPosition()) print("RED FP:", xy_cam.GetFocalPoint()) print("GREEN YZ:", yz_cam.GetPosition()) print("GREEN FP:", yz_cam.GetFocalPoint()) print("BLUE XZ:", xz_cam.GetPosition()) print("BLUE FP:", xz_cam.GetFocalPoint()) return [model_cam, xy_cam, yz_cam, xz_cam]
[docs] def set_camera_backcolor(self, color=None): """ define background colors of the renderers """ if color is None: color = (66 / 250, 66 / 250, 66 / 250) ren_color = [color for i in range(self.n_ren)] for i in range(self.n_ren): # set active camera for each renderer self.ren_list[i].SetActiveCamera(self.camera_list[i]) # set background color for each renderer self.ren_list[i].SetBackground(ren_color[i][0], ren_color[i][1], ren_color[i][2])
def _create_axes(self, camera, verbose=0, tick_vis=True): """ Create the axes boxes """ cube_axes_actor = vtk.vtkCubeAxesActor() cube_axes_actor.SetBounds(self.geo_data.extent) cube_axes_actor.SetCamera(camera) if verbose == 1: print(cube_axes_actor.GetAxisOrigin()) # set axes and label colors cube_axes_actor.GetTitleTextProperty(0).SetColor(1.0, 0.0, 0.0) cube_axes_actor.GetLabelTextProperty(0).SetColor(1.0, 0.0, 0.0) cube_axes_actor.GetTitleTextProperty(1).SetColor(0.0, 1.0, 0.0) cube_axes_actor.GetLabelTextProperty(1).SetColor(0.0, 1.0, 0.0) cube_axes_actor.GetTitleTextProperty(2).SetColor(0.0, 0.0, 1.0) cube_axes_actor.GetLabelTextProperty(2).SetColor(0.0, 0.0, 1.0) cube_axes_actor.DrawXGridlinesOn() cube_axes_actor.DrawYGridlinesOn() cube_axes_actor.DrawZGridlinesOn() if not tick_vis: cube_axes_actor.XAxisMinorTickVisibilityOff() cube_axes_actor.YAxisMinorTickVisibilityOff() cube_axes_actor.ZAxisMinorTickVisibilityOff() cube_axes_actor.SetXTitle("X") cube_axes_actor.SetYTitle("Y") cube_axes_actor.SetZTitle("Z") cube_axes_actor.SetXAxisLabelVisibility(1) cube_axes_actor.SetYAxisLabelVisibility(1) cube_axes_actor.SetZAxisLabelVisibility(1) # only plot grid lines furthest from viewpoint # ensure platform compatibility for the grid line options if sys.platform == "win32": cube_axes_actor.SetGridLineLocation(cube_axes_actor.VTK_GRID_LINES_FURTHEST) else: # rather use elif == "linux" ? but what about other platforms try: # apparently this can also go wrong on linux, maybe depends on vtk version? cube_axes_actor.SetGridLineLocation(vtk.VTK_GRID_LINES_FURTHEST) except AttributeError: pass return cube_axes_actor def update_surfaces_real_time(self, interp_data): lith_block, fault_block = gp.compute_model(interp_data) try: v_l, s_l = gp.get_surfaces(interp_data, lith_block[1], fault_block[1::2], original_scale=False) except IndexError: v_l, s_l = gp.get_surfaces(interp_data, lith_block[1], None, original_scale=False) return v_l, s_l # def load_posteriors(self, dbname): # import gempy.UncertaintyAnalysisPYMC2 # self.post = gempy.UncertaintyAnalysisPYMC2.Posterior(dbname)
[docs] @staticmethod def export_vtk_lith_block(geo_data, lith_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 """ from evtk.hl import gridToVTK import numpy as np # Dimensions nx, ny, nz = geo_data.resolution lx = geo_data.extent[1] - geo_data.extent[0] ly = geo_data.extent[3] - geo_data.extent[2] lz = geo_data.extent[5] - geo_data.extent[4] dx, dy, dz = lx / nx, ly / ny, lz / nz ncells = nx * ny * nz npoints = (nx + 1) * (ny + 1) * (nz + 1) # Coordinates x = np.arange(geo_data.extent[0], geo_data.extent[1] + 0.1, dx, dtype='float64') y = np.arange(geo_data.extent[2], geo_data.extent[3] + 0.1, dy, dtype='float64') z = np.arange(geo_data.extent[4], geo_data.extent[5] + 0.1, dz, dtype='float64') lith = lith_block.reshape((nx, ny, nz)) # Variables if not path: path = "./default" gridToVTK(path+'_lith_block', x, y, z, cellData={"Lithology": lith})
[docs] @staticmethod def export_vtk_surfaces(vertices, simplices, 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 """ import vtk for s_n in range(len(vertices)): # setup points and vertices Points = vtk.vtkPoints() Triangles = vtk.vtkCellArray() Triangle = vtk.vtkTriangle() for p in vertices[s_n]: Points.InsertNextPoint(p) # Unfortunately in this simple example the following lines are ambiguous. # The first 0 is the index of the triangle vertex which is ALWAYS 0-2. # The second 0 is the index into the point (geometry) array, so this can range from 0-(NumPoints-1) # i.e. a more general statement is triangle->GetPointIds()->SetId(0, PointId); for i in simplices[s_n]: Triangle.GetPointIds().SetId(0, i[0]) Triangle.GetPointIds().SetId(1, i[1]) Triangle.GetPointIds().SetId(2, i[2]) Triangles.InsertNextCell(Triangle) polydata = vtk.vtkPolyData() polydata.SetPoints(Points) polydata.SetPolys(Triangles) polydata.Modified() if vtk.VTK_MAJOR_VERSION <= 5: polydata.Update() writer = vtk.vtkXMLPolyDataWriter(); if not path: path = "./default_" writer.SetFileName(path+'_surfaces'+str(s_n)+'vtp') if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(polydata) else: writer.SetInputData(polydata) writer.Write()
def _create_color_lot(geo_data, cd_rgb): """Returns color [r,g,b] LOT for formation numbers.""" if "formation number" not in geo_data.interfaces or "formation number" not in geo_data.foliations: geo_data.set_formation_number() # if not, set formation numbers c_names = ["indigo", "red", "yellow", "brown", "orange", "green", "blue", "amber", "pink", "light-blue", "lime", "blue-grey", "deep-orange", "grey", "cyan", "deep-purple", "purple", "teal", "light-green"] lot = {} ci = 0 # use as an independent running variable because of fault formations # get unique formation numbers fmt_numbers = np.unique([val for val in geo_data.interfaces['formation number'].unique()]) # get unique fault formation numbers fault_fmt_numbers = np.unique(geo_data.interfaces[geo_data.interfaces["isFault"] == True]["formation number"]) # iterate over all unique formation numbers for i, n in enumerate(fmt_numbers): # if its a fault formation set it to black by default if n in fault_fmt_numbers: lot[n] = cd_rgb["black"]["400"] # if not, just go through else: lot[n] = cd_rgb[c_names[ci]]["400"] ci += 1 return lot