[Paraview] visualising a numpy array

[Cory Quammen]

On Wed, Apr 27, 2016 at 11:39 AM, Juha Jäykkä <juha.jaykka at gmail.com> wrote:
>> How are you getting the PETSc._DMDA_Vec_array into ParaView? The
>> Programmable Filter is also available in ParaView...
>
> I was going to use
>
> my_numpy_data = petsc_dmda_array[:]
> vtkobject = paraview.numpy_support.numpy_to_vtk(my_numpy_data.copy())
>
> That works fine, but then I got stuck with how to stick this into paraview's
> pipeline.
>
>> You'll need to set up a vtkImageData as the output of a Programmable
>> Source/Filter at the very least with the proper dimensions, then use
>> the vtk.numpy_interface.dataset_adapter tools to adapt your numpy
>> array as a VTK data array that you add to the vtkImageData's point
>> data.
>
> I would have thought I need a vtkStructuredGrid or vtkRectilinearGrid, but I
> don't seem able to figure out how to use them.

If your grid points are regularly spaced in each dimension, which I
think they are in this case, you can use a vtkImageData. See the
ParaViewGuide (http://www.paraview.org/paraview-guide/) pages 33 and
34 for pictorial representations of where vtkStructuredGrid or
vtkRectilinearGrid are needed.

> With vSG I can
>
> pts = vtk.vtkPoints()
> pdo = self.GetOutput()
> newArray = vtk.vtkDoubleArray()
> # this is how far I've got before...
> newArray.SetName("testdata")
> newArray.SetNumberOfComponents(1)
> pdo.GetPointData().AddArray(newArray)
>
> # stuff our numbers into pdo.GetPointData().GetArray(0).InsertNextValue(float)
> one by one,
> # followed by pts.InsertNextPoint(x,y,z), so need to get xyz as well
> xmin,xmax,Nx,ymin,ymax,Ny,zmin,zmax,Nz = -1,1,10,-1,2,20,-1,3,30
> loci =
> numpy.mgrid[xmin:xmax:1j*Nx,ymin:ymax:1j*Ny,zmin:zmax:1j*Nz].reshape(3,-1)
> values = numpy.random.random(loci[0].shape)
> for ind in range(loci[0].shape[0]):
>     pdo.GetPointData().GetArray(0).InsertNextValue(values[ind])
>     pts.InsertNextPoint(loci[0,ind],loci[1,ind],loci[2,ind])
> pdo.SetPoints(pts)
> exts = [xmin,xmax,ymin,ymax,zmin,zmax]
> pdo.SetExtent(exts)

Note that the extents are in the integer index space (e.g., i-j-k) of
the grid points rather than physical space, so exts should be
something like

[0, 9, 0, 19, 0, 29]

Also note that xmax, ymax, and zmax should be the highest possible
index, not the number of grid positions in that dimension.

> which seems to do the right thing except I can not make paraview plot anything
> using this. In the ProgrammableSource's Information tab I do see the right
> bounds and extents (but I am a bit confused about what the "dimension" means
> in extents: it says 3, 4, 5 for x, y, z, which makes me think it has used unit
> lattice spacing even though I would have thought it uses whatever it finds in
> pts. The "Data Arrays" box has "testdata" in it and has the correct value
> range. Also the number of points is correct (6000). I don't know where it gets
> its 24 cells from, though. (24 = 2*3*4, so I am a bit concerned it gets this
> implicitly from bounds again with unit spacing).
>
> But I don't understand what I'm doing wrong as I'm basically modifying
> examples from the wiki.
>
> Using vtkImageData I can only get a Number of Points: 60 (and no plot) and
> using vtkRectilinearGrid I get the same. I can fix the 60 -> 6000 by fixing
> the extents, but that still leaves the bounds and no plot on screen.
>
>> Parallel is tougher... You can do the full data set in serial and
>> distribute using the D3 filter, but obviously that has some downsides.
>
> Unacceptable ones: I cannot even imagine what that bottleneck would do to the
> code.
>
>> I believe it is possible to distribute the dataset using MPI in a
>> Programmable Source/Filter, but have never done so myself.
>
> Hmm... so where does the ProgrammableSource run? I was planning on running the
> whole code in parallel, but didn't really think about how to do that yet. Now
> that I start thinking about it, Catalyst is probably the way to do it, isn't
> it?

The Programmable Source is a VTK filter that runs on all nodes on
which you are running the ParaView server. Catalyst may indeed be
useful for you.

Cory

> Cheers,
> Juha
>
>> On Tue, Apr 26, 2016 at 8:49 AM, Juha Jäykkä <juha.jaykka at gmail.com> wrote:
>> > Dear list,
>> >
>> > What's the current best practice to visualise a (3D) numpy array with
>> > paraview? In particular, I'm wondering if
>> >
>> > vtkdata = numpy_support.numpy_to_vtk(num_array=NumPy_data.ravel(),
>> > deep=True, array_type=vtk.VTK_FLOAT)
>> >
>> > can somehow be directly visualised? In parallel.
>> >
>> > I'm aware of ProgrammableSource() and how to do it using that, but only
>> > from a file! What I'd very much like to do is avoid the file (if I cannot
>> > avoid the file, I'll just do xdmf and hdf5.)
>> >
>> > Just in case it matters, the numpy array is actually not really a numpy
>> > array, but petsc4py.PETSc._DMDA_Vec_array which of course can be "cast"
>> > as numpy.ndarray.
>> >
>> > Best regards,
>> > Juha
>> >
>> > _______________________________________________
>> > Powered by www.kitware.com
>> >
>> > Visit other Kitware open-source projects at
>> > http://www.kitware.com/opensource/opensource.html
>> >
>> > Please keep messages on-topic and check the ParaView Wiki at:
>> > http://paraview.org/Wiki/ParaView
>> >
>> > Search the list archives at: http://markmail.org/search/?q=ParaView
>> >
>> > Follow this link to subscribe/unsubscribe:
>> > http://public.kitware.com/mailman/listinfo/paraview
>



-- 
Cory Quammen
R&D Engineer
Kitware, Inc.