pythonpyvista

PyVista ray_trace for structured grids?

What I want to do

For a boundary meteorological application in complex terrain I would like to check whether a sun ray hits the topography so that the cell under consideration is shaded (see photo)

What I have done and what the problem is

PyVista has in form of a lightweight ray tracing procedure the perfect template for my need. It can be found here and here. These examples are based on sphere = pv.Sphere().

But I use a structured grid myMesh = pv.StructuredGrid() and this rises the error message : AttributeError: 'StructuredGrid' object has no attribute 'ray_trace'

My question

Is there a work-around so that I can use PyVista.ray_trace with a structured grid?

enter image description here

Solution

  • Thanks to the hint from @Andras Deak -- Слава Україні I was able to find a solution.

    Since I have spent some time on it, I give the whole code part. Maybe someone else will find it useful too.

    import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv

#---- 1. define a ray
start = [0.0,  0.0,  -1]
stop  = [0.0,  0.0,   1]
ray = pv.Line(start, stop)                                     


#---- 2. generate a 2-dimensional numpy grid (x,y,z)
nx, ny = 4, 4
ix = np.linspace(-1, 1, nx)
iy = np.linspace(-1, 1, ny)
xx, yy = np.meshgrid(ix, iy, indexing='ij')
zz = np.zeros_like(xx)
xyz = np.hstack((xx.reshape(-1,1), 
                 yy.reshape(-1,1), 
                 zz.reshape(-1,1)))


#---- 3. build the different pyvista data objects
cMesh = pv.PolyData(xyz)                                      # cloud of data
sMesh = pv.StructuredGrid(xx,yy,zz)                           # structured data
eMesh = sMesh.extract_surface()            # <<---- important MAGIC surface extracting


#---- 4. run the ray tracing/intersection
iPoints, ind = eMesh.ray_trace(start, stop, first_point=True)  # Runing the ray tracing gives the intersection points
ipMesh = pv.PolyData(iPoints)                                  # wrap the intersection points in a mesh


#---- 5. plot
p = pv.Plotter()
p.set_background('navy', top='teal')

p.add_mesh(ray,   color="orange",    line_width=5)           # add the ray
p.add_mesh(cMesh, color="red",       point_size=15)          # add the data as a cloud
p.add_mesh(sMesh, color="lightgrey", show_edges=True)        # add the structured grid
p.add_mesh(ipMesh,color="yellow",    point_size=30)          # add the computed intersetion points

p.camera.zoom(1.5)
p.show()

print('Intersection points:'); print(iPoints)

    enter image description here

    Update

    The ray_trace() computation seems to be very fast (the graphical display, which is hardware dependent, uses more ressources).

    enter image description here