Custom grid solution contains floats that do not cleanly convert to integer lithology identifiers

[benk-mira]

I have used the set_custom_grid method to compute solutions on a pre-existing grid. I am retrieving the solution from Solutions.raw_arrays.custom but I find that the solution is all floats and contains values that are not clearly in one of the categories:

Many of the floats are close to one of the integer values and converting to integer takes care of most of issues, but I am left with some bad values in the array:

I expected that the solution in the custom array would be integer valued.

Here is my code to reproduce the slice of data pictured above.

import numpy as np
import gempy as gp
from geoh5py import Workspace
from geoh5py.objects import BlockModel
import matplotlib.pyplot as plt
from scipy.spatial import cKDTree

color_generator = gp.core.data.ColorsGenerator()
extents = [311730, 318467, 6068686, 6076023, -750, 350]
x = np.arange(extents[0], extents[1], 50)
y = np.arange(extents[2], extents[3], 50)
z = np.arange(extents[4], extents[5], 50)
X, Y, Z = np.meshgrid(x, y, z)
grid = np.c_[X.flatten(), Y.flatten(), Z.flatten()]

# Create structural frame

elements = [
    gp.core.data.StructuralElement(
        name="unit3",
        color=next(color_generator),
        surface_points=gp.data.SurfacePointsTable.from_arrays(
            x=[314283.60, 314436.64, 314793.11, 314593.63, 314503.53, 314624.43],
            y=[6074135.65, 6072641.68, 6071800.08, 6071350.68, 6070537.54, 6069821.63,],
            z=[301.47, 301.57, 301.60, 301.64, 301.67, 301.72],
            names=["unit3"]*6,
        ),
        orientations=gp.data.OrientationsTable.initialize_empty(),
    ),
    gp.core.data.StructuralElement(
        name="unit2",
        color=next(color_generator),
        surface_points=gp.data.SurfacePointsTable.from_arrays(
            x=[313803.07, 314057.46, 314030.63],
            y=[6073334.87, 6071265.75, 6070534.66],
            z=[301.52, 301.63, 301.67],
            names=["unit2"]*3,
        ),
        orientations=gp.data.OrientationsTable.from_arrays(
            x=[313632.11],
            y=[ 6072104.38],
            z=[301.59],
            G_x=[0.96671406],
            G_y=[0.18791018],
            G_z=[0.17364818],
            names=["unit2"]
        ),
    ),
    gp.core.data.StructuralElement(
        name="unit1",
        color=next(color_generator),
        surface_points=gp.data.SurfacePointsTable.from_arrays(
            x=[316326.99, 315571.08],
            y=[6070539.09, 6072716.70],
            z=[301.67, 301.55],
            names=["unit1"] * 2,
        ),
        orientations=gp.data.OrientationsTable.from_arrays(
            x=[315905.12],
            y=[6071471.84],
            z=[301.62],
            G_x=[0.94177634],
            G_y=[0.28792992],
            G_z=[0.17364818],
            names=["unit1"]
        ),
    ),
]
structural_group = gp.core.data.StructuralGroup(
    name="Series1",
    elements=elements,
        structural_relation = gp.core.data.StackRelationType.ERODE,
    )
structural_frame = gp.core.data.StructuralFrame(
    structural_groups=[structural_group], color_gen=color_generator
)

# Compute model
geo_model = gp.create_geomodel(
    project_name="test",
    extent=extents,
    refinement=4,
    structural_frame=structural_frame,
)


gp.set_custom_grid(geo_model.grid, grid)
solution = gp.compute_model(geo_model)
model = solution.raw_arrays.custom

def plot_slice(ax, elevation):
    ind = grid[:, 2] == elevation
    ax.imshow(model.astype(np.int32)[ind].reshape(len(y), len(x)))

fig, ax = plt.subplots(figsize=(10,10))
plot_slice(ax, -50)

I'm running Windows11.

[benk-mira]

