Contour plot from data in a vtk file using Python
Question: I have a set of data stored in a VTK file which represents a cut through a
domain with scalar point data in an array. I am trying to produce a contour
plot of said scalar to make it look somewhat like the attached picture made
using ParaView. I'd rather stick with the vtk libraries than use something
else, like Matplotlib, as I think they produce better visualisations in
general.![Preview of what I want to
achieve](http://i.stack.imgur.com/4dQcU.jpg).

I have looked at several examples on-line but none of them work for me (no
errors are thrown, all I end up with is an empty render with just the
background), all I have been able to do is a surface plot of the data (e.g.:
[here](http://public.kitware.com/cgi-
bin/cvsweb.cgi/~checkout~/VTK/Examples/ImageProcessing/Python/Contours2D.py)
).

Here's the current version of the code I have (very similar to the one that
successfully produces the surface plot):

    
    
    # import data
    reader = vtk.vtkDataSetReader()
    reader.SetFileName('inputDataFiles/k_zCut.vtk')
    reader.ReadAllVectorsOn()
    reader.ReadAllScalarsOn()
    reader.Update()
    
    # access data
    data = reader.GetOutput()
    d = data.GetPointData()
    array=d.GetArray('k')
    
    # create the filter
    contours = vtk.vtkContourFilter()
    contours.SetInput(reader.GetOutput())
    contours.GenerateValues(5,1.,5.)
    
    # create the mapper
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInput(contours.GetOutput())
    mapper.ScalarVisibilityOff()
    mapper.SetScalarRange(1., 5.)
    
    # create the actor
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    
    # create a rendering window and renderer
    ren = vtk.vtkRenderer()
    ren.SetBackground( 0.329412, 0.34902, 0.427451 ) #Paraview blue
    
    # Assign actor to the renderer
    ren.AddActor(actor)
    
    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(ren)
    renWin.SetSize(750, 750)
    
    # create a renderwindowinteractor
    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)
    
    # render
    renWin.Render()
    
    # screenshot
    w2if = vtk.vtkWindowToImageFilter()
    w2if.SetInput(renWin)
    w2if.Update()
    w2if.SetMagnification(5.)
    
    writer = vtk.vtkPNGWriter()
    writer.SetFileName("screenshot.png")
    writer.SetInput(w2if.GetOutput())
    writer.Write()
    
    # Enable user interface interactor
    iren.Initialize()
    iren.Start()
    

Below you can see a shortened part of my input file. Any help will be much
appreciated.

    
    
    # vtk DataFile Version 2.0
    sampleSurface
    ASCII
    DATASET POLYDATA
    POINTS 34813 float
    0 0 0
    0 -0.000191589 0
    0.000264399 0.000157061 0
    0 0.000313389 0
    0.000264347 -0.000191923 0
    0 -0.000383178 0
    -0.000395709 0 0
    -0.000395709 0.000156695 0
    3.60174e-05 0.000486922 0
    0.000528387 0 0
    
    POLYGONS 69284 277136
    3 4105 4371 3861
    3 4102 3861 4371
    3 4656 4371 4373
    3 4105 4373 4371
    3 3624 3861 3390
    3 3621 3390 3861
    3 4105 3863 3861
    3 3624 3861 3863
    3 3188 3390 2990
    3 3187 2990 3390
    3 3624 3390 3391
    3 3188 3391 3390
    
    POINT_DATA 34813
    FIELD attributes 1
    k 1 34813 float
    0.849464 0.391519 1.52947 1.05206 0.391519 0.253736 1.39481 1.39481 0.636517 1.21019
    0.640193 0.114295 1.12557 0.644143 0.629569 0.114295 0.485032 0.477396 1.39961 0.0860201
    1.66665 1.24058 1.45939 0.483719 1.01318 0.163198 0.317574 0.792821 0.317125 0.658835
    


Answer: If you'd like to use a more "pythonic" interface to VTK, [consider using
`mayavi`/`tvtk`/`mlab`](http://docs.enthought.com/mayavi/mayavi/mlab.html).
(Either way, VTK is an excellent choice for this!)

`tvtk` is a slightly more pythonic, low-level, python binding to VTK with a
handful of really nice features (e.g. transparent usage of numpy arrays).

Mayavi and mlab give a more high-level interface to VTK.

The snippet of the data file that you showed is invalid as-is, but if we use a
similar one:

    
    
    # vtk DataFile Version 2.0
    Simple VTK file example
    ASCII
    
    DATASET POLYDATA
    POINTS 9 float
    3.0 0.0 0.0
    1.0 1.0 0.0
    0.0 3.0 0.0
    3.0 0.0 1.0
    1.0 1.0 1.0
    0.0 3.0 1.0
    3.0 2.0 2.0
    2.0 2.0 2.0
    2.0 3.0 2.0
    
    TRIANGLE_STRIPS 2 14
    6 0 3 1 4 2 5
    6 3 6 4 7 5 8
    
    POINT_DATA 9
    SCALARS nodal float
    LOOKUP_TABLE default
    0.0 0.1 0.0 0.3 0.6 0.3 0.8 1.0 0.8
    

Rendering contours on the surface can be as simple as:

    
    
    from mayavi import mlab
    
    source = mlab.pipeline.open('test.vtk')
    lines = mlab.pipeline.contour_surface(source)
    
    mlab.show()
    

![enter image description here](http://i.stack.imgur.com/H0l6j.png)

Or we can get a bit fancier:

    
    
    from mayavi import mlab
    
    # Make a figure with a black background
    fig = mlab.figure(bgcolor=(0,0,0))
    # Also see methods like: fig.scene.z_plus_view(), etc
    fig.scene.camera.azimuth(215)
    
    source = mlab.pipeline.open('test.vtk')
    # Show the surface, colored by the scalars
    surf = mlab.pipeline.surface(source)
    # Draw contours of the scalars on the surface
    lines = mlab.pipeline.contour_surface(source)
    
    mlab.show()
    

![enter image description here](http://i.stack.imgur.com/FGBhy.png)

