Slicing the Mesh Arbitrarily

[aluthfian]

Hi all, I just want to propose a feature addition (if there isn't any) for future releases of discretize: slicing the mesh using a straight line of arbitrary direction. I propose this feature as an addition to the already available discretize.TensorMesh.plot_slice. The draft code for this functionality is given below.

from discretize import TensorMesh, TreeMesh
from discretize.utils import active_from_xyz
import numpy as np
import pyvista as pv
from SimPEG.utils import mkvc, model_builder, plot2Ddata, mat_utils

def interpolate_points(start, end, N):
    """
    Interpolates N points along a straight line between start and end points.
    
    Parameters:
        start (tuple or list): The starting point (x, y).
        end (tuple or list): The ending point (x, y).
        N (int): The number of points to interpolate.
        
    Returns:
        numpy.ndarray: An array of interpolated points.
    """
    # Convert the start and end points to numpy arrays
    start = np.array(start)
    end = np.array(end)
    
    # Generate linearly spaced values between 0 and 1
    t_values = np.linspace(0, 1, N)
    
    # Interpolate points
    interpolated_points = np.outer(t_values, end - start) + start
    
    return interpolated_points

mesh_data = TreeMesh.read_UBC('msh_file_path')
model_data = mesh_data.read_model_UBC('mod_file_path')
mesh_data_vtk = mesh_data.to_vtk({'model':model_data})
mesh_data_vtk.set_active_scalars('model')

xy_coord = interpolate_points(line_start, line_end, N)
z_coord = np.linspace(lower_z, upper_z, num=len(xy_coord)) #lower_z<upper_z
xyz_coord = np.array([np.c_[xy_coord, depth*np.ones(xy_coord.shape[0])] for depth in np.flip(z_coord)])
xyz_coord = xyz_coord.reshape(xy_coord.shape[0]*z_coord.shape[0],3)

# Get voxel value
voxel_id = mesh_data_vtk.find_closest_cell(xyz_coord)
x_coord = np.array([mesh_data_vtk.get_cell(index).center[0] for index in np.unique(voxel_id)])
z_coord = np.array([mesh_data_vtk.get_cell(index).center[2] for index in np.unique(voxel_id)])
model_value_at_cellCenters = np.array([mesh_data_vtk.active_scalars[index] for index in np.unique(voxel_id)])
model_value_at_regXYZ = mesh_data_vtk.active_scalars[voxel_id]
model_value_at_regXYZ = model_value_at_regXYZ.reshape(xy_coord.shape[0],zo.shape[0])

# Plotting
fig, ax = plt.subplots()
ax.imshow(model_value_at_regXYZ,
         extent=(xyz_coord[:,0].min(), xyz_coord[:,0].max(), xyz_coord[:,2].min(), xyz_coord[:,2].max()),
         origin='upper', aspect='auto')

The final result will look like this:

Thanks for reading!

[micmitch]

Hi @aluthfian! I developed some code for doing this a couple of years back, which also uses the vtk model object. The last time I tried running this plotting code though I remember having problems loading the vtk model into pyvista. I haven't dug into this to try and sort it out. I'll try to revive my code and clean it up enough to share.

It would be very interesting to experiment with both of our plotting codes and combine the best of both. I seem to remember my approach being a little more complicated and requiring many more lines of code so your approach maybe cleaner. I know that in mine I opted to plot depth versus the along profile distance instead of x or y coordinates.

[micmitch]

I have definitely found the ability to plot arbitrarily oriented sections to be valuable. Here is an example plot from a recent paper.

[aluthfian]

@micmitch Good idea! In my earlier version of this, I used vtk's slice along line capability and then retrieve the points on the slice using ctp(). This was very slow (needs tens of seconds to run), memory intensive, and results in a longer code. The approach I proposed here is faster to run (a few seconds only) and keep the voxel-like shape intact for honesty and transparency (although interpolation option may be possible to implement so the plot can look "smooth").

Changing the x axis into distance instead of Northing/Easting should also be given as an option as well, like:

1. X or Y distance: Easting/Northing minus their minimum value.
2. Actual distance uses np.linalg.norm

I am interested on how you produced your smooth plot since I can only to use scipy.interpolation.griddata(..."linear"...) for my case (other interpolation methods were just failed).

(Sorry for the closing and reopening of the same issue, commenting on GitHub mobile and small screen is not that convenient!)