Xarray-to-vtk

[johnkit]

Using vtkNetCDFCFReader as described in #99

• to netcdf file to paraview
• to netcdf file to Python (PyVista Plotter for display)
• Xarray direct (no netcdf file)
• Xarray slice feature

[danlipsa]

to netcdf file to paraview

File > Open and choose NetCDF Reader
Attached is a paraview state file for opening air_temperature.nc also attached.
air_temperature.pvsm.zip
air_temperature.nc.zip

[danlipsa]

to netcdf file to Python (PyVista Plotter for display)
Xarray slice feature

Build VTK from:
https://gitlab.kitware.com/danlipsa/vtk/-/commits/xarray_quick
then execute the attached script like this:
PYTHONPATH=~/projects/pyvista:~/projects/vtk/xarray/build/lib/python3.10/site-packages/ python test.py

test.py.zip

Note PYTHONPATH points to the source code for pyvista and the location where you built VTK

Use n/N to increase/decrease time, l/L to increase/decrease longitude size

[johnkit]

Interim example:

def update(pl, xr_ds, time_index):
    ds = xr_ds.vtk.dataset(
        active_scalars='air', time_index=time_index,
        encoding={"air": {"dtype": float}})
    print(ds)
    mesh = pv.wrap(ds)
    pl.add_mesh(mesh, name="air")

time_index = 0
lon_remove = 0
pl = pv.Plotter()
pl.add_key_event('n', next_time_step)
pl.add_key_event('N', prev_time_step)
pl.add_key_event('l', increase_lon)
pl.add_key_event('L', decrease_lon)
xr_ds = xr.tutorial.load_dataset("air_temperature")
print(xr_ds)
lon_length=len(xr_ds['lon'])
update(pl, xr_ds, time_index)
pl.show()

[johnkit]

(also in #99)
I put a small zarr example (example.zarr.tgz) at https://data.kitware.com/#item/66d225af9eee4150438e24ae. (Github wouldn't let me upload it here for some reason.) Note that you have to untar before trying to load :)

[danlipsa]

Here is the WIP MR for the full xarray to vtk: https://gitlab.kitware.com/vtk/vtk/-/merge_requests/11513

[aashish24]

@danlipsa here is link to sample xarray datasets: https://github.com/pydata/xarray-data that we should try to load in your PR

[danlipsa]

The attached script works with the xarray MR: https://gitlab.kitware.com/vtk/vtk/-/merge_requests/11513
xarray_vtk.zip

Variables are shallow copied (unless the numpy array needs to be made contiguous: https://numpy.org/devdocs/reference/generated/numpy.ascontiguousarray.html). Coordinates are deep copied at the moment, because the are stored as double array in the netcdf reader so they require a bit more work - they need to be stored as the orginal type in the reader.

[danlipsa]

[danlipsa]

This is the script that produces these images from https://github.com/pydata/xarray-data
xarray_vtk.zip

[jourdain]

@danlipsa are everything available in the nightly VTK wheel?

Also do you have a small Python example using the reader from VTK?

[danlipsa]

Yes, everything is there. I used: https://gitlab.kitware.com/vtk/vtk/-/jobs/10654362
for linux.

I use
https://gitlab.kitware.com/vtk/vtk/-/blob/master/IO/NetCDF/Testing/Python/GetReader.py
to either get the reader from the file directly or zero-copy from the xarray.

In few cases you'll need to:

            ds_tree = xr.open_datatree(args.input)
            ds_xr = ds_tree[args.node].to_dataset()

and then get the dataset from the xarray as before.

[jourdain]

Where did you get the salt dataset?

[danlipsa]

All datasets are from: https://github.com/pydata/xarray-data/
The salt one is ROMS_example.nc

[jourdain]

Here is the API we have in our dataset builder that we would need to map to the VTK API. Some are rectilinear specific but for now I'm just listing everything so we can create an adapter that is more python friendly.

So far, by installing vtk nightly, I was able to load xarray dataset using the accessor. But because PyVista is not compatible with 9.4, I could not easily render the mesh within Jupyter.

Reader requirements

• Rectilinear specific
  • x, x_size, y, y_size, z, z_size (coord name + size) (? GetAllDimensions())
  • slice_extends
  • slices (set/get sub domain selection)
• Generic
  • t (set/get coord name), t_index, t_size, t_labels (list time in a human fashion)
  • available_coords (list 1d xarray coords)
  • arrays (set/get the fields to load/loaded)
  • available_arrays (list available arrays => reader.GetAllVariableArrayNames())

Questions: For the VTK reader how do I ...

1. Handle time? Use pipeline information pass?
2. Select fields to load? I can figure out, what is available, but how can I activate/deactivate specific field?
3. For compatible format, how can I query/set a subset selection (slice, range)

Next step:

Once we understand how do some of those operations, we should create a Python adapter or improve the C++ code to help that integration process.

Comments:

1. The xarray accessor should just do ds.vtk.reader [no ()] by defining the method as a property. Will need to fix
2. The helper get_accessor should either be internal or a @property
3. We may want another Python class that decorate an xarray into a vtk mesh algo/reader.

[danlipsa]

Questions: For the VTK reader how do I ...

1. Handle time? Use pipeline information pass?
2. Select fields to load? I can figure out, what is available, but how can I activate/deactivate specific field?
3. For compatible format, how can I query/set a subset selection (slice, range)

The reader returned is a vtkNetCDFCFReader so:

1. Indeed time is handled using the pipeline.
2. Many readers have this: https://vtk.org/doc/nightly/html/classvtkNetCDFReader.html#ab457b03b7de7908dd6437b14fcfc7538
3. vtkExtractGrid

[jourdain]

Some weird behavior that I'm noticing... When loading an XArray without enabling any field, the output is a vtkImageData. Once I activate a field, it became a vtkStructuredGrid.

I think the switch of dataset type is messing up the VTK pipeline. I'm going to talk with Berk about that and may implement a work around into our viewers/explorers.

[danlipsa]

Try to load ROMS_example.nc in paraview. You'll see that you can select Dimensions. (you should choose (s_rho, eta_rho, xi_rho)). The list you see is AllDimensions you asked about earlier: the list of all possible set of dimensions for all variables.

Once you select a dimension, all variables that have that dimension are loaded. Probably you'll need to do something similar in trame.

[danlipsa]

For different variables, you can have different types of output. If you don't specify and variable I assume it just uses the default which is image data.