GEOPY-2714: Refactor mask_by_extent utilities for surfaces

geoh5py/objects/block_model.py

@@ -25,7 +25,6 @@
 import numpy as np
 
 from ..data import Data
-from ..shared.utils import xy_rotation_matrix
 from .grid_object import GridObject
 
 
@@ -92,18 +91,13 @@ def centroids(self) -> np.ndarray:
             ]
         """
         if getattr(self, "_centroids", None) is None:
-            cell_center_u = self.local_axis_centers("u")
-            cell_center_v = self.local_axis_centers("v")
-            cell_center_z = self.local_axis_centers("z")
-            angle = np.deg2rad(self.rotation)
-            rot = xy_rotation_matrix(angle)
             u_grid, v_grid, z_grid = np.meshgrid(
-                cell_center_u, cell_center_v, cell_center_z
+                self.local_axis_centers("u"),
+                self.local_axis_centers("v"),
+                self.local_axis_centers("z"),
             )
-            xyz = np.c_[np.ravel(u_grid), np.ravel(v_grid), np.ravel(z_grid)]
-            self._centroids = np.dot(rot, xyz.T).T
-
-            self._centroids += np.asarray(self._origin.tolist())[None, :]
+            uvw = np.c_[np.ravel(u_grid), np.ravel(v_grid), np.ravel(z_grid)]
+            self._centroids = self.uvw_to_xyz(uvw)
 
         return self._centroids
 
@@ -149,6 +143,19 @@ def shaped_data_values(
             (self.shape[2], self.shape[0], self.shape[1]), order="F"
         ).transpose(2, 1, 0)
 
+    @property
+    def span(self) -> np.ndarray:
+        """
+        Upper and lower limits along u, v and w directions.
+        """
+        return np.vstack(
+            [
+                [self.u_cell_delimiters.min(), self.u_cell_delimiters.max()],
+                [self.v_cell_delimiters.min(), self.v_cell_delimiters.max()],
+                [self.z_cell_delimiters.min(), self.z_cell_delimiters.max()],
+            ]
+        )
+
     @property
     def u_cell_delimiters(self) -> np.ndarray:
         """

geoh5py/objects/grid2d.py

@@ -126,17 +126,10 @@ def centroids(self) -> np.ndarray:
             ]
         """
         if getattr(self, "_centroids", None) is None:
-            rotation_matrix = xy_rotation_matrix(np.deg2rad(self.rotation))
-            dip_matrix = yz_rotation_matrix(np.deg2rad(self.dip))
-
             u_grid, v_grid = np.meshgrid(self.cell_center_u, self.cell_center_v)
-            xyz = np.c_[np.ravel(u_grid), np.ravel(v_grid), np.zeros(self.n_cells)]
+            uvw = np.c_[np.ravel(u_grid), np.ravel(v_grid), np.zeros(self.n_cells)]
 
-            xyz_dipped = dip_matrix @ xyz.T
-            centroids = (rotation_matrix @ xyz_dipped).T
-            centroids += np.asarray(self._origin.tolist())[None, :]
-
-            self._centroids = centroids
+            self._centroids = self.uvw_to_xyz(uvw)
 
         return self._centroids
 
@@ -266,6 +259,19 @@ def shaped_data_values(
 
         return values.reshape(self.v_count, self.u_count)
 
+    @property
+    def span(self) -> np.ndarray:
+        """
+        Upper and lower limits along u, v and w directions.
+        """
+        return np.vstack(
+            [
+                np.sort([0, self.u_cell_size * self.u_count]),
+                np.sort([0, self.v_cell_size * self.v_count]),
+                [0, 0],
+            ]
+        )
+
     @property
     def u_cell_size(self) -> float:
         """

geoh5py/objects/grid_object.py

@@ -28,6 +28,7 @@
 
 from geoh5py.data import Data, DataAssociationEnum
 
+from ..shared.utils import xy_rotation_matrix, yz_rotation_matrix
 from .object_base import ObjectBase
 
 
@@ -64,6 +65,15 @@ def centroids(self) -> np.ndarray:
         Cell center locations in world coordinates of shape (n_cells, 3).
         """
 
+    @property
+    def extent(self) -> np.ndarray:
+        """
+        Compute outer extent of mesh span in world coordinates.
+        """
+        u, v, w = np.meshgrid(self.span[0, :], self.span[1, :], self.span[2, :])
+        xyz = self.uvw_to_xyz(np.c_[u.ravel(), v.ravel(), w.ravel()])
+        return np.c_[xyz.min(axis=0), xyz.max(axis=0)].T
+
     @property
     def n_cells(self) -> int:
         """
