Note

This page was generated from a jupyter notebook.

Writing with legacy vtk files

(GE Tucker, CU Boulder, first version June 2023)

The Visualization Toolkit (VTK) is an open-source software system for visualization. It offers two file formats: one based on XML, and the other—the so-called legacy VTK format—using a simple text-based format. These file formats are readable by visualization packages such as ParaView, so it is useful to be able to output VTK format for Landlab grids and fields. This tutorial demonstrates Landlab’s legacy VTK file-output capability.

The format

To quote from the VTK documentation:

“The legacy VTK file formats consist of five basic parts.

  1. The first part is the file version and identifier. This part contains the single line: vtk DataFile Version x.x. This line must be exactly as shown with the exception of the version number x.x, which will vary with different releases of VTK. (Note: the current version number is 3.0. Version 1.0 and 2.0 files are compatible with version 3.0 files.)

  2. The second part is the header. The header consists of a character string terminated by end-of-line character :nbsphinx-math:`n`. The header is 256 characters maximum. The header can be used to describe the data and include any other pertinent information.

  3. The next part is the file format. The file format describes the type of file, either ASCII or binary. On this line the single word ASCII or BINARY must appear.

  4. The fourth part is the dataset structure. The geometry part describes the geometry and topology of the dataset. This part begins with a line containing the keyword DATASET followed by a keyword describing the type of dataset.Then, depending upon the type of dataset, other keyword/data combinations define the actual data.

  5. The final part describes the dataset attributes. This part begins with the keywords POINT_DATA or CELL_DATA, followed by an integer number specifying the number of points or cells, respectively. (It doesn’t matter whether POINT_DATA or CELL_DATA comes first.) Other keyword/data combinations then define the actual dataset attribute values (i.e., scalars, vectors, tensors, normals, texture coordinates, or field data).”

Hex grid example

This example creates and outputs a tiny hex grid, along with two fields: topographic__elevation and surface_water__depth.

[ ]:

import numpy as np

import landlab.io.legacy_vtk as vtk
from landlab import HexModelGrid, RasterModelGrid

[ ]:

# Create a tiny grid with 1 core node and 6 boundary nodes
grid = HexModelGrid((3, 2))

# Add two fields with made-up values
topo = grid.add_zeros("topographic__elevation", at="node")
topo[3] = 1.0
grid.at_node["surface_water__depth"] = np.arange(grid.number_of_nodes)

# Write output in legacy VTK format
vtk_file = vtk.dump(grid)

Let’s see what the output looks like:

[ ]:

print(vtk_file)

Raster grid example

[ ]:

# Create a tiny grid with 1 core node and 6 boundary nodes
grid = RasterModelGrid((3, 3))

# Add two fields with made-up values
topo = grid.add_zeros("topographic__elevation", at="node")
topo[4] = 1.0
grid.at_node["surface_water__depth"] = np.arange(grid.number_of_nodes)

# Write output in legacy VTK format
vtk_file = vtk.dump(grid)

[ ]:

print(vtk_file)

Notice that this has saved the grid’s nodes and patches (VTK uses the terms points and cells). If you would like to save the dual grid, you can do this through the at keyword. The default (at="node") is to save the main grid but you can use at="corner" to save the dual grid.

[ ]:

print(vtk.dump(grid, at="corner"))

If your grid contains many fields, you may not want to save all of them. You can specify which fields to save through the include and exclude keywords. These operate much like UNIX filename pattern matching and are the same as described in the fields method of your grid. As an example, the following with exclude any fields that contain the string “surface_water”.

[ ]:

print(vtk.dump(grid, exclude="*surface_water*"))

The VTK format assumes points (i.e. either nodes or corners) are defined by x, y, and z coordinates. The default is to assign a value of 0 to all z coordinates. You can, however, change this behavior and assign your own data. For example, the following code uses the topo array for the z coordinate of each of the nodes.

[ ]:

print(vtk.dump(grid, z_coord=topo))

Generated by nbsphinx from a Jupyter notebook.