Source code for pygimli.viewer.showmesh

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Generic mesh visualization tools."""

import os
import sys
import time
import traceback

import numpy as np

from pygimli.viewer.mpl.colorbar import setMappableData

from .. core.logger import renameKwarg

try:
    import pygimli as pg
    from .showmatrix import showMatrix
    from .mpl import drawMesh, drawModel, drawField, drawSensors, drawStreams
    from .mpl import drawSelectedMeshBoundaries
    from .mpl import addCoverageAlpha
    from .mpl import updateAxes  # , registerShowPendingFigsAtExit
    from .mpl import createColorBar, updateColorBar
    from .mpl import CellBrowser
    from .mpl.colorbar import cmapFromName
except ImportError as e:
    print(e)
    traceback.print_exc(file=sys.stdout)
    pg.critical("ERROR: cannot import the library 'pygimli'."
                "Ensure that pygimli is in your PYTHONPATH ")


[docs] def show(obj=None, data=None, **kwargs): """Mesh and model visualization. Syntactic sugar to show a obj with data. Forwards to a known visualization for obj. Typical is :py:mod:`pygimli.viewer.showMesh` or :py:mod:`pygimli.viewer.mayaview.showMesh3D` to show most of the typical 2D and 3D content. See tutorials and examples for usage hints. An empty show call creates an empty ax window. Parameters ---------- obj: obj obj can be so far. * :gimliapi:`GIMLI::Mesh` or list of meshes * DataContainer * pg.core.Sparse[Map]Matrix data: iterable Optionally data to visualize. See appropriate show function. Keyword Arguments ----------------- **kwargs Additional kwargs forward to appropriate show functions. * ax : axe [None] Matplotlib axes object. Create a new if necessary. * fitView : bool [True] Scale x and y limits to match the view. Returns ------- Return the results from the showMesh* functions. Usually the axe object and a colorbar. See Also -------- showMesh """ if "axes" in kwargs: # remove me in 1.2 #20200515 print("Deprecation Warning: Please use keyword `ax` instead of `axes`") kwargs['ax'] = kwargs.pop('axes', None) # Empty call just to create an mpl axes # if obj is None and 'mesh' not in kwargs: if obj is None and 'mesh' not in kwargs.keys(): ax = kwargs.pop('ax', None) if ax is None: ax = pg.plt.subplots(figsize=kwargs.pop('figsize', None))[1] return ax, None # registerShowPendingFigsAtExit() # try to check if obj containes a mesh if hasattr(obj, 'mesh'): return pg.show(obj.mesh, obj, **kwargs) # try to interpret obj as ERT Data if isinstance(obj, pg.DataContainerERT): from pygimli.physics.ert import showERTData return showERTData(obj, vals=kwargs.pop('vals', data), **kwargs) # try to interpret obj as matrices if isinstance(obj, pg.core.MatrixBase) or (isinstance(obj, np.ndarray) and obj.ndim == 2): return showMatrix(obj, **kwargs) try: from scipy.sparse import spmatrix if isinstance(obj, spmatrix): return showMatrix(obj, **kwargs) except ImportError: pass # try to interprete obj as mesh or list of meshes mesh = kwargs.pop('mesh', obj) fitView = kwargs.get('fitView', True) if isinstance(mesh, list): ax = kwargs.pop('ax', None) ax, cBar = show(mesh[0], data, hold=True, ax=ax, fitView=fitView, **kwargs) for m in mesh[1:]: ax, cBar = show(m, data, ax=ax, hold=True, fitView=False, **kwargs) if fitView is not False: ax.autoscale(enable=True, axis='both', tight=True) ax.set_aspect('equal') return ax, cBar if isinstance(mesh, pg.Mesh): if isinstance(data, str): if data in mesh.dataKeys(): data = mesh[data] # elif 0: # maybe check x, y, z, cellMarker etc. else: pg.error(f"Could not retrieve data from key {data}") return None, None if mesh.dim() == 2: if pg.zero(pg.y(mesh)): pg.info("swap z<->y coordinates for visualization.") meshSwap = pg.Mesh(mesh) for n in meshSwap.nodes(): n.pos()[1] = n.pos()[2] return showMesh(meshSwap, data, **kwargs) return showMesh(mesh, data, **kwargs) elif mesh.dim() == 3: from .pv import showMesh3D return showMesh3D(mesh, data, **kwargs) else: pg.error("ERROR: Mesh not valid.", mesh) if isinstance(obj, pg.core.Boundary): ax = kwargs.pop('ax', None) drawSelectedMeshBoundaries(ax, [obj], **kwargs) return ax, None pg.error("Can't interprete obj: {0} to show.".format(obj)) return None, None
[docs] def showMesh(mesh, data=None, block=False, colorBar=None, label=None, coverage=None, ax=None, savefig=None, showMesh=False, showBoundary=None, markers=False, **kwargs): """2D Mesh visualization. Create an axis object and plot a 2D mesh with given node or cell data. Returns the axis and the color bar. The type of data determines the appropriate draw method. Parameters ---------- mesh: :gimliapi:`GIMLI::Mesh` 2D or 3D GIMLi mesh data: iterable [None] Optionally data to visualize. . None (draw mesh only) forward to :py:mod:`pygimli.viewer.mpl.drawMesh` or if no cells are given: forward to :py:mod:`pygimli.viewer.mpl.drawPLC` . [[marker, value], ...] List of Cellvalues per cell marker forward to :py:mod:`pygimli.viewer.mpl.drawModel` . float per cell -- model, patch forward to :py:mod:`pygimli.viewer.mpl.drawModel` . float per node -- scalar field forward to :py:mod:`pygimli.viewer.mpl.drawField` . iterable of type [float, float] -- vector field forward to :py:mod:`pygimli.viewer.mpl.drawStreams` . pg.PosVector -- vector field forward to :py:mod:`pygimli.viewer.mpl.drawStreams` . pg.core.stdVectorRVector3 -- sensor positions forward to :py:mod:`pygimli.viewer.mpl.drawSensors` block: bool [False] Force to open the Figure of your content and blocks the script until you close the current figure. Same like pg.show(); pg.wait() colorBar: bool [None], Colorbar Create and show a colorbar. If colorBar is a valid colorbar then only its values will be updated. label: str Set colorbar label. If set colorbar is toggled to True. [None] coverage: iterable [None] Weight data by the given coverage array and fadeout the color. ax: matplotlib.Axes [None] Instead of creating a new and empty ax, just draw into the given one. Useful to combine multiple plots into one figure. savefig: string Filename for a direct save to disc. showMesh: bool [False] Shows the mesh itself additional. showBoundary: bool [None] Highlight all boundaries with marker != 0. None means automatic. True for cell data and False for node data. marker: bool [False] Show cell markers and boundary marker. boundaryMarkers: bool [False] Highlight boundaries with marker !=0 and add Marker annotation. Applies :py:mod:`pygimli.viewer.mpl.drawBoundaryMarkers`. Dictionary "boundaryProps" can be added and will be forwarded to :py:mod:`pygimli.viewer.mpl.drawBoundaryMarkers`. Keyword Arguments ----------------- xlabel: str [None] Add label to the x axis ylabel: str [None] Add label to the y axis fitView: bool Fit the axes limits to the all content of the axes. Default True. boundaryProps: dict Arguments for plotboundar hold: bool [pg.hold()] Holds back the opening of the Figure. If set to True [default] nothing happens until you either force another show with hold=False or block=True or call pg.wait() or pg.plt.show(). If hold is set to False your script will open the figure and continue working. You can change global hold with pg.hold(bool). All remaining will be forwarded to the draw functions and matplotlib methods, respectively. Examples -------- >>> import pygimli as pg >>> import pygimli.meshtools as mt >>> world = mt.createWorld(start=[-10, 0], end=[10, -10], ... layers=[-3, -7], worldMarker=False) >>> mesh = mt.createMesh(world, quality=32, area=0.2, smooth=[1, 10]) >>> _ = pg.viewer.showMesh(mesh, markers=True) Returns ------- ax : matplotlib.axes cBar : matplotlib.colorbar """ renameKwarg('cmap', 'cMap', kwargs) cMap = kwargs.pop('cMap', 'viridis') cBarOrientation = kwargs.pop('orientation', 'horizontal') replaceData = kwargs.pop('replaceData', False) if ax is None: ax, _ = pg.show(figsize=kwargs.pop('figsize', None), **kwargs) # adjust limits only when axis is empty fitViewDefault = True # if (ax.lines or ax.collections or ax.patches): # fitViewDefault = False # else: # plt.subplots() resets locale setting to system default .. this went # horrible wrong for german 'decimal_point': ',' pg.checkAndFixLocaleDecimal_point(verbose=False) hold = kwargs.pop('hold', pg.viewer.mpl.utils.__holdAxes__) if block is True: hold = True lastHoldStatus = pg.viewer.mpl.utils.__holdAxes__ pg.viewer.mpl.hold(val=hold) gci = None validData = False if markers: kwargs["boundaryMarkers"] = kwargs.get("boundaryMarkers", True) if mesh.cellCount() > 0: uniquemarkers, uniqueidx = np.unique( np.array(mesh.cellMarkers()), return_inverse=True) label = label or "Cell markers" cMap = cmapFromName("Set3", ncols=len(uniquemarkers)) #cMap = pg.plt.colormaps.get_cmap("Set3", len(uniquemarkers)) kwargs["logScale"] = False kwargs["cMin"] = -0.5 kwargs["cMax"] = len(uniquemarkers) - 0.5 data = np.arange(len(uniquemarkers))[uniqueidx] if data is None: showMesh = True mesh.createNeighborInfos() if showBoundary is None: showBoundary = True elif isinstance(data, pg.core.stdVectorRVector3): drawSensors(ax, data, **kwargs) elif isinstance(data, pg.PosVector): drawStreams(ax, mesh, data, **kwargs) else: # check for map like data=[[marker, val], ....] if isinstance(data, list) and \ isinstance(data[0], list) and isinstance(data[0][0], int): data = pg.solver.parseMapToCellArray(data, mesh) if hasattr(data[0], '__len__') and not \ isinstance(data, np.ma.core.MaskedArray) and not \ isinstance(data[0], str): # [u,v] x N if len(data) == 2: data = np.array(data).T # N x [u,v] if data.shape[1] == 2: drawStreams(ax, mesh, data, **kwargs) # N x [u,v,w] elif data.shape[1] == 3: # if sum(data[:, 0]) != sum(data[:, 1]): # drawStreams(ax, mesh, data, **kwargs) drawStreams(ax, mesh, data[:, :2], **kwargs) else: # Try animation frames x N data = np.asarray(data) if data.ndim == 2: if data.shape[1] == mesh.cellCount() or \ data.shape[1] == mesh.nodeCount(): return showAnimation(mesh, data, cMap=cMap, ax=ax, **kwargs) pg.warn("No valid stream data:", data.shape, data.ndim) showMesh = True # elif min(data) == max(data): # or pg.core.haveInfNaN(data): # pg.warn("No valid data: ", min(data), max(data), # pg.core.haveInfNaN(data)) # showMesh = True else: if bool(colorBar) is not False: colorBar = True if len(data) == mesh.cellCount(): if showBoundary is None: showBoundary = True def _drawField(ax, mesh, data, kwargs): # kwargs as reference here to set defaults valid outside too validData = True if len(data) == mesh.cellCount(): kwargs['nCols'] = kwargs.pop('nCols', 256) gci = drawModel(ax, mesh, data, **kwargs) elif len(data) == mesh.nodeCount(): kwargs['nLevs'] = kwargs.pop('nLevs', 5) kwargs['nCols'] = kwargs.pop('nCols', kwargs['nLevs']-1) gci = drawField(ax, mesh, data, **kwargs) else: pg.error("Data size invalid") print("Data: ", len(data), min(data), max(data), pg.core.haveInfNaN(data)) print("Mesh: ", mesh) validData = False drawMesh(ax, mesh) return gci, validData try: if label is None: label = "" if replaceData and hasattr(mesh, 'gci') and ax in mesh.gci: gci = mesh.gci[ax] if 'TriContourSet' in str(type(gci)): ax.clear() gci, validData = _drawField(ax, mesh, data, kwargs) updateAxes(ax, force=True) else: setMappableData(gci, data, cMin=kwargs.get('cMin', None), cMax=kwargs.get('cMax', None),) updateAxes(ax, force=True) return ax, gci.colorbar else: gci, validData = _drawField(ax, mesh, data, kwargs) # Cache mesh and scalarmappable to make replaceData work if not hasattr(mesh, 'gci'): mesh.gci = {} mesh.gci[ax] = gci if cMap is not None and gci is not None: gci.set_cmap(cmapFromName(cMap)) # gci.cmap.set_under('k') except BaseException as e: print(e) traceback.print_exc(file=sys.stdout) if mesh.cellCount() == 0: showMesh = False if mesh.boundaryCount() == 0: pg.viewer.mpl.drawPLC(ax, mesh, showNodes=True, fillRegion=False, showBoundary=False, **kwargs) showBoundary = False # ax.plot(pg.x(mesh), pg.y(mesh), '.', color='black') else: kwargs['orientation'] = cBarOrientation pg.viewer.mpl.drawPLC(ax, mesh, **kwargs) if showMesh: if gci is not None and hasattr(gci, 'set_antialiased'): gci.set_antialiased(True) gci.set_linewidth(0.3) gci.set_edgecolor(kwargs.pop('color', "0.1")) #drawMesh(ax, mesh, lw=0.3, **kwargs) #else: drawMesh(ax, mesh, lw=0.3, **kwargs) # pg.viewer.mpl.drawSelectedMeshBoundaries(ax, # mesh.boundaries(), # color=kwargs.pop('color', "0.1"), # linewidth=kwargs.pop('lw', 0.3)) if bool(showBoundary) is True: b = mesh.boundaries(mesh.boundaryMarkers() != 0) pg.viewer.mpl.drawSelectedMeshBoundaries(ax, b, color=(0.0, 0.0, 0.0, 1.0), linewidth=1.4) if kwargs.pop("boundaryMarkers", False): pg.viewer.mpl.drawBoundaryMarkers( ax, mesh, clipBoundaryMarkers=kwargs.pop('clipBoundaryMarkers', False), **kwargs.pop('boundaryProps', {})) fitView = kwargs.pop('fitView', fitViewDefault) if fitView is not False: ax.autoscale(enable=True, axis='both', tight=True) ax.set_aspect(kwargs.pop('aspect', 'equal')) cBar = None if label is not None and colorBar is None: colorBar = True # pg._r(validData, gci) if validData: labels = ['cMin', 'cMax', 'nCols', 'nLevs', 'logScale', 'levels'] subkwargs = {key: kwargs[key] for key in labels if key in kwargs} subkwargs['cMap'] = cMap if isinstance(colorBar, bool): if colorBar is True: subkwargs['label'] = label subkwargs['orientation'] = cBarOrientation cBar = createColorBar(gci, size=kwargs.pop('size', 0.2), pad=kwargs.pop('pad', None), **subkwargs, onlyColorSet=not colorBar, ) elif colorBar is not False: cBar = updateColorBar(colorBar, **subkwargs) if markers and cBar is not None: ticks = np.arange(len(uniquemarkers)) cBar.set_ticks(ticks) labels = [] for marker in uniquemarkers: labels.append(str((marker))) cBar.set_ticklabels(labels) if coverage is not None: if isinstance(coverage, (float, int)): gci.set_alpha(coverage) elif len(data) == len(coverage) == mesh.cellCount(): addCoverageAlpha(gci, coverage, dropThreshold=kwargs.pop('dropThreshold', 0.4)) else: pg.error('Coverage needs to be either of type float or an array', 'with the same length as data and mesh.cellCount().') # addCoverageAlpha(gci, pg.core.cellDataToPointData(mesh, # coverage)) if not hold or block is not False and \ pg.plt.get_backend().lower() != "agg": if data is not None: if len(data) == mesh.cellCount(): CellBrowser(mesh, data, ax=ax) pg.plt.show(block=block) try: pg.plt.pause(0.01) except BaseException: pass pg.viewer.mpl.updateAxes(ax) pg.viewer.mpl.hold(val=lastHoldStatus) if savefig: print('saving: ' + savefig + ' ...', end="") if '.' not in savefig: savefig += '.pdf' ax.figure.savefig(savefig, bbox_inches='tight') # rc params savefig.format=pdf print('.. done') return ax, cBar
def showBoundaryNorm(mesh, normMap=None, **kwargs): """Show mesh boundaries normals. Show the mesh and draw a black line along the normal direction of all boundaries. If you provide a boundary marker vs. norm direction map, then only these norms are drawn. Parameters ---------- mesh : :gimliapi:`GIMLI::Mesh` 2D or 3D GIMLi mesh normMap : list list of [boundary marker, [norm]] pairs. e.g. [[1, [0.0,1.0]], ... ] **kwargs : Will be forwarded to the draw functions and matplotlib methods, respectively. Returns ------- ax : matplotlib.ax """ ax = kwargs.pop('ax', None) col = kwargs.pop('color', 'Black') if normMap: for pair in normMap: bounds = mesh.findBoundaryByMarker(pair[0]) for b in bounds: c1 = b.center() if (pair[1][0] != 0) or (pair[1][1] != 0): ax.arrow(c1[0], c1[1], pair[1][0], pair[1][1], head_width=0.1, head_length=0.3, color=col, **kwargs) else: ax.plot(c1[0], c1[1], 'o', color=col) return ax = show(mesh, hold=True, ax=ax)[0] for b in mesh.boundaries(): c1 = b.center() c2 = c1 + b.norm() ax.plot([c1[0], c2[0]], [c1[1], c2[1]], color=col, **kwargs) time.sleep(0.05) return ax __Animation_Keeper__ = None def showAnimation(mesh, data, ax=None, **kwargs): """Show timelapse mesh data. Time will be annotated if the mesh contains a valid 'times' data array. Note, there can be only one animation per time. Best viewed in a notebook, because of caching and better animation control elements. TODO ---- * 3D * allow for multiple animations per script Parameters ---------- mesh: :gimliapi:`GIMLI::Mesh` 2D GIMLi mesh data: [NxM] iterable Data matrix for N frames with suitable data size. Keyword Args ------------ dpi : int[96] Movie resolution. Rest is forwarded to :py:func:pygimli.viewer.show: """ import matplotlib.animation plt = pg.plt plt.rcParams["animation.html"] = "jshtml" plt.rcParams['figure.dpi'] = kwargs.pop('dpi', 96) plt.rcParams['animation.embed_limit'] = 50 figsize = kwargs.pop("figsize", None) flux = kwargs.pop('flux', None) plt.ioff() plc = kwargs.pop("plc", None) pg.show(mesh, data[0], ax=ax, figsize=figsize, **kwargs) if flux is not None: pg.show(mesh, flux[0], ax=ax) try: times = mesh['times'] except Exception as e: times = None p = pg.utils.ProgressBar(len(data)) def animate(t): p.update(t) ax.clear() pg.show(mesh, data[t], ax=ax, **kwargs) if flux is not None: pg.show(mesh, flux[t], ax=ax) if plc is not None: pg.viewer.mpl.drawMesh(ax, plc, fillRegion=False, fitView=False) if times is not None and len(times) > t: # ax.text(0.02, 0.02, f't={pg.pf(times[t])}', ax.text(0.01, 1.01, f't={pg.utils.prettyTime(times[t])}', horizontalalignment='left', verticalalignment='bottom', transform=ax.transAxes, color='k', fontsize=8) if pg.isNotebook() is False: global __Animation_Keeper__ __Animation_Keeper__ = matplotlib.animation.FuncAnimation(ax.figure, animate, frames=len(data)) return __Animation_Keeper__