diff --git a/packages/loopstructural_visualisation/CHANGELOG.md b/packages/loopstructural_visualisation/CHANGELOG.md new file mode 100644 index 000000000..132ac2c19 --- /dev/null +++ b/packages/loopstructural_visualisation/CHANGELOG.md @@ -0,0 +1,141 @@ +# Changelog + +## [0.1.17](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.16...v0.1.17) (2025-08-14) + + +### Bug Fixes + +* linting ([87257f4](https://github.com/Loop3D/loopstructural-visualisation/commit/87257f43b8915061c7f8e3334ded449feaad35d4)) +* update for new stratigraphic colum in LoopStructural ([56def61](https://github.com/Loop3D/loopstructural-visualisation/commit/56def61d639492d5d60ccc9bfb74c54e9bafadd6)) + +## [0.1.16](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.15...v0.1.16) (2025-05-27) + + +### Bug Fixes + +* Change scalar_bar arg to show_scalar_bar ([949b3eb](https://github.com/Loop3D/loopstructural-visualisation/commit/949b3eb23c3bf38b67a9cf1f9a332aa3defdac7b)) + +## [0.1.15](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.14...v0.1.15) (2025-04-17) + + +### Bug Fixes + +* release 0.1.15 ([2e4f5c9](https://github.com/Loop3D/loopstructural-visualisation/commit/2e4f5c9628177916b7174dd48e8a321534a116b1)) + +## [0.1.14](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.13...v0.1.14) (2025-02-20) + + +### Bug Fixes + +* scaling of vector fields is now consistent for different scale models. ([49450bb](https://github.com/Loop3D/loopstructural-visualisation/commit/49450bb0514e58f3da5ee3af98bd307a3fc8c4d8)) + +## [0.1.13](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.12...v0.1.13) (2025-01-15) + + +### Bug Fixes + +* Stratigraphic column visualiser display in correct order ([7dcd5b8](https://github.com/Loop3D/loopstructural-visualisation/commit/7dcd5b8c06b637d33ac2e7c5e02bb52a1991cfdd)) +* Stratigraphic column visualiser display in correct order ([7dcd5b8](https://github.com/Loop3D/loopstructural-visualisation/commit/7dcd5b8c06b637d33ac2e7c5e02bb52a1991cfdd)) + +## [0.1.12](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.11...v0.1.12) (2024-12-18) + + +### Bug Fixes + +* add import check so trame is not required ([d6aa391](https://github.com/Loop3D/loopstructural-visualisation/commit/d6aa391257604026261c45d22e6cc64792fc669c)) + +## [0.1.11](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.10...v0.1.11) (2024-12-18) + + +### Bug Fixes + +* adding disk option to fault visualisation ([1fd1a9d](https://github.com/Loop3D/loopstructural-visualisation/commit/1fd1a9d3843e1266f4af03d6572558fa65f53527)) +* require ls >1.6.4 and add fault ellipsoid ([e71b0ab](https://github.com/Loop3D/loopstructural-visualisation/commit/e71b0abd1b4e242582214c57d485cf2939b90899)) +* update bug with trame folder ([edd3188](https://github.com/Loop3D/loopstructural-visualisation/commit/edd318850fc7e97ef4c5e37d43511f9476046324)) + +## [0.1.10](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.9...v0.1.10) (2024-10-23) + + +### Bug Fixes + +* catch case where name is passed but is None ([52ebf0f](https://github.com/Loop3D/loopstructural-visualisation/commit/52ebf0f0fe6ff253e782f8bc2832717df573e0fe)) + +## [0.1.9](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.8...v0.1.9) (2024-10-23) + + +### Bug Fixes + +* add subdir to package ([2663673](https://github.com/Loop3D/loopstructural-visualisation/commit/2663673e8bb55d28b754ab9eeb9d6ed277b2e2d8)) +* ignore subplots and fix invalid actor names ([0d0bc20](https://github.com/Loop3D/loopstructural-visualisation/commit/0d0bc2064b85bb83fa2c3a830d48598f9ab12ef1)) +* random colour was calling invalid rng function ([0182b46](https://github.com/Loop3D/loopstructural-visualisation/commit/0182b46a956a31b533934bc857d64eb5a1dc946c)) +* return cmap not array of colours ([6d0a408](https://github.com/Loop3D/loopstructural-visualisation/commit/6d0a408a8486c05c0f39ab5411f46253c3d92694)) + +## [0.1.8](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.7...v0.1.8) (2024-10-22) + + +### Bug Fixes + +* avoiding cycle ([30c8142](https://github.com/Loop3D/loopstructural-visualisation/commit/30c8142d36629db9675bf93900d5ac8477745e19)) + +## [0.1.7](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.6...v0.1.7) (2024-08-06) + + +### Bug Fixes + +* adding v ([696f407](https://github.com/Loop3D/loopstructural-visualisation/commit/696f407447d140e1e9e6b70150537fd013edf4e3)) +* removing print ([95b007d](https://github.com/Loop3D/loopstructural-visualisation/commit/95b007d5d1eda8bbbd3fb9df4c0171e8128dd3ec)) +* restrict to LS>=1.6.0 and adding conda build ([fcc9ca3](https://github.com/Loop3D/loopstructural-visualisation/commit/fcc9ca3cd458195e157a0bbefc1f5b985e300b84)) + +## [0.1.6](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.5...v0.1.6) (2024-06-06) + + +### Bug Fixes + +* update for ls scalar field changes ([553cd89](https://github.com/Loop3D/loopstructural-visualisation/commit/553cd89f978f2a025ae446cc3e4a6f40cb24e167)) + +## [0.1.5](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.4...v0.1.5) (2024-06-04) + + +### Bug Fixes + +* addin old loopstructural plots ([e716096](https://github.com/Loop3D/loopstructural-visualisation/commit/e7160967272fabb172c43a5c7a476540bdbd2715)) + +## [0.1.4](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.3...v0.1.4) (2024-06-04) + + +### Bug Fixes + +* adding ls logger ([0914d43](https://github.com/Loop3D/loopstructural-visualisation/commit/0914d439ef1606dfb57a3f1e283b991a5149d9f8)) +* adding random colour generator and rbga to hex ([ccc3994](https://github.com/Loop3D/loopstructural-visualisation/commit/ccc3994f44d34d069dcb958b4f45e6a50a652cce)) +* adding threshold to block model ([8c3c0c6](https://github.com/Loop3D/loopstructural-visualisation/commit/8c3c0c62ab1a6565892098cf9fc591e4a817f408)) +* updating to .vtk as a method ([b5ab78c](https://github.com/Loop3D/loopstructural-visualisation/commit/b5ab78c2503a59578133ad75bd329a87c10846f8)) + +## [0.1.3](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.2...v0.1.3) (2024-05-31) + + +### Bug Fixes + +* add wrapper to use pyvista plane slicer ([7bbce47](https://github.com/Loop3D/loopstructural-visualisation/commit/7bbce472083308e76828beaa867170ba6ebdafb9)) +* make sure pyvista has from_regular_faces ([6a1d402](https://github.com/Loop3D/loopstructural-visualisation/commit/6a1d402f3637d986b7e56e316f030e0ff3629e1f)) +* rename scale bar to scalar_bar ([7025441](https://github.com/Loop3D/loopstructural-visualisation/commit/7025441c74ed6518debe63fa279924d3873183bf)) + +## [0.1.2](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.1...v0.1.2) (2024-05-28) + + +### Bug Fixes + +* update pypi to remove username and pw ([c8d1b43](https://github.com/Loop3D/loopstructural-visualisation/commit/c8d1b430443a2c454b5aeeff3f2099f8d5a3a978)) + +## [0.1.1](https://github.com/Loop3D/loopstructural-visualisation/compare/v0.1.0...v0.1.1) (2024-05-28) + + +### Bug Fixes + +* trigger release ([fc6a70f](https://github.com/Loop3D/loopstructural-visualisation/commit/fc6a70fc18ecd43fedd7a3f3d8d717ae790800e7)) + +## 0.1.0 (2024-05-28) + + +### Bug Fixes + +* adding _2d_viewer file ([099c43c](https://github.com/Loop3D/loopstructural-visualisation/commit/099c43c693fb7f5e3c3f1cee587e80c12694411b)) diff --git a/packages/loopstructural_visualisation/LICENSE b/packages/loopstructural_visualisation/LICENSE new file mode 100644 index 000000000..0b251d90d --- /dev/null +++ b/packages/loopstructural_visualisation/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 Loop3D + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/packages/loopstructural_visualisation/README.md b/packages/loopstructural_visualisation/README.md new file mode 100644 index 000000000..91ab10f9b --- /dev/null +++ b/packages/loopstructural_visualisation/README.md @@ -0,0 +1,7 @@ +# loopstructural-visualisation + +A LoopStructural interface for pyvista's Plotter class. + + +To install `pip install loopstructuralvisualisation` or for a jupyter notebook environment (including the pyvista[jupyuter] dependencies) `pip install loopstructuralvisualisation[jupyter]` + diff --git a/packages/loopstructural_visualisation/pyproject.toml b/packages/loopstructural_visualisation/pyproject.toml new file mode 100644 index 000000000..933c7ad5f --- /dev/null +++ b/packages/loopstructural_visualisation/pyproject.toml @@ -0,0 +1,50 @@ +[build-system] +requires = ["setuptools"] +build-backend = "setuptools.build_meta" + +[project] +name = "loopstructuralvisualisation" +description = "3D geological modelling visualisation for LoopStructural" +dynamic = ["version"] +requires-python = ">=3.9" +authors = [{ name = "Lachlan Grose", email = "lachlan.grose@monash.edu" }] +readme = "README.md" +license = { text = "MIT" } +keywords = [ + "earth sciences", + "geology", + "3-D modelling", + "structural geology", + "uncertainty", +] +classifiers = [ + "Development Status :: 5 - Production/Stable", + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering :: Information Analysis", + "License :: OSI Approved :: MIT License", + "Operating System :: Microsoft :: Windows", + "Operating System :: POSIX", + "Operating System :: MacOS", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", +] +dependencies = ["numpy>=1.18", "pyvista>=0.42", "LoopStructural>=1.6.17"] + +[project.optional-dependencies] +jupyter = ["pyvista[jupyter]"] +all = ["pyvista[all]"] +tests = ["pytest"] + +[project.urls] +Documentation = "https://Loop3d.org/LoopStructural/" +"Bug Tracker" = "https://github.com/loop3d/loopstructural-visualisation/issues" +"Source Code" = "https://github.com/loop3d/loopstructural-visualisation" + +[tool.setuptools.dynamic] +version = { attr = "loopstructuralvisualisation.version.__version__" } + +[tool.setuptools.packages.find] +where = ["src"] +include = ["loopstructuralvisualisation", "loopstructuralvisualisation.*"] diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_2d_viewer.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_2d_viewer.py new file mode 100644 index 000000000..23cb01101 --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_2d_viewer.py @@ -0,0 +1,369 @@ +import matplotlib.pyplot as plt +import numpy as np + +from LoopStructural.utils import getLogger + +logger = getLogger(__name__) +from LoopStructural.modelling.features import FeatureType + + +class Loop2DView: + """ """ + + def __init__(self, model=None, bounding_box=np.zeros((2, 2)), nsteps=None, ax=None, **kwargs): + """ + + Parameters + ---------- + origin - lower left + maximum - upper right + nsteps - number of cells + kwargs + """ + + self.xx = None + self.yy = None + + self._bounding_box = bounding_box + self._nsteps = nsteps + if self._nsteps is not None and self._bounding_box is not None: + self._update_grid() + if model is not None: + # make sure self._nsteps is 2d + self.model = model + self.ax = ax + if self.ax is None: + fig, self.ax = plt.subplots(1, figsize=(10, 10)) + self.ax.set_aspect("equal", adjustable="box") + + # set plot limits to model bounding box + self._xmin = self.bounding_box[0, 0] + self._xmax = self.bounding_box[1, 0] + self._ymin = self.bounding_box[0, 1] + self._ymax = self.bounding_box[1, 1] + self._update_plot_limits() + + @property + def model(self): + return self._model + + @model.setter + def model(self, model): + if model is not None: + bb = np.array([model.origin[:2], model.maximum[:2]]) + self.bounding_box = bb # model.bounding_box + self.nsteps = model.nsteps[:2] + self._model = model + self._update_grid() + + @property + def nsteps(self): + return self._nsteps + + @nsteps.setter + def nsteps(self, nsteps): + if len(nsteps) != 2: + logger.error("Can't update nsteps, needs to be 2D") + return + self._nsteps = nsteps + self._update_grid() + + @property + def bounding_box(self): + return self._bounding_box + + @bounding_box.setter + def bounding_box(self, bounding_box): + self._bounding_box = bounding_box + self._update_grid() + + @property + def xmin(self): + return self._xmin + + @xmin.setter + def xmin(self, xmin): + self._xmin = xmin + self._update_plot_limits() + + @property + def xmax(self): + return self._xmax + + @xmax.setter + def xmax(self, xmax): + self._xmax = xmax + self._update_plot_limits() + + @property + def ymin(self): + return self._ymin + + @ymin.setter + def ymin(self, ymin): + self._ymin = ymin + self._update_plot_limits() + + @property + def ymax(self): + return self._ymax + + @ymax.setter + def ymax(self, ymax): + self._ymax = ymax + self._update_plot_limits() + + def _update_plot_limits(self): + self.ax.set_xlim([self._xmin, self._xmax]) + self.ax.set_ylim([self._ymin, self._ymax]) + + def _update_grid(self): + """Internal function to update the current grid when the bounding box + or number of steps changes + """ + if self.nsteps is None or self.bounding_box is None: + return + x = np.linspace(self.bounding_box[0, 0], self.bounding_box[1, 0], self.nsteps[0]) + y = np.linspace(self.bounding_box[0, 1], self.bounding_box[1, 1], self.nsteps[1]) + self.xx, self.yy = np.meshgrid(x, y, indexing="ij") + self.xx = self.xx.flatten() + self.yy = self.yy.flatten() + + def add_data(self, feature, val=True, grad=True, unfault=False, dip=True, **kwargs): + """ + Adds the data associated to the feature to the plot + Parameters + ---------- + feature : GeologicalFeature + the feature whose data you want to add + val : bool + whether to add value data + grad : bool + whether to add gradient data + unfault : bool + plot points in their restored location + dip : bool + whether to annotate the dip, default False + + kwargs are passed to matplotlib functions and draw strike + Returns + ------- + + """ + # logger.warning("Plotting restored data locations") + ori_data = [] + gradient_data = feature.builder.get_gradient_constraints() + if unfault: + gradient_data = feature.interpolator.get_gradient_constraints() + + if gradient_data.shape[0] > 0: + ori_data.append(gradient_data) + norm_data = feature.builder.get_norm_constraints() + if unfault: + norm_data = feature.interpolator.get_norm_constraints() + + if norm_data.shape[0] > 0: + ori_data.append(norm_data) + cmap = kwargs.pop("cmap", "rainbow") + # if single colour then specify kwarg, otherwise use point value + if val: + value_data = np.copy(feature.builder.get_value_constraints()) + if unfault: + value_data = np.copy(feature.interpolator.get_value_constraints()) + + value_data[:, :3] = self.model.rescale(value_data[:, :3], inplace=False) + point_colour = kwargs.pop("point_colour", None) + if point_colour is None: + self.ax.scatter( + value_data[:, 0], + value_data[:, 1], + c=value_data[:, 3], + vmin=feature.min(), + vmax=feature.max(), + cmap=cmap, + ) + if point_colour is not None: + self.ax.scatter(value_data[:, 0], value_data[:, 1], c=point_colour) + if grad: + symb_colour = kwargs.pop("symb_colour", "black") + symb_scale = kwargs.pop("symb_scale", 1.0) + gradient_data = np.hstack(ori_data) + gradient_data[:, :3] = self.model.rescale(gradient_data[:, :3], inplace=False) + gradient_data[:, 3:5] /= np.linalg.norm(gradient_data[:, 3:5], axis=1)[:, None] + t = gradient_data[:, [4, 3]] * np.array([1, -1]).T + n = gradient_data[:, 3:5] + t *= symb_scale + n *= 0.5 * symb_scale + p1 = gradient_data[:, [0, 1]] - t + p2 = gradient_data[:, [0, 1]] + t + # plt.scatter(val[:,0],val[:,1],c='black') + self.ax.plot([p1[:, 0], p2[:, 0]], [p1[:, 1], p2[:, 1]], symb_colour) + p1 = gradient_data[:, [0, 1]] + p2 = gradient_data[:, [0, 1]] + n + self.ax.plot([p1[:, 0], p2[:, 0]], [p1[:, 1], p2[:, 1]], symb_colour) + if dip: + dip_v = np.rad2deg(np.arccos(gradient_data[:, 5])).astype(int) + for d, xy, v in zip(dip_v, gradient_data[:, :2], gradient_data[:, 3:6]): + self.ax.annotate(d, xy, xytext=xy + v[:2] * symb_scale * 0.1, fontsize="small") + + def add_fault_ellipse(self, faults=None, **kwargs): + from matplotlib.patches import Ellipse + + for f in self.model.stratigraphic_column["faults"].values(): + center = self.model.rescale(f["FaultCenter"]) + e = Ellipse( + (center[0], center[1]), + f["HorizontalRadius"] * 2, + f["InfluenceDistance"] * 2, + 360 - f["FaultDipDirection"], + facecolor="None", + edgecolor="k", + ) + + self.ax.add_patch(e) + + def add_scalar_field(self, feature, z=0, **kwargs): + """ + Plot the scalar field value on a map + + Parameters + ---------- + feature : GeologicalFeature + which feature to plot on the map + z : double/np.array + height + kwargs + + Returns + ------- + + """ + zz = np.zeros(self.xx.shape) + zz[:] = z + v = feature.evaluate_value( + self.model.scale(np.array([self.xx, self.yy, zz]).T, inplace=False) + ) + return self.ax.imshow( + v.reshape(self.nsteps).T, + extent=[ + self.bounding_box[0, 0], + self.bounding_box[1, 0], + self.bounding_box[0, 1], + self.bounding_box[1, 1], + ], + vmin=feature.min(), + vmax=feature.max(), + origin="lower", + **kwargs, + ) + + def add_contour(self, feature, values, z=0, mask=None, **kwargs): + """Add an isoline of a scalar field to the map + + Parameters + ---------- + feature : GeologicalFeature + the feature to isosurface + values : list + list of values to contour + z : double/np.array, optional + elevation of map, by default 0 + """ + zz = np.zeros(self.xx.shape) + zz[:] = z + v = feature.evaluate_value( + self.model.scale(np.array([self.xx, self.yy, zz]).T, inplace=False) + ) + if mask: + maskv = mask(self.model.scale(np.array([self.xx, self.yy, zz]).T, inplace=False)) + v[~maskv] = np.nan + return self.ax.contour( + v.reshape(self.nsteps).T, + extent=[ + self.bounding_box[0, 0], + self.bounding_box[1, 0], + self.bounding_box[0, 1], + self.bounding_box[1, 1], + ], + origin="lower", + levels=values, + **kwargs, + ) + + def add_model(self, z=0, cmap=None): + """Plot the model onto a map + + Parameters + ---------- + z : int/numpy array, optional + height of the map surface (could also be a dem), by default 0 + cmap : str/matplotlib colourmap, optional + specify a colour map, by default 'tab20' + """ + if cmap is None: + import matplotlib.colors as colors + + colours = [] + boundaries = [] + data = [] + for g in self.model.stratigraphic_column.keys(): + if g == "faults": + continue + for v in self.model.stratigraphic_column[g].values(): + data.append((v["id"], v["colour"])) + colours.append(v["colour"]) + boundaries.append(v["id"]) # print(u,v) + cmap = colors.ListedColormap(colours) + + zz = np.zeros_like(self.xx) + zz[:] = z # self.bounding_box[1,2] + pts = np.vstack([self.xx.flatten(), self.yy.flatten(), zz.flatten()]) + if self.model is None: + logger.error("Mapview needs a model assigned to plot model on map") + return + vals = self.model.evaluate_model(pts.T, scale=True) + return self.ax.imshow( + vals.reshape(self.nsteps).T, + extent=[ + self.bounding_box[0, 0], + self.bounding_box[1, 0], + self.bounding_box[0, 1], + self.bounding_box[1, 1], + ], + origin="lower", + cmap=cmap, + ) + + def add_fault_displacements(self, z=0, cmap="rainbow"): + + zz = np.zeros_like(self.xx) + zz[:] = z # self.bounding_box[1,2] + pts = np.vstack([self.xx.flatten(), self.yy.flatten(), zz.flatten()]) + if self.model is None: + logger.error("Mapview needs a model assigned to plot model on map") + return + vals = self.model.evaluate_fault_displacements(pts.T, scale=True) + return self.ax.imshow( + vals.reshape(self.nsteps).T, + extent=[ + self.bounding_box[0, 0], + self.bounding_box[1, 0], + self.bounding_box[0, 1], + self.bounding_box[1, 1], + ], + origin="lower", + cmap=cmap, + ) + + def add_faults(self, **kwargs): + for f in self.model.features: + if f.type == FeatureType.FAULT: + # create a function to return true if displacement > 0 + def mask(x): + val = f.displacementfeature.evaluate_value(x) + val[np.isnan(val)] = 0 + maskv = np.zeros(val.shape).astype(bool) + maskv[np.abs(val) > 0.001] = 1 + return maskv + + self.add_contour(f, 0, mask=mask, **kwargs) diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_3d_viewer.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_3d_viewer.py new file mode 100644 index 000000000..093abcf7d --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_3d_viewer.py @@ -0,0 +1,708 @@ +import pyvista as pv +import numpy as np +import re + +from LoopStructural.datatypes import VectorPoints, ValuePoints +from LoopStructural.modelling.features import BaseFeature, StructuralFrame + +from LoopStructural.modelling.features.fault import FaultSegment +from LoopStructural.datatypes import BoundingBox +from LoopStructural import GeologicalModel +from LoopStructural.utils import getLogger +from typing import Callable, Union, Optional, List + +logger = getLogger(__name__) + + +class Loop3DView(pv.Plotter): + def __init__(self, model=None, background='white', *args, **kwargs): + """Loop3DView is a subclass of pyvista. Plotter that is designed to + interface with the LoopStructural geological modelling package. + + Parameters + ---------- + model : GeologicalModel, optional + A loopstructural model used as reference for some methods, by default None + background : str, optional + colour for the background, by default 'white' + """ + if 'shape' in kwargs: + logger.warning('shape argument is not used in Loop3DView') + kwargs.pop('shape') + super().__init__(*args, **kwargs) + self.set_background(background) + self.model = model + self.objects = {} + + def subplot(self, *args, **kwargs): + logger.warning('subplot is not supported in Loop3DView') + return self + + def add_mesh(self, *args, **kwargs): + if 'name' not in kwargs or kwargs['name'] is None: + name = 'unnamed_object' + kwargs['name'] = name + logger.warning( + f'No name provided, using {name}. Pass name argument to add_mesh to remove this error' + ) + kwargs['name'] = kwargs['name'].replace(' ', '_') + kwargs['name'] = re.sub(r'[^a-zA-Z0-9_$]', '_', kwargs['name']) + if kwargs['name'][0].isdigit(): + kwargs['name'] = 'ls_' + kwargs['name'] + if kwargs['name'][0] == '_': + kwargs['name'] = 'ls' + kwargs['name'] + kwargs['name'] = self.increment_name(kwargs['name']) + if '__opacity' in kwargs['name']: + raise ValueError('Cannot use __opacity in name') + if '__visibility' in kwargs['name']: + raise ValueError('Cannot use __visibility in name') + if '__control_visibility' in kwargs['name']: + raise ValueError('Cannot use __control_visibility in name') + return super().add_mesh(*args, **kwargs) + + def increment_name(self, name): + parts = name.split('_') + if len(parts) == 1: + name = name + '_1' + while name in self.actors: + parts = name.split('_') + try: + parts[-1] = str(int(parts[-1]) + 1) + except ValueError: + parts.append('1') + name = '_'.join(parts) + return name + + def _check_model(self, model: GeologicalModel) -> GeologicalModel: + """helper method to assign a geological model""" + if model is None: + model = self.model + if model is None: + raise ValueError("No model provided") + return model + + def _get_vector_scale(self, scale: Optional[Union[float, int]]) -> float: + autoscale = 1.0 + if self.model is not None: + # automatically scale vector data to be 5% of the bounding box length + autoscale = self.model.bounding_box.length.max() * 0.05 + if scale is None: + scale = autoscale + else: + scale = scale * autoscale + if scale > 10 * autoscale: + logger.warning( + "Vector scale magnitude greater than half of the model bounding box length, is this correct?" + ) + + return scale + + def plot_surface( + self, + geological_feature: BaseFeature, + value: Optional[Union[float, int]] = None, + paint_with: Optional[BaseFeature] = None, + colour: Optional[str] = "red", + cmap: Optional[str] = None, + opacity: Optional[float] = None, + vmin: Optional[float] = None, + vmax: Optional[float] = None, + pyvista_kwargs: dict = {}, + show_scalar_bar: bool = False, + slicer: bool = False, + name: Optional[str] = None, + bounding_box: Optional[BoundingBox] = None, + ): + """Add an isosurface of a geological feature to the model + + Parameters + ---------- + geological_feature : BaseFeature + The geological feature to plot + value : Optional[Union[float, int, List[float]]], optional + isosurface value, or list of values, by default average value of feature + paint_with : Optional[BaseFeature], optional + Paint the surface with the value of another geological feature, by default None + colour : Optional[str], optional + colour of the surface, by default "red" + cmap : Optional[str], optional + matplotlib colourmap, by default None + opacity : Optional[float], optional + opacity of the surface, by default None + vmin : Optional[float], optional + minimum value of the colourmap, by default None + vmax : Optional[float], optional + maximum value of the colourmap, by default None + pyvista_kwargs : dict, optional + other parameters passed to Plotter.add_mesh, by default {} + name : Optional[str], optional + name of the object, by default None + slicer : bool, optional + If an interactive plane slicing tool should be added, by default False + show_scalar_bar : bool, optional + Whether to show the scalar bar, by default False + """ + + if name is None: + name = geological_feature.name + '_surfaces' + name = self.increment_name(name) # , 'surface') + + surfaces = geological_feature.surfaces(value, bounding_box=bounding_box) + meshes = [] + for surface in surfaces: + s = surface.vtk() + if paint_with is not None: + clim = [paint_with.min(), paint_with.max()] + if vmin is not None: + clim[0] = vmin + if vmax is not None: + clim[1] = vmax + pyvista_kwargs["clim"] = clim + pts = np.copy(surface.vertices) + if self.model is not None: + pts = self.model.scale(pts) + scalars = paint_with(pts) + s["values"] = scalars + s.set_active_scalars("values") + colour = None + meshes.append(s) + mesh = pv.MultiBlock(meshes).combine() + actor = None + try: + + if slicer: + actor = self.add_mesh_clip_plane( + mesh, + color=colour, + cmap=cmap, + opacity=opacity, + name=name, + **pyvista_kwargs, + ) + else: + actor = self.add_mesh( + mesh, + color=colour, + cmap=cmap, + opacity=opacity, + name=name, + **pyvista_kwargs, + ) + + except ValueError: + logger.warning("No surfaces to plot") + if paint_with is not None and not show_scalar_bar: + self.remove_scalar_bar('values') + return actor + + def plot_scalar_field( + self, + geological_feature: BaseFeature, + cmap: str = "viridis", + vmin: Optional[float] = None, + vmax: Optional[float] = None, + opacity: Optional[float] = None, + pyvista_kwargs: dict = {}, + show_scalar_bar: bool = False, + slicer: bool = False, + name: Optional[str] = None, + bounding_box: Optional[BoundingBox] = None, + ): + """Plot a volume with the scalar field as the property + calls feature.scalar_field() to get the scalar field and + then pyvista add_mesh(feature.scalar_field().vtk()) + + Parameters + ---------- + geological_feature : BaseFeature + The geological feature to plot the scalar field of + cmap : str, optional + matplotlib colourmap to use, by default "viridis" + vmin : Optional[float], optional + minimum value for cmap, by default None + vmax : Optional[float], optional + max value for cmap, by default None + opacity : Optional[float], optional + opacity of the object, by default None + pyvista_kwargs : dict, optional + additional kwargs sent to add_mesh, by default {} + show_scalar_bar : bool, optional + whether to show or hide the scalar bar, by default False + slicer : bool, optional + whether to plot using a plane slicer widget, by default False + name : Optional[str], optional + name for the object to appear in the object list, by default None + + Returns + ------- + pv.Actor + a reference to the actor that is added to the mesh + """ + + if name is None: + name = geological_feature.name + '_scalar_field' + name = self.increment_name(name) # , 'scalar_field') + + volume = geological_feature.scalar_field(bounding_box=bounding_box).vtk() + if vmin is not None: + pyvista_kwargs["clim"][0] = vmin + if vmax is not None: + pyvista_kwargs["clim"][1] = vmax + if slicer: + actor = self.add_mesh_clip_plane( + volume, cmap=cmap, opacity=opacity, name=name, **pyvista_kwargs + ) + else: + actor = self.add_mesh(volume, cmap=cmap, opacity=opacity, name=name, **pyvista_kwargs) + if not show_scalar_bar: + self.remove_scalar_bar(geological_feature.name) + return actor + + def plot_block_model( + self, + cmap=None, + model=None, + pyvista_kwargs={}, + show_scalar_bar: bool = False, + slicer: bool = False, + threshold: Optional[Union[float, List[float]]] = None, + name: Optional[str] = None, + ): + """Plot a voxel model where the stratigraphic id is the active scalar. + It will use the colours defined in the stratigraphic column of the model + unless a cmap is provided. + Min/max range of cmap are defined by the min/max values of the stratigraphic ids or if + clim is provided in pyvista_kwargs + + Parameters + ---------- + cmap : str, optional + matplotlib cmap string, by default None + model : GeologicalModel, optional + the model to pass if it is not the active geologicalmodel, by default None + pyvista_kwargs : dict, optional + additional arguments to be passed to pyvista add_mesh, by default {} + show_scalar_bar : bool, optional + whether show/hide the scalar bar, by default False + slicer : bool, optional + If an interactive plane slicing tool should be added, by default False + threshold : Optional[Union[float, List[float]]], optional + Whether to threshold values of the stratigraphy. Uses same syntax as pyvista threshold., by default None + """ + model = self._check_model(model) + if name is None: + name = 'block_model' + name = self.increment_name(name) # , 'block_model') + block, codes = model.get_block_model() + block = block.vtk() + block.set_active_scalars('stratigraphy') + actor = None + if cmap is None: + cmap = self._build_stratigraphic_cmap(model) + if "clim" not in pyvista_kwargs: + pyvista_kwargs["clim"] = (np.min(block['stratigraphy']), np.max(block['stratigraphy'])) + if threshold is not None: + if isinstance(threshold, float): + block = block.threshold(threshold) + elif isinstance(threshold, (list, tuple, np.ndarray)) and len(threshold) == 2: + block = block.threshold((threshold[0], threshold[1])) + if slicer: + actor = self.add_mesh_clip_plane(block, cmap=cmap, name=name, **pyvista_kwargs) + else: + actor = self.add_mesh(block, cmap=cmap, name=name, **pyvista_kwargs) + + if not show_scalar_bar: + self.remove_scalar_bar('stratigraphy') + return actor + + def plot_fault_displacements( + self, + fault_list: Optional[List[FaultSegment]] = None, + bounding_box: Optional[BoundingBox] = None, + model=None, + cmap="rainbow", + pyvista_kwargs={}, + show_scalar_bar: bool = False, + name: Optional[str] = None, + ): + """Plot the dispalcement magnitude for faults in the model + on a voxel block + + Parameters + ---------- + fault_list : _type_, optional + list of faults to plot the model, by default None + bounding_box : _type_, optional + _description_, by default None + model : _type_, optional + _description_, by default None + cmap : str, optional + _description_, by default "rainbow" + pyvista_kwargs : dict, optional + _description_, by default {} + show_scalar_bar : bool, optional + _description_, by default False + """ + if name is None: + name = 'fault_displacement' + name = self.increment_name(name) # , 'fault_displacement_map') + if fault_list is None: + model = self._check_model(model) + fault_list = model.faults + if bounding_box is None: + model = self._check_model(model) + bounding_box = model.bounding_box + pts = bounding_box.regular_grid() + displacement_value = np.zeros(pts.shape[0]) + for f in fault_list: + disp = f.displacementfeature.evaluate_value(bounding_box.vtk().points) + displacement_value[~np.isnan(disp)] += disp[~np.isnan(disp)] + volume = bounding_box.vtk() + volume['displacement'] = displacement_value + actor = self.add_mesh(volume, cmap=cmap, **pyvista_kwargs) + if not show_scalar_bar: + self.remove_scalar_bar('displacement') + return actor + + def plot_model_surfaces( + self, + strati: bool = True, + faults: bool = True, + cmap: Optional[str] = None, + model: Optional[GeologicalModel] = None, + fault_colour: str = "black", + pyvista_kwargs: dict = {}, + show_scalar_bar: bool = False, + name: Optional[str] = None, + ): + """Plot the surfaces of the model + + Parameters + ---------- + strati : bool, optional + should stratigraphy surfaces be plotted, by default True + faults : bool, optional + should faults be plotted, by default True + cmap : Optional[str], optional + What cmap to use for the stratigraphy ids, by default None + model : Optional[GeologicalModel], optional + a GeologicalModel, if not provided will use self.model, by default None + fault_colour : str, optional + colour for the fault surfaces, by default "black" + pyvista_kwargs : dict, optional + Additional kwargs to send to add_mesh, by default {} + show_scalar_bar : bool, optional + whether to add the scalar bar, by default False + name : Optional[str], optional + name to add objects to object list with, by default None + + Returns + ------- + pv.Actor + The actor that is added to the scene + """ + model = self._check_model(model) + + actors = [] + if strati: + strati_surfaces = [] + surfaces = model.get_stratigraphic_surfaces() + if cmap is None: + cmap = model.stratigraphic_column.cmap().colors + for s in surfaces: + strati_surfaces.append(s.vtk()) + if name is None: + object_name = 'model_surfaces' + else: + object_name = f'{name}_model_surfaces' + object_name = self.increment_name(object_name) # , 'model_surfaces') + actors.append( + self.add_mesh( + pv.MultiBlock(strati_surfaces).combine(), + cmap=cmap, + name=object_name, + **pyvista_kwargs, + ) + ) + if not show_scalar_bar: + self.remove_scalar_bar() + if faults: + fault_list = model.get_fault_surfaces() + for f in fault_list: + if name is None: + object_name = f'{f.name}_surface' + if name is not None: + object_name = f'{name}_{f.name}_surface' + object_name = self.increment_name(object_name) # , 'fault_surfaces') + actors.append( + self.add_mesh(f.vtk(), color=fault_colour, name=object_name, **pyvista_kwargs) + ) + return actors + + def plot_vector_field( + self, + geological_feature: BaseFeature, + scale: Optional[float] = None, + name: Optional[str] = None, + geom='arrow', + scalars: Optional[np.ndarray] = None, + normalise: bool = False, + scale_function: Optional[Callable[[np.ndarray], np.ndarray]] = None, + pyvista_kwargs: dict = {}, + bounding_box: Optional[BoundingBox] = None, + ) -> pv.Actor: + """Plot a vector field + + Parameters + ---------- + geological_feature : BaseFeature + Geological feature to plot the vector field of + scale : float, optional + magnitude scale for the glyphs, by default 1.0 + name : Optional[str], optional + name for the viewer object list, by default None + pyvista_kwargs : dict, optional + additional kwargs to pass to add_mesh, by default {} + + Returns + ------- + pv.Actor + actor that is added to the scene + """ + if name is None: + name = geological_feature.name + '_vector_field' + name = self.increment_name(name) # , 'vector_field') + vectorfield = geological_feature.vector_field(bounding_box=bounding_box) + scale = self._get_vector_scale(scale) + return self.add_mesh( + vectorfield.vtk( + scale=scale, + geom=geom, + normalise=normalise, + scalars=scalars, + scale_function=scale_function, + ), + name=name, + **pyvista_kwargs, + ) + + def plot_data( + self, + feature: Union[BaseFeature, StructuralFrame], + value: bool = True, + vector: bool = True, + scale: Optional[Union[float, int]] = None, + geom: str = "arrow", + name: Optional[str] = None, + scalars: Optional[np.ndarray] = None, + normalise: bool = True, + pyvista_kwargs: dict = {}, + ) -> List[pv.Actor]: + """Add the data associated with a feature to the plotter + + Parameters + ---------- + feature : Union[BaseFeature, StructuralFrame] + feature to add data from + value : bool, optional + whether to add value data, by default True + vector : bool, optional + whether to plot vector data, by default True + scale : Union[float, int], optional + vector scale, by default 1 + geom : str, optional + vector glyph, by default "arrow" + name : Optional[str], optional + name to use in object list, by default None + normalise: bool, optional + normalise the vectors to be unit norm, by default True + pyvista_kwargs : dict, optional + additional kwargs to pass to pyvista add_mesh, by default {} + + Returns + ------- + List[pv.Actor] + list of actors added to the pv plotter + + Notes + ------ + When plotting a vector the bounding box is used to scale the vectors. By default + the length of the arrows will be 5% of the bounding box. The scale parameter is a + multiplier for this value. If you sent normalise to False the vectors will not be normalised + + """ + if issubclass(type(feature), BaseFeature): + feature = [feature] + logger.info(f"Scale vectors by {scale}") + scale = self._get_vector_scale(scale) + logger.info(f"Vector scale is {scale}") + actors = [] + bb = self.model.bounding_box if self.model is not None else None + for f in feature: + for d in f.get_data(): + if isinstance(d, ValuePoints): + if value: + if name is None: + object_name = d.name + '_values' + else: + object_name = f'{d.name}_values_{name}' + object_name = self.increment_name(object_name) # , 'values') + actors.append( + self.add_mesh( + d.vtk(scalars=scalars), name=object_name, **pyvista_kwargs + ) + ) + if isinstance(d, VectorPoints): + if vector: + if name is None: + object_name = d.name + '_vectors' + else: + object_name = f'{d.name}_vectors_{name}' + object_name = self.increment_name(object_name) # , 'vectors') + actors.append( + self.add_mesh( + d.vtk( + geom=geom, + scale=scale, + scalars=scalars, + bb=bb, + tolerance=None, + normalise=normalise, + ), + name=name, + **pyvista_kwargs, + ) + ) + return actors + + def plot_fold(self, folded_feature: BaseFeature, pyvista_kwargs={}): + + # folded_feature. + pass + + def plot_fault( + self, + fault: FaultSegment, + surface: bool = True, + slip_vector: bool = True, + displacement_scale_vector: bool = True, + fault_volume: bool = True, + vector_scale: Optional[Union[float, int]] = None, + name: Optional[str] = None, + geom: str = "arrow", + pyvista_kwargs: dict = {}, + bounding_box: Optional[BoundingBox] = None, + ) -> List[pv.Actor]: + """Plot a fault including the surface, slip vector and displacement volume + + Parameters + ---------- + fault : FaultSegment + the fault to plot + surface : bool, optional + flag for the 0.0 surface, by default True + slip_vector : bool, optional + flag for scaled vector field, by default True + displacement_scale_vector : bool, optional + _description_, by default True + fault_volume : bool, optional + fault displacement scalar field, by default True + vector_scale : Union[float, int], optional + scale factor for vectors, by default 200 + name : Optional[str], optional + name of the object for pyvista, by default None + pyvista_kwargs : dict, optional + additional kwargs for the pyvista plotter, by default {} + + Returns + ------- + List[pv.Actor] + list of actors added to the plot + """ + actors = [] + if surface: + if name is None: + surface_name = fault.name + '_surface' + else: + surface_name = f'{fault.name}_surface_{name}' + surface_name = self.increment_name(surface_name) + surf = fault.surfaces([0], bounding_box=bounding_box)[0] + actors.append(self.add_mesh(surf.vtk(), name=surface_name, **pyvista_kwargs)) + if slip_vector: + if name is None: + vector_name = fault.name + '_vector' + else: + vector_name = f'{fault.name}_vector_{name}' + vector_name = self.increment_name(vector_name) + + vectorfield = fault.vector_field(bounding_box=bounding_box) + vector_scale = self._get_vector_scale(vector_scale) + actors.append( + self.add_mesh( + vectorfield.vtk(scale=vector_scale, normalise=False), + name=vector_name, + **pyvista_kwargs, + ) + ) + if fault_volume: + if name is None: + volume_name = fault.name + '_volume' + else: + volume_name = f'{fault.name}_volume_{name}' + + volume = fault.displacementfeature.scalar_field(bounding_box=bounding_box) + + volume = volume.vtk().threshold([-1.0, 1.0]) + if geom == "arrow": + geom = pv.Arrow() + elif geom == "disc": + geom = pv.Disc() + geom = geom.rotate_y(90) + else: + raise ValueError(f"Unknown glyph type {geom}") + actors.append(self.add_mesh(volume, name=volume_name, **pyvista_kwargs)) + if len(actors) == 0: + logger.warning(f"Nothing added to plot for {fault.name}") + return actors + + def plot_fault_ellipsoid( + self, fault: FaultSegment, name: Optional[str] = None, pyvista_kwargs: dict = {} + ) -> pv.Actor: + """Plot the fault ellipsoid + + Parameters + ---------- + fault : FaultSegment + the fault to plot + name : Optional[str], optional + name of the object for pyvista, by default None + pyvista_kwargs : dict, optional + additional kwargs for the pyvista plotter, by default {} + + Returns + ------- + pv.Actor + actor added to the plot + """ + if name is None: + name = fault.name + '_ellipsoid' + name = self.increment_name(name) + ellipsoid = fault.fault_ellipsoid() + return self.add_mesh(ellipsoid, name=name, **pyvista_kwargs) + + def rotate(self, angles: np.ndarray): + """Rotate the camera by the given angles + order is roll, azimuth, elevation as defined by + pyvista + + Parameters + ---------- + angles : np.ndarray + roll, azimuth, elevation + """ + self.camera.roll += angles[0] + self.camera.azimuth += angles[1] + self.camera.elevation += angles[2] + + def display(self): + self.show(interactive=False) diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/__init__.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/__init__.py new file mode 100644 index 000000000..cb83b3b04 --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/__init__.py @@ -0,0 +1,13 @@ +from LoopStructural.utils import getLogger + +from ._3d_viewer import Loop3DView +from ._rotation_angle import RotationAnglePlotter +from ._2d_viewer import Loop2DView +from ._stratigraphic_column import StratigraphicColumnView + +logger = getLogger(__name__) + +try: + from ._register_loop_ui import * # this replaces default pyvista trame ui +except ImportError: + logger.warning("Could not import trame ui") diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_register_loop_ui.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_register_loop_ui.py new file mode 100644 index 000000000..8ec129b4a --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_register_loop_ui.py @@ -0,0 +1,5 @@ +import pyvista +from .trame import initialize + + +pyvista.trame.jupyter.initialize = initialize diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_rotation_angle.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_rotation_angle.py new file mode 100644 index 000000000..a02cd5dda --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_rotation_angle.py @@ -0,0 +1,88 @@ +import matplotlib.pyplot as plt +import numpy as np + +from LoopStructural.utils import getLogger + +logger = getLogger(__name__) + + +class RotationAnglePlotter: + def __init__(self, feature, axis=True): + """ """ + self.fig, self.ax = plt.subplots(2, 2, figsize=(30, 15)) + self.ax[0][0].set_ylim(-90, 90) + self.ax[1][0].set_ylim(-90, 90) + self.feature = feature + self.feature.builder.up_to_date() + + def plot(self, x, y, ix, iy, symb, **kwargs): + """ + + Parameters + ---------- + x : np.array + vector of x + y + ix + iy + symb + + Returns + ------- + + """ + return self.ax[iy][ix].plot(x, y, symb, **kwargs) + + def default_titles(self): + self.ax[0][0].set_title("Fold Axis S-Plot") + self.ax[0][1].set_title("Fold Axis S-Variogram") + self.ax[1][0].set_title("Fold Limb S-Plot") + self.ax[1][1].set_title("Fold Limb S-Variogram") + + self.ax[1][1].set_xlabel("Variogram Steps") + self.ax[1][1].set_ylabel("Fold Limb S-Variogram") + self.ax[1][0].set_ylabel("Fold Limb Rotation Angle") + self.ax[1][0].set_xlabel("Fold Frame Axial Surface Field") + + self.ax[0][1].set_xlabel("Variogram Steps") + self.ax[0][1].set_ylabel("Fold Axis S-Variogram") + self.ax[0][0].set_ylabel("Fold Axis Rotation Angle") + self.ax[0][0].set_xlabel("Fold Frame Axis Direction Field") + + def add_fold_limb_data(self, symb="bo", **kwargs): + fold_frame = self.feature.builder.fold.fold_limb_rotation.fold_frame_coordinate + rotation = self.feature.fold.fold_limb_rotation.rotation_angle + return self.plot(fold_frame, rotation, 0, 1, symb, **kwargs) + + def add_fold_limb_curve(self, symb="r-", **kwargs): + x = np.linspace( + self.feature.fold.foldframe[0].min(), + self.feature.fold.foldframe[0].max(), + 100, + ) + return self.plot(x, self.feature.builder.fold.fold_limb_rotation(x), 0, 1, symb, **kwargs) + + def add_axis_svariogram(self, symb="bo", **kwargs): + svariogram = self.feature.builder.fold.fold_axis_rotation.svario + if svariogram: + svariogram.calc_semivariogram() + return self.plot(svariogram.lags, svariogram.variogram, 1, 0, symb, **kwargs) + + def add_limb_svariogram(self, symb="bo", **kwargs): + svariogram = self.feature.builder.fold.fold_limb_rotation.svario + if svariogram: + svariogram.calc_semivariogram() + return self.plot(svariogram.lags, svariogram.variogram, 1, 1, symb, **kwargs) + + def add_fold_axis_data(self, symb="bo", **kwargs): + fold_frame = self.feature.builder.fold.fold_axis_rotation.fold_frame_coordinate + rotation = self.feature.builder.fold.fold_axis_rotation.rotation_angle + return self.plot(fold_frame, rotation, 0, 0, symb, **kwargs) + + def add_fold_axis_curve(self, symb="r-", **kwargs): + x = np.linspace( + self.feature.builder.fold.foldframe[1].min(), + self.feature.builder.fold.foldframe[1].max(), + 100, + ) + return self.plot(x, self.feature.builder.fold.fold_axis_rotation(x), 0, 0, symb, **kwargs) diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_stratigraphic_column.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_stratigraphic_column.py new file mode 100644 index 000000000..9dca3483c --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/_stratigraphic_column.py @@ -0,0 +1,81 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import cm +from matplotlib.patches import Polygon +from matplotlib.collections import PatchCollection +from LoopStructural.utils import rng + + +class StratigraphicColumnView: + def __init__(self, model, ax=None, cmap=None, labels=None): + self.model = model + self.ax = ax + self.cmap = cmap + self.labels = labels + + def plot(self): + n_units = 0 # count how many discrete colours (number of stratigraphic units) + xmin = 0 + ymin = 0 + ymax = 1 + xmax = 1 + fig = None + if self.ax is None: + fig, self.ax = plt.subplots(figsize=(2, 10)) + patches = [] # stores the individual stratigraphic unit polygons + + total_height = 0 + prev_coords = [0, 0] + + # iterate through groups, skipping faults + + for g in reversed(self.model.stratigraphic_column.get_groups()): + for u in g.units: + n_units += 1 + + ymax = total_height + ymin = ymax - (u.thickness) + + if not np.isfinite(ymin): + ymin = prev_coords[1] - (prev_coords[1] - prev_coords[0]) * (1 + rng.random()) + + total_height = ymin + + prev_coords = (ymin, ymax) + + polygon_points = np.array([[xmin, ymin], [xmax, ymin], [xmax, ymax], [xmin, ymax]]) + patches.append(Polygon(polygon_points)) + xy = (0, ymin + (ymax - ymin) / 2) + if self.labels: + self.ax.annotate(self.labels[u], xy) + else: + self.ax.annotate(u.name, xy) + + if self.cmap is None: + import matplotlib.colors as colors + + colours = [] + boundaries = [] + data = [] + for g in self.model.stratigraphic_column.get_groups(): + if g == "faults": + continue + for u in g.units: + data.append((u.id, u.colour)) + colours.append(u.colour) + boundaries.append(u.id) # print(u,v) + cmap = colors.ListedColormap(colours) + else: + cmap = cm.get_cmap(self.cmap, n_units - 1) + p = PatchCollection(patches, cmap=cmap) + + colors = np.arange(len(patches)) + p.set_array(np.array(colors)) + + self.ax.add_collection(p) + + self.ax.set_ylim(total_height - (total_height - prev_coords[0]) * 0.1, 0) + + self.ax.axis("off") + + return fig diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/api.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/api.py new file mode 100644 index 000000000..035859c72 --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/api.py @@ -0,0 +1,39 @@ +from ._3d_viewer import Loop3DView + + +def plot_block_model(model, filename=None, **kwargs): + """ + Plot the model using pyvista + + Parameters + ---------- + model : LoopStructuralModel + The model to plot + kwargs : dict + Keyword arguments to pass to the plot + """ + + p = Loop3DView(model) + p.plot_block_model(**kwargs) + if filename is not None: + p.screenshot(filename) + p.show() + + +def plot_surface(model, geological_feature, **kwargs): + """ + Plot a surface using pyvista + + Parameters + ---------- + model : LoopStructuralModel + The model to plot + geological_feature : BaseFeature + The feature to plot + kwargs : dict + Keyword arguments to pass to the plot + """ + + p = Loop3DView(model) + p.plot_surface(geological_feature, **kwargs) + return p diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/trame/__init__.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/trame/__init__.py new file mode 100644 index 000000000..8cf564742 --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/trame/__init__.py @@ -0,0 +1,40 @@ +from pyvista.trame.ui import UI_TITLE +from pyvista.trame.ui import get_viewer +from .ui.vuetify3 import LoopViewer as Viewer +from .. import Loop3DView + + +def initialize( + server, + plotter, + mode=None, + default_server_rendering=True, + collapse_menu=False, + **kwargs, +): # numpydoc ignore=PR01,RT01 + """Generate the UI for a given plotter.""" + state = server.state + state.trame__title = UI_TITLE + if issubclass(type(plotter), Loop3DView): + # only use the loopviewer if the plotter is a Loop3DView + viewer = Viewer(plotter, server=server) + + else: + # if pyvista use trame.ui.Viewer + viewer = get_viewer( + plotter, + server=server, + suppress_rendering=mode == "client", + ) + with viewer.make_layout(server, template_name=plotter._id_name) as layout: + viewer.layout = layout + viewer.ui( + mode=mode, + default_server_rendering=default_server_rendering, + collapse_menu=collapse_menu, + **kwargs, + ) + if issubclass(type(plotter), Loop3DView): + viewer.object_menu() + + return viewer diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/trame/ui/vuetify3.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/trame/ui/vuetify3.py new file mode 100644 index 000000000..c9511ad35 --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/trame/ui/vuetify3.py @@ -0,0 +1,109 @@ +# ruff: noqa: D102 +"""PyVista Trame Viewer class for a Vue 3 client. + +This class, derived from `pyvista.trame.ui.base_viewer`, +is intended for use with a trame application where the client type is "vue3". +Therefore, the `ui` method implemented by this class utilizes the API of Vuetify 3. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from trame.ui.vuetify3 import SinglePageWithDrawerLayout +from trame.widgets import vuetify3 as vuetify + +import pyvista +from pyvista.trame.ui.vuetify3 import Viewer + + +if TYPE_CHECKING: # pragma: no cover + from trame_client.ui.core import AbstractLayout + + +class LoopViewer(Viewer): + def __init__(self, *args, **kwargs): + """Overwrite the pyvista trame layout to use a singlepage layout + and add an object visibility menu to the drawer + """ + super().__init__(*args, **kwargs) + + def make_layout(self, *args, **kwargs) -> AbstractLayout: + + return SinglePageWithDrawerLayout(*args, **kwargs) + + def ui(self, *args, **kwargs): + with self.layout as layout: + layout.title.set_text("LoopStructural Viewer") + with self.layout.content: + + return super().ui(*args, **kwargs) + + def toggle_visibility(self, **kwargs): + """Toggle the visibility of an object in the plotter. + this is a slot called by the state change, the kwargs are the current state + so we need to check the keys and update accordingly + """ + for k in kwargs.keys(): + object_name = k.split("__visibility")[0] + if object_name in self.plotter.actors: + self.plotter.actors[object_name].visibility = kwargs[k] + self.update() + # self.actors[k].visibility = not self.actors[k].visibility + + def set_opacity(self, **kwargs): + """Set the opacity of an object in the plotter. + this is a slot called by the state change, the kwargs are the current state + so we need to check the keys and update accordingly + """ + for k in kwargs.keys(): + object_name = k.split("__opacity")[0] + if object_name in self.plotter.actors: + self.plotter.actors[object_name].prop.opacity = kwargs[k] + self.update() + # self.actors[k].visibility = not self.actors[k].visibility + + def object_menu(self): + with self.layout.drawer as drawer: + with vuetify.VCard(): + + for k, a in self.plotter.actors.items(): + if type(a) is not pyvista.plotting.actor.Actor: + continue + drawer.server.state[f"{k}__visibility"] = True + drawer.server.state[f"{k}__control_visibility"] = False + drawer.server.state[f"{k}__opacity"] = a.prop.opacity + drawer.server.state.change(f"{k}__visibility")(self.toggle_visibility) + drawer.server.state.change(f"{k}__opacity")(self.set_opacity) + with vuetify.VRow( + classes='pa-0 ma-0 align-center fill-height', + style='flex-wrap: nowrap', + ): + with vuetify.VCol(): + + vuetify.VCheckbox( + label=k, + classes="ma-0 pa-0", + v_model=(f"{k}__visibility"), + # click=(self.toggle_visibility("test")), + ) + with vuetify.VCol(): + vuetify.VBtn( + icon='mdi-dots-horizontal', + click=(f'{k}__control_visibility=!{k}__control_visibility'), + # click=(self.toggle_visibility("test")), + ) + with vuetify.VCard(classes="ma-0 pa-0", v_show=f'{k}__control_visibility'): + + with vuetify.VRow( + classes="ma-0 pa-0 d-flex align-center", + ): + + vuetify.VSlider( + v_model=(f"{k}__opacity"), + label="Opacity", + min=0, + max=1, + step=0.1, + thumb_label=True, + ) diff --git a/packages/loopstructural_visualisation/src/loopstructuralvisualisation/version.py b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/version.py new file mode 100644 index 000000000..86205cbac --- /dev/null +++ b/packages/loopstructural_visualisation/src/loopstructuralvisualisation/version.py @@ -0,0 +1 @@ +__version__ = "0.1.17" diff --git a/pyproject.toml b/pyproject.toml index e2bd24cba..1e87b78d3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -90,6 +90,8 @@ members = ["packages/*"] [tool.uv.sources] loop-common = { workspace = true } loop-interpolation = { workspace = true } +loopstructuralvisualisation = { workspace = true } +LoopStructural = { workspace = true } [tool.pytest.ini_options] addopts = "--import-mode=importlib"