[Paraview] visualising a numpy array

[Cory Quammen]

Hi Juha,

This is a bit beyond my current numpy-foo. Hopefully I'll get a chance
to look at it in more depth in the next couple months.

If you figure anything else out, I and others on the list would surely
appreciate knowing what you find.

Thanks,
Cory

On Fri, Apr 29, 2016 at 2:19 PM, Cory Quammen <cory.quammen at kitware.com> wrote:
> 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.



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