..So I have discovered that if I replace the integer conversion (.astype(np.int32)) with proper rounding (.round(0)) the image is much neater. I'm not sure this strategy is perfect though.. Is it expected that the raw_arrays.custom array is float valued? And if so, what's the recommended solution for this issue?

[javoha]

Hi @benk-mira,
thanks for your question. First of all I was able to reproduce your model. I am not 100% sure how gempy handles custom grids at the moment. I think you are correct, that the geo_model.Solutions.raw_arrays.custom should be in an integer format (or at least have a lith_block attribute that is in that format). I think this is just an oversight in gempy version 3. Your rounding solution seems to be the right approach, but I will need to double check that.

Apart from that I have two questions regarding your model: You create the structural frame but I think your elements are not ordered correctly. It still seems to work in this example but best order them correctly to not mess up the algorithm.
Second: Your custom grid is a regular grid that you can easily create within gempy. This would allow you too also use the gempy plotting with the right coordinates etc. . I just did:

geo_model = gp.create_geomodel( project_name="test", extent=extents, resolution=[len(x), len(y), len(z)], structural_frame=structural_frame )

With gempy-viewer you then get:

Let me know if this helps or if you have further questions.
Cheers,
Jan

[benk-mira]

Thanks for getting back to me - Yes you're right - I have adapted my example here from something I cooked up elsewhere and messed up the order. I used a regular grid to simplify things for the example, but I am interested in the custom grid primarily to be able to pass octree mesh centers in and get models I can use for geophysical simulation with SimPEG.

I do have one more question. While looking into this issue, I couldn't seem to find a straightforward way to access Gempys internal octree cells centers and the computed solution on those centers. Any suggestions for this?

I'll use the round solution for now, but look for an update in the future. My initial results with Gempy have been very positive overall, and I'm looking forward to exploring more!

Cheers,
Ben

[javoha]

The octree result is a little hidden: You can find the coordinates in geo_model.solutions.octrees_output[i].grid_centers.values (i indicating the chosen level) and the computed result in geo_model.solutions.octrees_output[i].outputs_centers[0].final_block.

The following code (using pyvista, gives you a plot of the centers with their corresponding lith block value:

# Compute model
geo_model = gp.create_geomodel(
    project_name="test",
    extent=extents,
    refinement=7,
    structural_frame=structural_frame
)

gp.compute_model(geo_model)


gpv.plot_2d(geo_model, direction="z")
gpv.plot_3d(geo_model)

# Create object containing back transformed octree grids per level
corner_list = []
center_list = []

for i in range(len(geo_model.solutions.octrees_output) - 1):
    corner_list.append(
        geo_model.input_transform.apply_inverse(geo_model.solutions.octrees_output[i].grid_corners.values))
    center_list.append(
        geo_model.input_transform.apply_inverse(geo_model.solutions.octrees_output[i].grid_centers.values))

octree_levels = np.array(corner_list, dtype=object)
octree_centers = np.array(center_list, dtype=object)

# Plotting octree level with corresponding lith_block (refinement-1)
level_show = 6

plotter = pv.Plotter()

# Plot each level of octree refinement
for level in range(level_show):
    # centers
    plotter.add_mesh(pv.PolyData(octree_centers[level]),
                     render_points_as_spheres=True,
                     point_size=5,
                     scalars=geo_model.solutions.octrees_output[level].outputs_centers[0].final_block,
                     cmap="viridis")

    # corners
    # plotter.add_mesh(pv.PolyData(octree_levels[level]),
    #                  render_points_as_spheres=True,
    #                  point_size=5,
    #                  scalars=geo_model.solutions.octrees_output[level].outputs_corners[0].final_block,
    #                  cmap="viridis",)



# Set the bounds and grid of the plotter
plotter.show_bounds(bounds=geo_model.grid.extent,
                    location="furthest",
                    grid=True)

# Display the interactive plot
plotter.show()

Gives somethin like this for your model: