# Copyright (c) 2024 Mira Geoscience Ltd.
#
# This file is part of geoh5py.
#
# geoh5py is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# geoh5py 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 Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with geoh5py. If not, see <https://www.gnu.org/licenses/>.
from __future__ import annotations
import uuid
import numpy as np
from .grid_object import GridObject
from .object_base import ObjectType
[docs]
class DrapeModel(GridObject):
"""
Drape (curtain) model object made up of layers and prisms.
"""
__TYPE_UID = uuid.UUID("{C94968EA-CF7D-11EB-B8BC-0242AC130003}")
def __init__(self, object_type: ObjectType, **kwargs):
self._layers: np.ndarray | None = None
self._prisms: np.ndarray | None = None
super().__init__(object_type, **kwargs)
object_type.workspace._register_object(self)
[docs]
@classmethod
def default_type_uid(cls) -> uuid.UUID:
return cls.__TYPE_UID
@property
def centroids(self):
"""
:obj:`numpy.array` of :obj:`float`,
shape (:obj:`~geoh5py.objects.drape_model.Drapemodel.n_cells`, 3):
Cell center locations in world coordinates.
.. code-block:: python
centroids = [
[x_1, y_1, z_1],
...,
[x_N, y_N, z_N]
]
"""
if getattr(self, "_centroids", None) is None:
if self.layers is None:
raise AttributeError(
"Attribute 'layers' must be defined before accessing 'centroids'."
)
if self.prisms is None:
raise AttributeError(
"Attribute 'prisms' must be defined before accessing 'centroids'."
)
self._centroids = np.vstack(
[
np.ones((int(val), 3)) * self.prisms[ii, :3]
for ii, val in enumerate(self.prisms[:, 4])
]
)
tops = np.hstack(
[
np.r_[
cells[2],
self.layers[int(cells[3]) : int(cells[3] + cells[4] - 1), 2],
]
for cells in self.prisms.tolist()
]
)
self._centroids[:, 2] = (tops + self.layers[:, 2]) / 2.0
return self._centroids
@property
def layers(self) -> np.ndarray | None:
"""
:obj:`numpy.array`, shape(*, 3): Layers in the drape model with columns: X
(prism index), K (depth index), elevation (cell bottom)).
shape(*, 3) organized into blocks representing each prism in the model.
.. code-block:: python
layers = [
[x_1, k_1, z_11],
[x_1, k_2, z_12],
...
[x_1, k_N, z_1N],
.
.
.
[x_M, k_1, z_M1],
[x_M, k_2, z_M2],
...
[x_M, k_N, z_MM]
]
"""
if self._layers is None and self.on_file:
self._layers = self.workspace.fetch_array_attribute(self, "layers")
if self._layers is not None:
return np.asarray(self._layers.tolist())
return None
@layers.setter
def layers(self, xyz: np.ndarray):
if any(np.diff(np.unique(xyz[:, 0])) != 1):
msg = "Prism index (first column) must be monotonically increasing."
raise ValueError(msg)
if xyz.shape[1] != 3:
msg = f"Array of layers must be of shape (*, 3). Array of shape {xyz.shape} provided."
raise ValueError(msg)
self._layers = np.asarray(
np.core.records.fromarrays(
xyz.T.tolist(),
dtype=[("I", "<i4"), ("K", "<i4"), ("Bottom elevation", "<f8")],
)
)
self.workspace.update_attribute(self, "layers")
@property
def n_cells(self):
if self._prisms is not None:
return int(self._prisms["Layer count"].sum())
return None
@property
def prisms(self) -> np.ndarray | None:
"""
:obj:`numpy.array`, shape(*, 5) detailing the assembly of :obj:
`geoh5py.objects.drape_model.Drapemodel.layers` within the trace
of the drape model.
Columns: Easting, Northing, Elevation (top),
layer index (first), layer count.
.. code-block:: python
prisms = [
[e_1, n_1, z_1, l_1, c_1],
...,
[e_N, n_N, z_N, l_N, c_N]
]
"""
if self._prisms is None and self.on_file:
self._prisms = self.workspace.fetch_array_attribute(self, "prisms")
if self._prisms is not None:
return np.array(self._prisms.tolist())
return None
@prisms.setter
def prisms(self, xyz: np.ndarray):
assert (
xyz.shape[1] == 5
), f"Array of prisms must be of shape (*, 5). Array of shape {xyz.shape} provided."
self._prisms = np.asarray(
np.core.records.fromarrays(
xyz.T.tolist(),
dtype={
"names": [
"Top easting",
"Top northing",
"Top elevation",
"First layer",
"Layer count",
],
"formats": ["<f8", "<f8", "<f8", "<i4", "<i4"],
},
)
)
self.workspace.update_attribute(self, "prisms")