Voxelize model with known block size

[javoha]

Hello everyone,

lately I have been trying to create a voxelized ("block") model. The data I want to visualize is based on x,y,z data that was grided before, but does not cover a whole cuboid space. When the orginial grid that the data is retrieved from is regularly spaced, I am able to get a nice voxelized model using PVGeo.filters.VoxelizePoints() . However, in my real application case the orginal grid is not cubic, but instead I have different length in different directions. I do know the resolution of the original grid (thus also the voxel size I am aiming for). When using PVGeo.filters.VoxelizePoints() I get a result that doesnt cover the whole space - there are gaps between voxels and maybe also some overlapping.

Maybe I am missing a simple solution (maybe even jsut in pyvista) for this as I know the original resolution. Some help or ideas would be highly appreciated.

Thanks a lot for any help and support and of course for your time!

[banesullivan]

It's tough to tell what's going on without the data... Can you share a file and sceenshots?

[banesullivan]

FYI, the pyvista support forum gets more visibility than this repo: https://github.com/pyvista/pyvista-support

[javoha]

Hi @banesullivan,

thanks for taking a look. Should I move the issue to Pyvista or just keep it here for now?

Maybe for some context - I am on working in gempy. My aim is to populate single units of a gempy model with property values (e.g. density) and display the result as a voxel model. Here a simple test model of a foreset kind of setting that I am using (2.5D):

I extract the grid data for the red unit (domain of interest) and populate it with density values using SGS or Kriging. If my orginal gempy model had a evenly spaced resolution (e.g extent of model [1000,1000,1000] and resolution [100,100,100] in this example, I get a nice result using the PVGeo.filters.VoxelizePoints():

If my grid grid is not evenly spaced though (here extend still [1000,1000,1000] but res [100,50,100]), the resulting voxel volume does not fill the whole space, as seen here:

Probably I am missing a simple pyvista solution to this, just wasnt able to find a nice one yet.

Again thanks a lot for checking it out, I am grateful for any help or hints!

[banesullivan]

Keep the issue here for now. There is likely a simpler solution with PyVista. Can you show me/send me whatever input you're passing to PVGeo.filters.VoxelizePoints?

e.g.

points.save('points.vtk')

[banesullivan]

PVGeo's VoxelizePoints does not currently handle non-uniform grid cells

[javoha]

The unstrucutrued grid for pyvista is created based on x,y, z and property value (called field) using the workflow from GSTools. Here a screenshot of the resulting unstructured grid for the second example shown before:

This is also what I feed the voxelizer:

grid = PVGeo.filters.VoxelizePoints().apply(pc)

I also attached the points file as you suggested:

points

Thanks again for helping.

[banesullivan]

If my grid grid is not evenly spaced though (here extend still [1000,1000,1000] but res [100,50,100]), the resulting voxel volume does not fill the whole space, as seen here:

Ah, actually I just realized you know the grid spacing ahead of time. VoxelizePoints cannot figure out the non-uniform spacing on its own but it can handle non-uniform spacing if it is informed via the dx, dy, dz, and estimate arguments.

Try:

grid = PVGeo.filters.VoxelizePoints(dx=10, dy=20, dz=10, estimate=False).apply(pc)

FYI, this is not documented. I only know because I wrote the algorithm.

Also, I found the resolution of the points by:

>>> spacing = lambda arr: np.unique(np.diff(np.unique(arr)))
>>> spacing(pc.points[:,0]), spacing(pc.points[:,1]), spacing(pc.points[:,2])
(pyvista_ndarray([10.]), pyvista_ndarray([20.]), pyvista_ndarray([10.]))

---

[javoha]

Perfect, works like a charm. Thanks for your help. I suppose this means there is also not a much easier solution just using pyvista?

[banesullivan]

I thought there might be by making a StucturedGrid of these data but they are indeed unstructured, so this is the best I can think of for now

[RichardScottOZ]

Thanks!

I have a similar sort of problem to do with density - where I want to make a density model for which I know the resolution I want - e.g. 30000 25000 75 - and have a bunch of surfaces of rock units to use (and make) - voxelise? Not sure of the best/most efficient way to do it.

Basically over the same area - e.g. state of SA - but the surfaces data will be in whatever coordinate system - and coarser resolution vertically, Blocks will be cubes though! e.g. 50x50x50.

I merged the surfaces that I was working on just to see what it would look like - whether it is better to do them all separately?

[RichardScottOZ]

Eventually into a floating point Zarr datastore as such

In 2D

[RichardScottOZ]

There'll be more surfaces to go into the above, just haven't made it yet. 76 Z being like 1300m above ground and 2500 under in this case.

[RichardScottOZ]

Making a gempy version at whatever resolution to look at will be good, too @javoha

[RichardScottOZ]

One of the various 100+ surfaces looks like this

[RichardScottOZ]

I 'high density number' voxelised that to have a look

[banesullivan]

@RichardScottOZ, I'm not following what you're hoping to do... I'd suggest opening a new support topic in the PyVista forum

[banesullivan]

Closing this as the original issue using PVGeo.filters.VoxelizePoints is resolved

[RichardScottOZ]

Ok, will do, thanks!

[RichardScottOZ]

Done in https://github.com/pyvista/pyvista-support/issues/404