Python script for plotting the evolution of charts such as Paraview does

Question

I want to code a python script that generates a plot like the one shown to the right in the next screenshot from Paraview:

I have a series of files generated with the command foamToVTK:

Is there in VTK any function similar to PlotOverLine method of Paraview?

Answer 1

Solution with vtkplotter:

from vtkplotter import *

img = load(datadir+'vase.vti').imagedata()

# write a test file, return a vtkMultiBlockData
mblock = write([img], "multiblock.vtm")

# load back from file into a list of meshes/volumes
mobjs = load("multiblock.vtm")

p1, p2 = (10,10,10), (80,80,80)
pl = probeLine(mobjs[0], p1, p2).lw(3)

vals = pl.getPointArray()
xvals = pl.points()[:,0]

plt = plotxy([xvals, vals],
             lc="r",       # line color
             marker="*",   # marker style
             mc="dr",      # marker color
             ms=0.6,       # marker color
            )

show([(img, pl), plt], N=2, sharecam=False, axes=1, bg='w')

Answer 2

Inspired by @Nico Vuaille and https://stackoverflow.com/a/21990296/4286662 answers, I came across another different solution based on probe filters. This solution remains faithfully to the Paraview results. This is the working code:

from vtk.util import numpy_support as VN
from matplotlib import pyplot as plt
import vtk
import numpy as np
from os import walk, path, system
import pandas as pd

def getFilenames():
    ### Get the file names of each step of the simulation
    (dirpath, dirnames, filenames) = next(walk('VTK'))
    ids=[]
    for dir in dirnames:
        ids.append(int(dir.split("_")[1]))
    ids = sorted(ids)
    basename = dirnames[0].split("_")[0]

    return ids, basename, dirpath

def createLine(p1, p2):
    # Create the line along which you want to sample
    line = vtk.vtkLineSource()
    line.SetResolution(numPoints)
    line.SetPoint1(p1)
    line.SetPoint2(p2)
    line.Update()
    return line

def probeOverLine(line,reader):
    ### Interpolate the data from the VTK-file on the created line.
    probe = vtk.vtkProbeFilter()
    probe.SetInputConnection(line.GetOutputPort())
    probe.SetSourceConnection(reader.GetOutputPort())
    probe.Update()
    ### Get the velocity from the probe
    return VN.vtk_to_numpy(probe.GetOutput().GetPointData().GetArray('U'))


### Initialization of variables
cnt=1
fig= plt.figure()
numPoints=100

ids, basename, dirpath = getFilenames()

### Iteration of time steps
for id in ids[1:]:
    ### Read values from the file of this time step
    filename = "%s/%s_%d/internal.vtu" % (dirpath, basename, id)
    reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(filename)
    reader.Update()

    if cnt==1:
        ### Get points for a line in the center of the object
        bounds = reader.GetOutput().GetBounds()
        p1 = [(bounds[1] - bounds[0]) / 2.0 + bounds[0], bounds[2], 0]
        p2 = [(bounds[3] - bounds[2]) / 2.0 + bounds[2], bounds[3], 0]
        line = createLine(p1, p2)

    plotOverLine = probeOverLine(line, reader)

    df = pd.DataFrame(plotOverLine, columns=['X', 'Y', 'Z'])
    plt.clf()
    plt.title('Frame %d' % cnt)
    plt.plot(df)
    plt.legend(df.columns)
    axes = plt.gca()
    plt.draw()
    plt.pause(.05)

    cnt+=1

plt.show()

With the following result:

Answer 3

The ParaView PlotOverLine is in vtk as vtkProbeFilter.

Minimal working example:

import vtk

# Original data
source = vtk.vtkRTAnalyticSource()

# the line to plot over
line = vtk.vtkLineSource()

# filter
probeFilter = vtk.vtkProbeFilter()
probeFilter.SetInputConnection(line.GetOutputPort())
probeFilter.SetSourceConnection(source.GetOutputPort())

# rendering
plot = vtk.vtkXYPlotActor()
plot.AddDataSetInputConnection(probeFilter.GetOutputPort())
plot.SetTitle("My plot")
window = vtk.vtkRenderWindow()
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(window)
renderer = vtk.vtkRenderer()
renderer.AddActor2D(plot)
window.AddRenderer(renderer)
window.Render()
interactor.Start()

You can find a more complex (multi plot, colors ...) c++ example here: https://lorensen.github.io/VTKExamples/site/Cxx/Annotation/XYPlot/ The python API is the same.

Answer 4

I figured out a solution to this problem. It is not likely the optimal one, though. For this solution, I firstly converted the mesh to a vtkUnstructuredGrid (in this case using a resolution of 256 points.)

Following I include the code I used to carry this out:

import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy
from os import walk, path, system
import pandas as pd

### Create the VTK files
system("foamToVTK")

### Initialization of variables
cnt=1
fig= plt.figure()
npts = 256  #dimensions of the grid

### Get the file names of each step of the simulation
(dirpath, dirnames, filenames) = next(walk('VTK'))
ids=[]
for dir in dirnames:
    ids.append(int(dir.split("_")[1]))
ids = sorted(ids)
basename = dirnames[0].split("_")[0]

### Iteration of time steps
for id in ids[1:]:

    ### Read values from the file of this time step
    filename = "%s/%s_%d/internal.vtu" % (dirpath, basename, id)
    reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(filename)
    reader.Update()

    ### Get the coordinates of nodes in the mesh using VTK methods
    nodes_vtk_array= reader.GetOutput().GetPoints().GetData()
    vtk_array = reader.GetOutput().GetPointData().GetArray('U') #Velocity (3 dimensions)
    numpy_array = vtk_to_numpy(vtk_array)
    nodes_nummpy_array = vtk_to_numpy(nodes_vtk_array)
    x,y,z= nodes_nummpy_array[:,0] , nodes_nummpy_array[:,1] , nodes_nummpy_array[:,2]
    xmin, xmax = min(x), max(x)
    ymin, ymax = min(y), max(y)

    ### Define grid
    xi = np.linspace(xmin, xmax, npts)
    yi = np.linspace(ymin, ymax, npts)

    ### Grid the data
    interpolated = griddata((x, y), numpy_array, (xi[None,:], yi[:,None]), method='cubic')

    ### Create the list of points
    plotOverLine=[]
    for point in range(len(interpolated[0])):
        plotOverLine.append(interpolated[127][point])

    ### Update and plot the chart for this time step
    df = pd.DataFrame(plotOverLine, columns=['X', 'Y', 'Z'])
    plt.clf()
    plt.title('Frame %d' % cnt)
    plt.plot(df)
    plt.legend(df.columns)
    axes = plt.gca()
    axes.set_ylim([-15,10])
    plt.draw()
    plt.pause(.05)

For each time step, it updates and shows a plot like this: