you need to change the 'array name' in this line points = point_data.GetAbstractArray('array name') to one of the names that is printed out.
Sorry on my version I changed this, but I pasted the original above.
In my version it reads:
vel = vtk_to_numpy(point_data.GetAbstractArray('Velocity'))
points = vtk_to_numpy(data.GetPoints().GetData())how big are your files? can you post one?
Hi Justin,
I actually had a similar question. Could you help me with it? Thanks!
The website does not allow me to post vtu file, but I got my vtu using this mfx file. vtu.mfx (13.6 KB)
Then I used the following code to extract data
import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(“BACKGROUND_IC_0005.vtu”)
reader.Update()
data = reader.GetOutput()
point_data = data.GetPointData()
for i in range(point_data.GetNumberOfArrays()):
print(i, point_data.GetArrayName(i), ‘Components:’, point_data.GetArray(i).GetNumberOfComponents())
But I didn’t see anything printing out. Did I misunderstand somewhere on how to use the code?
Thanks,
Jack
There’s no point data in that file (it contains cell data only). It’s useful to open an interactive Python session for data exploration:
% python
>>> import vtk
>>> reader = vtk.vtkXMLUnstructuredGridReader()
>>> reader.SetFileName("BACKGROUND_IC_0005.vtu")
>>> reader.Update()
>>> data = reader.GetOutput()
>>> data.GetPointData().GetNumberOfArrays()
0
>>> data.GetCellData().GetNumberOfArrays()
1
Also note that although you cannot post VTU files, you can post zipfiles, so just use ‘zip’ to create a zipfile if you need to upload data files.
– Charles
Thanks a lot Charles! I see where the problem is now. I got the cell array data using triangles=vtk_to_numpy(data.GetCells().GetData()). However, I am very new to vtk and still got really confused by the result. Could you give me some hint on what the result mean? I noticed that the my result of vtk_to_numpy(data.GetCells().GetData()) was a series of numbers with many 4s. Each 4 was followed by 4 integers. What information is it providing? Thank you!
I checked out the kitware tutorials that you sent to me, but I didn’t find how GetCells() or GetPoints() work. What are the points and cells? In my understanding, point data stores information about where the points are. Namely it stores arrays of tuples that represent coordinates of points. Then cell data can be parsed into tuples of integers, and each integer represents the index of a point in the point data. Did I misunderstand how vtk works?
Thanks a lot,
Jack
Jack -
You are calling vtk_to_numpy on the wrong object. You need to call GetCellData, not GetCells.
Try this:
import vtk
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName("BACKGROUND_IC_0005.vtu")
reader.Update()
output = reader.GetOutput()
data = output.GetCellData()
print(data.GetNumberOfArrays())
arr = data.GetArray(0)
print(arr)
print(arr.GetTuple(0))
print(vtk_to_numpy(arr))
these values are just the x-velocity components that you’re writing in the VTK file. If you had selected more variables the tuples would contain multiple values. CORRECTION: each variable is stored in its own array.
– Charles
1 Like
This is how this histogram viewer works in the gui, if it is a vtu file, use read_griddata, if it is a vtp file, use read_partdata. The functions will return a dict of variables names with some other info including the actual array, data. Then just use the vtk_to_numpy function on that data array.
import vtk
def read_griddata(path):
ugrid_reader = (
vtk.vtkXMLPUnstructuredGridReader()
if path.suffix == "pvtu"
else vtk.vtkXMLUnstructuredGridReader()
)
ugrid_reader.SetFileName(path)
ugrid_reader.Update()
data = ugrid_reader.GetOutput()
cell_data = data.GetCellData()
array_info = {}
for i in range(cell_data.GetNumberOfArrays()):
array = cell_data.GetArray(i)
n_comp = array.GetNumberOfComponents()
name = cell_data.GetArrayName(i)
array_info[name] = {
'magnitude': array.GetRange(-1),
'components': n_comp,
'range': [array.GetRange(i) for i in range(n_comp)],
'data': cell_data.GetAbstractArray(name),
}
return array_info
def read_partdata(path):
particle_reader = (
vtk.vtkXMLPPolyDataReader()
if path.suffix == ".pvtp"
else vtk.vtkXMLPolyDataReader()
)
particle_reader.SetFileName(path)
particle_reader.Update()
data = particle_reader.GetOutput()
point_data = data.GetPointData()
array_info = {}
n_tuples = None
for i in range(point_data.GetNumberOfArrays()):
array = point_data.GetArray(i)
n_comp = array.GetNumberOfComponents()
n_tuples = array.GetNumberOfTuples()
name = point_data.GetArrayName(i)
array_info[name] = {
'number_of_tuples': n_tuples,
'components': n_comp,
'range': [array.GetRange(i) for i in range(n_comp)],
'magnitude': array.GetRange(-1),
'data': point_data.GetAbstractArray(name),
}
return array_info
if __name__ == "__main__":
vtp_arrays = read_partdata('/path/to/output.vtp')
1 Like
Yeah that works. I can retrieve the data now. Thanks a lot Charles! One more follow up question though. If I selected multiple variables, how can I figure out which value corresponds to which variable in each tuple?
Thanks,
Jack
I was mistaken when I said VTK would store a tuple. Each variable goes into its own array, and every array has a name.
Use GetArrayName:
>>> print(data.GetNumberOfArrays())
3
>>> print(data.GetArrayName(0))
U_G
Note that this is contained in the code Justin posted above.
Also if you’re trying to learn VTK, at an interactive Python prompt you can type
>>> help(data)
(where data is a VTK data object) you will get a list of all the methods available with a description of each. Or consult a VTK manual.
I see. That works now. Thanks Charles! I was trying out Justin’s code and it gave me an error saying NameError: name ‘GridArray’ is not defined. I guess I was missing some library. But your solution works now.
OK cool I will use the help() function.
Thank you,
Jack
Yes, GridArray and ParticleArray are not part of VTK, they are defined elsewhere. The code Justin posted was incomplete, just intended as a code sample. Sorry for the confusion!
Sorry. You can replace those “data classes” with regular dictionaries. I edited the code above.
2 Likes
Hi, Justin, May I ask a question? I have used the code obtaining the information about “Diameter” , but how can I obtain the coordinate of each particle from vtp files? The coordinate is not a component name.
you could store the particle x y z coordinates as 3 distinct DES_USR_VAR that can thus be saved as vtk along with velocity, diameter etc
Thanks very much. But in this case, I will recalculate all the cases.
The points (particle locations) are stored in the vtp files differently then the variable arrays. You can access the points like this:
import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy
particle_reader= vtk.vtkXMLPolyDataReader()
particle_reader.SetFileName(fname)
particle_reader.Update()
data = particle_reader.GetOutput()
# get the positions!
particle_positions = vtk_to_numpy(data.GetPoints().GetData())
FYI, I learned most of this from reading the VTK docs and the examples.
Thanks very much for your valuable instructions. It really helps me.
@onlyjus How can I get the velocity of the particles? I can get the positions of the particles by the piece of code you provided, but I am not able to get the velocity?
If you are saving the particle velocity in the VTK files, you can access that variable for each particle.
To see the array names, you can do something like this (untested):
particle_reader = (
vtk.vtkXMLPPolyDataReader()
if path.suffix == ".pvtp"
else vtk.vtkXMLPolyDataReader()
)
particle_reader.SetFileName(path)
particle_reader.Update()
data = particle_reader.GetOutput()
point_data = data.GetPointData()
for i in range(point_data.GetNumberOfArrays()):
print(point_data.GetArrayName(i))
Then to get the actual data array, do something like this (with the array name you want, printed with the above code):
array = point_data.GetAbstractArray(name)