HOWTO: convert volumes (Face3D from dxf) to proportions in a regular grid (e.g., proportions of clays)

[antoinecollet5]

Hi !

What I am trying to achieve:

I have a dxf composed of around 70k Face3D entities. These represent volumes of clays. I can load it with tools like ezdxf without problem.

Now I have a regular grid (sugarbox kind with Parallelepiped grid cells) which is referenced in the same space (basically a numpy array with centroid coords x-y-z). What I want to achieve is to iterate over the Face3D and compute the intersection with the grid. Basically, I would like to compute for each grid cell the proportion of coverage by these volumes, i.e., the proportion of clays in each grid cells. Or at least, does that cell contains clay (a boolean).

I know how to do that in 2d using libs like shapely, but what would be the right tool for such a task in 3D ? Can gempy helps since I have seen it deals with both vectors and raster-like arrays.

Note: I have tried with proprietary tools like IsatisNEO, but I guess their algs are not that optimized and even after hours, the job is not complete. I guess I can do a bit better by working with subgrids (using mask based on the bounding box of each Face3D), even for tools in pure python.

Thank you in advance for the help.

[javoha]

Hi @antoinecollet5,
thanks for your question. Let me just reformulate your problem to see if I got it right: You have a set of clay volumes, that are defined by Face3D entities (which I guess are surface meshes around these volumes?). You would either (base version) like to know for each point of your regular grid if it lies within a clay volume or (advanced version) know for each grid cell (volume) how much of it (percentage) lies within a caly volume.
Is that correct?

If it is: How many clay volumes do you have (70k faces does not mean 70k clay bodies, right?)

While GemPy should in theory be able to solve the base version of this problem by modeling the clay volumes as elements and then computing a solution for your grid cells, in reality it will probably not work, as gempy works best for layered elements. While it is possible to brute force other shapes, this gets out of hand (and impossible) very quickly for a huge number of unconnected volumes (e.g. clay lenses).
Your problem to me seems more like a "computer graphics" problem of calculating the intesections between meshes and grids. We use pyvista for the 3D visualization and they might have what you need (e.g. something like this, if though this is not the exact thing).

Let me know if I got your problem correctly, if not maybe I am able to help if you tell me where I went wrong.

Cheers,
Jan

[antoinecollet5]

Hi @javoha !

Thanks for the insightful answer.
You got the problem right 100%. I have to admit that I don't know exactly how many volumes I have. I will investigate that. Thank you for the tips with pyvista. I will look into it starting from paraview first.

Another option I got from the ezdxf community is to go for cadquery.

I will post a solution here if I found something performant as it might interest the community of gempy.

Cheers,
Antoine

[javoha]

Sounds great, I would definitely like to see the solution if/when you find it.