Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Contour plot from data in a vtk file using Python

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.

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 ).

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
like image 479
Artur Avatar asked Feb 01 '14 15:02

Artur


People also ask

How do you plot a contour plot?

MatPlotLib with Python Contour plots (sometimes called Level Plots) are a way to show a three-dimensional surface on a two-dimensional plane. It graphs two predictor variables X Y on the y-axis and a response variable Z as contours. These contours are sometimes called the z-slices or the iso-response values.

What is contour plot in Ansys?

A contour plot allows you to visualize three-dimensional data in a two-dimensional plot. You insert a contour plot by selecting Contour Plot in the Traces group. The keyboard shortcut is CTRL+5. You cannot switch between a contour plot and a 3D plot. The contour plot region is shown below.


1 Answers

If you'd like to use a more "pythonic" interface to VTK, consider using mayavi/tvtk/mlab. (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

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

like image 53
Joe Kington Avatar answered Oct 07 '22 01:10

Joe Kington