@@ -135,6 +145,31 @@ def shape(self) -> np.ndarray:
         Cell center locations in world coordinates.
         """
 
+    @property
+    @abstractmethod
+    def span(self) -> np.ndarray:
+        """
+        Upper and lower limits along u, v and w directions.
+        """
+
+    def uvw_to_xyz(self, coordinates: np.ndarray) -> np.ndarray:
+        """
+        Apply rotation to the input coordinates.
+
+        :param coordinates: Array of coordinates along the axes, transformed to world coordinates.
+
+        :return:
+        """
+
+        rotation_matrix = xy_rotation_matrix(np.deg2rad(getattr(self, "rotation", 0)))
+        dip_matrix = yz_rotation_matrix(np.deg2rad(getattr(self, "dip", 0)))
+
+        xyz_dipped = dip_matrix @ coordinates.T
+        xyz = (rotation_matrix @ xyz_dipped).T
+        xyz += np.asarray(self._origin.tolist())[None, :]
+
+        return xyz
+
     def validate_cell_mask(self, cell_mask: np.ndarray | None):
         """
         Validate cell mask array, which is the same as validate_mask for grid objects.

geoh5py/objects/octree.py

@@ -209,6 +209,19 @@ def shape(self) -> tuple[np.int32, np.int32, np.int32]:
         """
         return self.u_count, self.v_count, self.w_count
 
+    @property
+    def span(self) -> np.ndarray:
+        """
+        Upper and lower limits along u, v and w directions.
+        """
+        return np.vstack(
+            [
+                np.sort([0, self.u_cell_size * self.u_count]),
+                np.sort([0, self.v_cell_size * self.v_count]),
+                np.sort([0, self.w_cell_size * self.w_count]),
+            ]
+        )
+
     @property
     def u_cell_size(self) -> float:
         """

geoh5py/objects/vp_model.py

@@ -318,18 +318,14 @@ def centroids(self) -> np.ndarray:
         """
         if getattr(self, "_centroids", None) is None:
             # Create a grid of u, v coordinates based on the cell sizes
-            rotation_matrix = xy_rotation_matrix(np.deg2rad(self.rotation))
             v_grid, u_grid = np.meshgrid(
                 np.cumsum(np.ones(self._v_count) * self.v_cell_size)
                 - self.v_cell_size / 2,
                 np.cumsum(np.ones(self._u_count) * self.u_cell_size)
                 - self.u_cell_size / 2,
             )
-            xyz = (
-                rotation_matrix
-                @ np.c_[np.ravel(u_grid), np.ravel(v_grid), self.prisms[:, 0]].T
-            ).T
-            xyz += np.asarray(self._origin.tolist())[None, :]
+            uvw = np.c_[np.ravel(u_grid), np.ravel(v_grid), self.prisms[:, 0]]
+            xyz = self.uvw_to_xyz(uvw)
 
             # Instantiate an array of centroids, one for each prism
             indices = (self.layers[:, 0] * self._v_count + self.layers[:, 1]).astype(
@@ -359,6 +355,16 @@ def centroids(self) -> np.ndarray:
 
         return self._centroids
 
+    @property
+    def extent(self) -> np.ndarray:
+        """
+        Compute outer extent of mesh span in world coordinates.
+        """
+        extent = super().extent
+        extent[0, 2] = self.layers[:, 2].min()
+        extent[1, 2] = self.prisms[:, 0].max()
+        return extent
+
     @property
     def n_cells(self) -> int:
         """
@@ -373,6 +379,19 @@ def shape(self) -> np.ndarray:
         """
         return np.array([self._u_count, self._v_count, self.layers.shape[0]])
 
+    @property
+    def span(self) -> np.ndarray:
+        """
+        Upper and lower limits along u, v and w directions.
+        """
+        return np.vstack(
+            [
+                np.sort([0, self.u_cell_size * self.u_count]),
+                np.sort([0, self.v_cell_size * self.v_count]),
+                [self.layers[:, 2].min(), self.prisms[:, 0].max()],
+            ]
+        )
+
     @property
     def u_cell_size(self) -> float:
         """
@@ -395,6 +414,24 @@ def u_count(self) -> np.int32:
         """
         return self._u_count
 
+    def uvw_to_xyz(self, coordinates: np.ndarray) -> np.ndarray:
+        """
+        Apply rotation to the input coordinates.
+
+        Deal with clockwise rotation angle of VPMesh
+
+        :param coordinates: Array of coordinates along the axes, transformed to world coordinates.
+
+        :return: Array of coordinates along the axes, transformed to world coordinates.
+        """
+
+        rotation_matrix = xy_rotation_matrix(-np.deg2rad(getattr(self, "rotation", 0)))
+
+        xyz = (rotation_matrix @ coordinates.T).T
+        xyz += np.asarray(self._origin.tolist())[None, :]
+
+        return xyz
+
     @property
     def v_cell_size(self) -> float:
         """

