21

I've been given a legacy format vtk file (I think its an unstructured grid) and I'd like to read it in with python and output a .npy file instead, since I know how to deal with that.

The file is a dump from ATHENA and so has density, velocity, magnetic field along with the coordinates.

I'm very much a procedural programmer, so all these objects are confusing...

Ben Jackel
  • 429
  • 1
  • 3
  • 6

6 Answers6

21

Here is the solution that I came up with, the trick was turning on ReadAllVectorsOn().

import numpy
from vtk import vtkStructuredPointsReader
from vtk.util import numpy_support as VN

reader = vtkStructuredPointsReader()
reader.SetFileName(filename)
reader.ReadAllVectorsOn()
reader.ReadAllScalarsOn()
reader.Update()

data = reader.GetOutput()

dim = data.GetDimensions()
vec = list(dim)
vec = [i-1 for i in dim]
vec.append(3)

u = VN.vtk_to_numpy(data.GetCellData().GetArray('velocity'))
b = VN.vtk_to_numpy(data.GetCellData().GetArray('cell_centered_B'))

u = u.reshape(vec,order='F')
b = b.reshape(vec,order='F')

x = zeros(data.GetNumberOfPoints())
y = zeros(data.GetNumberOfPoints())
z = zeros(data.GetNumberOfPoints())

for i in range(data.GetNumberOfPoints()):
        x[i],y[i],z[i] = data.GetPoint(i)

x = x.reshape(dim,order='F')
y = y.reshape(dim,order='F')
z = z.reshape(dim,order='F')
Guillaume Jacquenot
  • 11,217
  • 6
  • 43
  • 49
Ben Jackel
  • 429
  • 1
  • 3
  • 6
12

meshio (a project of mine) knows the VTK format, so you could simply

pip install meshio

and then

import meshio
mesh = meshio.read('file.vtk')
# mesh.points, mesh.cells, ...
Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
  • can we save to .csv or .txt using meshio? – npm Sep 11 '19 at 03:11
  • csv and txt are no mesh formats, so no, not out of the box. You can always read the data in and just dump it though. – Nico Schlömer Sep 11 '19 at 06:56
  • @npm The coordinates and field data can be readily determined from the mesh: coords, data = mesh.points, mesh.point_data['']. It is then easy to write these to file in whatever format you wish – Dai May 17 '21 at 19:44
10

Here's a script that reads polygon data into numpy arrays from a VTK file using the VTK Python SDK:

import sys

import numpy
import vtk

reader = vtk.vtkPolyDataReader()
reader.SetFileName(sys.argv[1])
reader.Update()

polydata = reader.GetOutput()

for i in range(polydata.GetNumberOfCells()):
   pts = polydata.GetCell(i).GetPoints()    
   np_pts = numpy.array([pts.GetPoint(i) for i in range(pts.GetNumberOfPoints())])
   print np_pts
jterrace
  • 64,866
  • 22
  • 157
  • 202
  • I have been able to get to the GetOutput() stage using vtkDataSetReader() instead of vtkPolyDataReader, but I am still confused about how to get the information out. GetNumberOfPoints and GetNumberOfCells seems to give me sensible numbers for how big I think the arrays should be, but I still don't know how to pull all the variables out. Is there a way to get information on what exactly the vtk has in it, and in what form? In an understandable way? – Ben Jackel Aug 10 '12 at 19:34
  • If you pick one of the files here http://people.sc.fsu.edu/~jburkardt/data/vtk/vtk.html I can try and help extract the format. There are so many different formats. – jterrace Aug 10 '12 at 19:40
  • reader.IsFileStructuredPoints() returns 1, the other options return 0, so I'm going out on a limb and saying its a legacy vtk of Structured Points. – Ben Jackel Aug 10 '12 at 20:45
  • You probably want to use vtkStructuredPointsReader then? – jterrace Aug 10 '12 at 20:57
5

Here's a short snippet for reading points from legacy VTK files:

import numpy as np
import vtk

filename = 'model.vtk'

reader = vtk.vtkGenericDataObjectReader()
reader.SetFileName(filename)
reader.Update()

points = np.array( reader.GetOutput().GetPoints().GetData() )

The variable points is an (N,2) or (N,3) array, where N is the number of points.

ToddP
  • 652
  • 13
  • 18
4

Have you tried using paraview? (http://www.paraview.org/) It can give you a visual idea of what is going on behind the scenes and can output the file in a number of different ways. I would suggest this as I don't have a clue what your data is like. http://www.vtk.org/Wiki/VTK/Examples/Python may also have an example that may fit the bill for you. Personally, I'd have a play with paraview and go from there.

g.stevo
  • 745
  • 3
  • 10
  • 2
    I've used paraview before, or rather an add on to it called VISIT. However, I need to analyse what's in the file and do things like fft etc. and so simply visualizing it is not enough. – Ben Jackel Aug 01 '12 at 13:16
4

It should be mentioned that in its latest release, the yt project http://yt-project.org/ includes support for ATHENA, meaning that by all means this is way to analyze the simulation data using python.

Tim
  • 49
  • 1
  • 1