I have an array which I wish to add as another data field to an existing vtk file - but im struggling to get it to work. I’m using python, specificaly the build of python and vtk that comes with paraview 5.6.This is what i have so far:
Then I have a for loop to add the data values to the array:
for value in range(np.shape(MyArray)[0]):
MyVTKObject.GetPointData().GetArray(“NewDataField”).SetTuple(value,MyArray[value])
I can then access the field using something like
MyVTKObject.GetPointData().GetArray(NewDataField).GetTuple(n)
However, i think i am doing something wrong, because it does not show as an additional data field in Paraview, and is not save when i save the vtk object to file. Any ideas?
How are you saving that object? STL’s for example cannot store data arrays. Best use .vtp files, which you can write using the vtk.vtkXMLPolyDataWriter.
The way you create the float array looks correct. Also you mentioned that your check with GetTuple() was successful. So the problem must be caused elsewhere.
For comparison, here some code I found in my library to add array data.
# data: n rows, d cols
# source: poly data, with n points
data = np.asarray(data)
data = np.atleast_2d(data)
nPoints = source.GetNumberOfPoints()
nRows, nCols = data.shape
assert(nPoints==nRows)
array = vtk.vtkDoubleArray()
array.SetName("arrayName")
if nDims > 1:
# For a multi-dimensional array.
array.SetNumberOfComponents(nRows)
array.SetNumberOfTuples(nCols)
for x in zip(range(nElms), data.T):
array.SetTuple(*x)
else:
# For a one-dimensional array.
array.SetNumberOfValues(nElms)
for x in zip(range(nElms), data.T):
array.SetValue(*x)
dataset = source.GetPointData()
dataset.AddArray(array)
# In case you want to make the array the current one.
isScalar = array.GetNumberOfComponents() == 1
if isScalar:
dataset.SetActiveScalars(array.GetName())
else:
dataset.SetActiveVectors(array.GetName())
Make sure you call Modified() on the array/pointdata/polydata object.
Also note that setting elements one by one would be extremely inefficient. In general, you should never iterate through all points of a mesh or all voxel of a volume in Python, but you should always operate on entire arrays or matrices. You can get/set/modify the entire array at once by mapping it to a numpy array.
See complete implementation for both mapping a numpy array to a VTK array and calling the necessary Modify() methods here:
Thumbs up for the use of numpy containers. A pity that the vtk examples use that slow paradigm, so that the python folks are misled into writing unnecessarily slow code… I didn’t know that you can modify the returned numpy-array directly, that’s a nice!
Just want to point out that this question is about vtkPointData (or vtkCellData), as returned by poly.GetPointData() or poly.GetCellData(). Your pasted code acts on the points themselves, as returned by poly.GetPoints(). No big deal, just a couple of modifications required.
@Shahid My sample code simplifies to the following:
data = np.asarray(data)
data = np.atleast_2d(data)
nPoints = source.GetNumberOfPoints()
assert(nPoints==len(data))
from vtk.util.numpy_support import numpy_to_vtk
# Setting deep=True is a bit safer, but you lose a bit of performance.
array = numpy_to_vtk(data, deep=True)
array.SetName("arrayName")
dataset = source.GetPointData()
dataset.AddArray(array)
dataset.Modified() # As suggested by Andras
I’m not sure if allocating an array with numpy and pass ownership to VTK is safe, as you can have so many different kind of numpy arrays, with different memory layouts and different memory management policies. If you deep-copy the array content then of course it is safe but then the two arrays are completely independent duplicates and they are not kept in sync. Therefore, I always allocate memory in VTK if I want VTK to manage the object. Then the only thing to pay attention to is not to perform any operations that would change the memory layout of the array.
Thanks - the numpy utilities are very useful. Didnt know they exisited.
In the end my problem was more to do with the client server architecture of Paraview - I have solved it by writing the vtk using vtk.vtkXMLPolyDataWriter (before i was using one of paraviews functions) and reloading the vtk.