geoh5py/ui_json/validation.py

@@ -373,6 +373,7 @@ def dependency_type_validation(
 
     dependency = json_dict[name]["dependency"]
     dependency_form = json_dict[dependency]
+
     if "optional" not in dependency_form and not isinstance(data[dependency], bool):
         raise UIJsonError(
             f"Dependency {dependency} must be either optional or of boolean type."
@@ -400,6 +401,7 @@ def mesh_type_validation(name: str, data: dict[str, Any], json_dict: dict[str, A
     """
 
     mesh_types = json_dict[name]["mesh_type"]
+
     obj = data[name]
     if not isinstance(obj, tuple(mesh_types)):
         raise UIJsonError(f"Object's mesh type must be one of {mesh_types}.")
@@ -415,6 +417,7 @@ def parent_validation(name: str, data: dict[str, Any], json_dict: dict[str, Any]
     """
 
     form = json_dict[name]
+
     child = data[name]
     parent = data[form["parent"]]
 

tests/block_model_test.py

@@ -153,6 +153,32 @@ def test_create_block_model_data(tmp_path):
         )
 
 
+def test_block_extents(tmp_path):
+    origin = np.random.randn(3)
+    rotation = np.random.randint(0, 90, 1)
+
+    with Workspace.create(tmp_path / f"{__name__}.geoh5") as workspace:
+        mesh = BlockModel.create(
+            workspace,
+            name="test",
+            u_cell_delimiters=np.array([-0.5, 0.0, 0.5]),
+            v_cell_delimiters=np.array([-0.5, 0.0, 0.5]),
+            z_cell_delimiters=np.array([-0.5, 0.0, 0.5]),
+            origin=origin,
+            rotation=rotation,
+        )
+
+    angle = np.deg2rad(rotation)
+    delta = np.sin(angle) + np.cos(angle)
+    assert np.allclose(
+        mesh.extent,
+        np.c_[
+            origin - np.r_[delta, delta, 1] * 0.5,
+            origin + np.r_[delta, delta, 1] * 0.5,
+        ].T,
+    )
+
+
 def test_shaped_data_values():
     with Workspace() as workspace:
         block_model = BlockModel.create(

tests/grid_2d_test.py

@@ -186,6 +186,36 @@ def test_grid2d_to_geoimage(tmp_path):
             converter.key_to_data(grid, [0, 1])
 
 
+def test_grid2d_extents(tmp_path):
+    origin = np.random.randn(3)
+    rotation = np.random.randint(0, 90, 1)
+
+    with Workspace.create(tmp_path / f"{__name__}.geoh5") as workspace:
+        mesh = Grid2D.create(
+            workspace,
+            origin=origin,
+            u_cell_size=np.r_[1.0],
+            v_cell_size=np.r_[1.0],
+            u_count=1,
+            v_count=1,
+            rotation=rotation,
+            dip=np.random.randint(0, 90, 1)[0],
+        )
+
+    azm = np.deg2rad(rotation)
+    dip = np.deg2rad(mesh.dip)
+    assert np.allclose(
+        mesh.extent,
+        np.c_[
+            origin + np.r_[-np.sin(azm) * np.cos(dip), 0, 0],
+            origin
+            + np.r_[
+                np.cos(azm), (np.sin(azm) + np.cos(azm) * np.cos(dip)), np.sin(dip)
+            ],
+        ].T,
+    )
+
+
 def test_shaped_data_values():
     n_u, n_v = 5, 4
     with Workspace() as workspace:

tests/octree_test.py

@@ -185,3 +185,30 @@ def test_change_octree_cells(tmp_path: Path):
             ValueError, match="New octree_cells array must have the same shape"
         ):
             base_mesh.octree_cells = octree_cells
+
+
+def test_octree_extents(tmp_path):
+    origin = np.random.randn(3)
+    rotation = np.random.randint(0, 90, 1)
+    with Workspace.create(tmp_path / f"{__name__}.geoh5") as workspace:
+        mesh = Octree.create(
+            workspace,
+            origin=origin,
+            u_count=32,
+            v_count=16,
+            w_count=8,
+            u_cell_size=1.0,
+            v_cell_size=2.0,
+            w_cell_size=4.0,
+            rotation=rotation,
+        )
+
+    angle = np.deg2rad(rotation)
+    assert np.allclose(
+        mesh.extent,
+        np.c_[
+            origin + np.r_[-np.sin(angle) * 32, 0, 0],
+            origin
+            + np.r_[np.cos(angle) * 32, (np.sin(angle) + np.cos(angle)) * 32, 32],
+        ].T,
+    )