surface-apps/tests/utils_test.py at f7c3315fb9679a95d8ca69d5666811ae286aee1e · MiraGeoscience/surface-apps · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
#  Copyright (c) 2024-2025 Mira Geoscience Ltd.                                '
#                                                                              '
#  This file is part of surface-apps package.                                  '
#                                                                              '
#  surface-apps is distributed under the terms and conditions of the MIT License
#  (see LICENSE file at the root of this source code package).                 '
# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

import numpy as np
from geoh5py import Workspace
from geoh5py.objects import Curve, Points

from surface_apps.iso_surfaces.utils import interp_to_grid


def get_points(workspace):
    x = np.linspace(0, 10, 11)
    y = np.linspace(0, 10, 11)
    z = np.linspace(-8, 2, 11)
    x_grid, y_grid, z_grid = np.meshgrid(x, y, z)
    vertices = np.c_[x_grid.flatten(), y_grid.flatten(), z_grid.flatten()]
    points = Points.create(workspace, name="my points", vertices=vertices)

    return points


def get_curve(workspace):
    pts = get_points(workspace)
    locs = pts.vertices
    cells = np.c_[np.arange(len(locs) - 1), np.arange(1, len(locs))]
    cell_distances = np.linalg.norm(locs[cells[:, 1]] - locs[cells[:, 0]], axis=1)
    too_far = cell_distances > 2
    cells = cells[~too_far]
    curve = Curve.create(
        workspace,
        name="my curve",
        vertices=locs,
        cells=cells,
    )

    return curve


def test_interp_points_to_grid(tmp_path):
    ws = Workspace(tmp_path / "test.geoh5")
    pts = get_points(ws)
    values = np.zeros(pts.n_vertices)

    right_side = pts.vertices[:, 0] > 5
    values[right_side] = 1

    above_zero = pts.vertices[:, 2] > 0
    values[above_zero] = np.nan

    data = pts.add_data({"data": {"values": values}})
    grid, gridded_data = interp_to_grid(pts, data, resolution=0.5, max_distance=1)
    y_grid, x_grid, z_grid = np.meshgrid(grid[1], grid[0], grid[2])
    pts = Points.create(
        ws,
        name="my grid",
        vertices=np.c_[x_grid.flatten(), y_grid.flatten(), z_grid.flatten()],
    )
    pts.add_data({"my grid data": {"values": gridded_data.flatten()}})

    assert all(gridded_data[x_grid > 6] == 1)
    assert all(gridded_data[x_grid < 4] == 0)
    assert grid[2].max() == 0


def test_interp_curve_vertex_data_to_grid(tmp_path):
    ws = Workspace(tmp_path / "test.geoh5")
    crv = get_curve(ws)
    values = np.zeros(crv.n_vertices)

    right_side = crv.vertices[:, 0] > 5
    values[right_side] = 1

    above_zero = crv.vertices[:, 2] > 0
    values[above_zero] = np.nan

    data = crv.add_data({"vertex data": {"values": values, "association": "VERTEX"}})
    grid, gridded_data = interp_to_grid(crv, data, resolution=0.5, max_distance=1)
    y_grid, x_grid, z_grid = np.meshgrid(grid[1], grid[0], grid[2])
    pts = Points.create(
        ws,
        name="my grid",
        vertices=np.c_[x_grid.flatten(), y_grid.flatten(), z_grid.flatten()],
    )
    pts.add_data({"my grid data": {"values": gridded_data.flatten()}})

    assert all(gridded_data[x_grid > 6] == 1)
    assert all(gridded_data[x_grid < 4] == 0)
    assert grid[2].max() == 0


def test_interp_curve_cell_data_to_grid(tmp_path):
    ws = Workspace(tmp_path / "test.geoh5")
    crv = get_curve(ws)
    values = np.zeros(crv.n_cells)

    right_side = crv.vertices[crv.cells[:, 1], 0] > 5
    values[right_side] = 1

    above_zero = crv.vertices[crv.cells[:, 1], 2] > 0
    values[above_zero] = np.nan

    data = crv.add_data({"cell data": {"values": values, "association": "CELL"}})
    grid, gridded_data = interp_to_grid(crv, data, resolution=0.5, max_distance=1)
    y_grid, x_grid, z_grid = np.meshgrid(grid[1], grid[0], grid[2])
    pts = Points.create(
        ws,
        name="my other grid",
        vertices=np.c_[x_grid.flatten(), y_grid.flatten(), z_grid.flatten()],
    )
    pts.add_data({"my other grid data": {"values": gridded_data.flatten()}})

    assert all(gridded_data[x_grid > 6] == 1)
    assert all(gridded_data[x_grid < 4] == 0)
    assert grid[2].max() == -0.5