MPAS-Tools has kind of a hodge-podge of unrelated visualization tools. This is largely because MPAS-Tools is conceived of primarily as a repository for tools for manipulating MPAS meshes during their construction and initialization.
ParaView VTK Extractor¶
Data on MPAS meshes can be viewed directly in ParaView. However, those of us who work with output from MPAS components on a regular basis have found that the built-in capability has a number of significant limitations. These include:
ParaView visualizes data on MPAS cell centers by connecting them to form triangles. The data is linearly interpolated within each triangle and geometry associated with the parts of MPAS cells that are adjacent to boundaries (e.g. land in the ocean and sea-ice components) are not visualized.
MPAS output may not contain the mesh information needed by ParaView. This is typical of monthly output files, which do not include the mesh data for more efficient file storage.
MPAS output files may contain significantly more data that is desired for visualization. The data may contain full 3D fields, whereas a few slices at different depths may suffice for visualization. There may also be many more variables in the output than are of interest. This may be a problem if the datasets are quite large and therefore prohibitive to transfer to another machine where visualization may be more practical.
For these reasons and more, we created a Python tool for extracting data on
MPAS meshes in files in VTK format that can be read in by
ParaView. The VTK files provide the data as “cell”, rather than “point” data in
ParaVeiw, meaning each cell has the correct polygonal geometry and will be
visualized with a single, constant value for a given field. The extractor,
which is available either through the
mpas_tools.viz.paraview_extractor.extract_vtk() function or the
paraview_vtk_field_extractor.py command-line interface, can also be used to
visualize data on MPAS vertices (visualized on triangles) and edges (visualized
The extractor also allows the user to select a subset of the indices along
dimensions other than
can be useful for selecting, for example, only the top layer of a 3D field.
This capability can also be used to ignore fields with a given dimension by
selecting no indices along that dimension.
By default, time-independent fields on cells are written to a file
and time-dependent fields on cells are written to a series of files
vtk_files/timeDependentFieldsOnCells.pvd vtk_files/time_series/timeDependentFieldsOnCells.0.vtp vtk_files/time_series/timeDependentFieldsOnCells.1.vtp ...
and similarly for edges and vertices. Time-independent fields can be
included in each time step of the time-dependent fields for with
combine=True. This allows time-dependent and -independent fields
to be combined in filters within Paraview at the expense of considerable
additional storage space.
Since the extractor is fairly complex with a lot of possible use cases, we will describe its functionality in several examples.
Extracting a Temperature Time-series¶
To extract fields across multiple files, passing in a regular expression
for the filename pattern, for example
filename_pattern="hist.mpaso.timeSeriesStatsMonthly.*.nc". Since the mesh
data is often excluded from time series output, you may need to supply a
different file with mesh information, such as an initial condition or a restart
In the following example, we extract the temperature field from a time series of daily output from the MPAS-Ocean model:
from mpas_tools.viz.paraview_extractor import extract_vtk extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list='timeDaily_avg_activeTracers_temperature', dimension_list=['nVertLevels=0'], mesh_filename='init.nc', out_dir='vtk_files', xtime='xtime_startDaily')
The top index (
nVertLevels=0) is extracted from the temperature time series.
The output directory is
'vtk_files'. MPAS output typically includes time
information in string format in a variable called
xtime. In this instance,
we have time-averaged data and the output instead includes a start time
xtime_startDaily that we need to supply instead. (We could also supply the
xtime_endDaily if that were preferred for some reason.)
The result is:
vtk_files/fieldsOnCells.pvd vtk_files/time_series/fieldsOnCells.0.vtp ...
These files can be opened in ParaView with:
$ paraview vtk_files/fieldsOnCells.pvd
This will open all of the files in the
vtk_files/time_series directory as
they are requested within ParaView. See ParaView Macros for some tips
on filling in the continents and displaying the current date in ParaView.
The same extraction could be accomplished with the command-line tool as follows:
$ paraview_vtk_field_extractor.py \ -f "analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc" \ -v timeDaily_avg_activeTracers_temperature -d nVertLevels=0 \ -m init.nc -o vtk_files --xtime=xtime_startDaily
Extracting Multiple Fields¶
In this next example, we will extract
areaCell in addition to temperature.
First, we will extract it into a separate VTK file for time-independent
variables, then we will demonstrate combining it with the time-dependent data.
In the first instance, we add
areaCell to the
combine=False, the default, to indicate that we want to
keep time-independent and time-dependent variables separate:
from mpas_tools.viz.paraview_extractor import extract_vtk extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list=['timeDaily_avg_activeTracers_temperature', 'areaCell'], dimension_list=['nVertLevels=0'], mesh_filename='init.nc', combine=False, out_dir='vtk_files2', xtime='xtime_startDaily')
The result is:
vtk_files2/staticFieldsOnCells.vtp vtk_files2/timeDependentFieldsOnCells.pvd vtk_files2/time_series/timeDependentFieldsOnCells.0.vtp ...
We can open both
vtk_files/timeDependentFieldsOnCells.pvd in the same ParaVeiw sesson and
plot them as we like. But we cannot perform calculations involving both
temperature and cell area very easily. If this were a necessity, it might be
convenient to combine them into the same files:
from mpas_tools.viz.paraview_extractor import extract_vtk extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list=['timeDaily_avg_activeTracers_temperature', 'areaCell'], dimension_list=['nVertLevels=0'], mesh_filename='init.nc', combine=True, out_dir='vtk_files3', xtime='xtime_startDaily')
Now, the result is that
areaCell is included in the time-series files
vtk_files3/fieldsOnCells.pvd vtk_files3/time_series/fieldsOnCells.0.vtp ...
Extracting “All” Fields¶
Sometimes, you want all of the fields from your input files to be extracted.
For this purpose there is a special “variable” called
'all' that gets
translated into the full list of available variables. More often, you want to
extract all the variables on cells, edges or vertices, so there are special
“variables” for this, too:
'allOnVertices'. By default, only variables from the files found by
filename_pattern are expanded by these special variables. If you also want
to include variables from the mesh file, you need to specify
The following example extracts all the variables on cells for both the
time-series and the mesh data. It specifies
maxEdges= so that variables
cellsOnCell) that include this dimension are
from mpas_tools.viz.paraview_extractor import extract_vtk extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list=['allOnCells'], dimension_list=['nVertLevels=0', 'nVertLevelsP1=0', 'maxEdges='], mesh_filename='init.nc', combine=True, include_mesh_vars=True, out_dir='vtk_files4', xtime='xtime_startDaily')
In the previous examples, we saw a basic example of indexing the “extra”
dimensions (i.e. dimensions other than
nEdges) from MPAS output. Here, we show some slightly more involved
Indices for extra dimensions can either be supplied at runtime at a prompt (if
dimension_list=None) or via a list of strings with the dimensions and
associated indices. For each extra dimension, you can specify nothing for the
indices (an empty string, meaning skip any fields with this dimension), a single
index, a comma-separated list of indices, or a range of indices (separated by 1
or 2 colons). For example,
dimension_list=['maxEdges=', 'nVertLeves=0:10:2', 'nParticles=0,2,4,6,8']
will ignore any fields with dimension
maxEdges, extract every other
layer from the first 10 vertical levels (each into its own field) and
extract the 5 specified particles.
An index array can also be specified in this way (and these can be mixed with integer indices in a comma-separated list but not in a colon-separated range):
This will extract fields from the first vertical level and the vertical level
with index given by
maxLevelCell (the deepest layer in each ocean column).
Here is a more complete example that extracts the temperature, salinity and layer thickness at the sea surface and seafloor.
from mpas_tools.viz.paraview_extractor import extract_vtk extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list=['timeDaily_avg_activeTracers_temperature', 'timeDaily_avg_activeTracers_salinity', 'timeDaily_avg_layerThickness'], dimension_list=['nVertLevels=0,maxLevelCell'], mesh_filename='init.nc', out_dir='vtk_files5', xtime='xtime_startDaily')
The resulting fields are named:
timeDaily_avg_activeTracers_temperature_0 timeDaily_avg_activeTracers_temperature_maxLevelCell timeDaily_avg_activeTracers_salinity_0 timeDaily_avg_activeTracers_salinity_maxLevelCell timeDaily_avg_layerThickness_0 timeDaily_avg_layerThickness_maxLevelCell
Time can also be indexed like the other dimensions, but it is not passed to the
dimension_list argument but instead to the
time argument. The time
index string can have any of the following formats:
''- no times are to be extracted (probably not useful for
'n'- the index n is to be extracted
'm,n,p'- the list of indices is to be extracted
'm:n'- all indices from m to n are to be extracted (including m but excluding n, in the typical python indexing convention)
'm:n:s'- all indices from m to n are to be extracted (including m but excluding n, in the typical python indexing convention) with stride s between indices
In this example, we extract every 6 days from the daily data set starting with
the beginning fo the data set and continuing to the end (by not specifying the
from mpas_tools.viz.paraview_extractor import extract_vtk extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list=['timeDaily_avg_activeTracers_temperature'], dimension_list=['nVertLevels=0'], mesh_filename='init.nc', time='0::6', out_dir='vtk_files6', xtime='xtime_startDaily')
Some MPAS files, for example mesh files and initial conditions, contain a
Time dimension but no
xtime variable. The extractor will complain about
this unless you specify
ignore_time=True. In this case, only the first time
index is used and all fields are considered to be time-independent, ending
The extractor can produce files in lon/lat coordinates instead of 3D Cartesian
lonlat=True. Polygons near the prime meridian (0 or 360 degrees
longitude) will end up on one side or the other based on the location of the
cell center. This leads to a “ragged” edge a the prime meridian, particularly
for coarse-resolution meshes:
Below, we will provide a method for handling this issue in ParaView.
Here, we extract the temperature field as in Extracting a Temperature Time-series, but this time in lon/lat coordinates.
from mpas_tools.viz.paraview_extractor import extract_vtk extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list='timeDaily_avg_activeTracers_temperature', dimension_list=['nVertLevels=0'], mesh_filename='init.nc', lonlat=True, out_dir='vtk_files7', xtime='xtime_startDaily')
In ParaView, the data lies approximately 0 and 360 degrees long the x axis but with some polygons extending partially beyond these bounds. Since we typically wish to see the data between -180 and 180 degrees longitude, the proposed fix will take care of both the longitude range and the ragged edges in one go.
First, we make a duplicate copy of the data, translated by -360 degrees. Open
a Transform Filter in ParaView and enter
-360 in the first Translate cell
(the x axis). Uncheck “Show Box” and hit “Apply”. The original data will
disappear but you can simply click the eye icon next to it to make it reappear.
Next, we want to combine the original and translated versions of the mesh into a single dataset. Use the shift key to select both, then open the Group Datasets Filter, then hit “Apply”.
Finally, we will crop the grouped dataset to the range of -180 to 180 degrees:
Select the final of the three “GroupDataset1” items,
open a Clip Filter,
select “Plane” as the Clip Type,
set the Origin to -180, 0, 0,
set the Normal to -1, 0, 0
This has clipped the extra data off the left edge. Now for the right edge:
Select the “Clip1” item,
open a Clip Filter again,
select “Plane” as the Clip Type,
set the Origin to 180, 0, 0,
set the Normal to 1, 0, 0,
uncheck “Show Plane”
Now, the data should have clean edges. If you want, you can put a plane behind it to fill in the land:
The extractor includes optional support for extracting geometry appropriate
for displaying variables at the depth of a topographic feature (typically
the top or bottom of the domain) for MPAS components with a spatially
variable top or bottom index (e.g.
maxLevelCell in MPAS-Ocean). This is
accomplished with arguments such as:
Fields on cells are sampled at the topographic index and the geometry
includes polygons corresponding to edges so that vertical faces between
adjacent cells can be displayed. Fields are extracted as normal except
that they are sampled as point data rather than cell data, allowing
computations in ParaView to display the topography. A mask field is also
included, indicating which parts of edge polygons correspond to the boundary
of the domain (
boundaryMask == 1) and which parts of cell and edge
polygons are interior (
boundaryMask == 0).
In the following, we make sure to use
include_mesh_vars=True because we need
bottomDepth from the mesh file
to be included in the time-dependent output files. We are not interested in
variables with dimensions
maxEdges so we remove those
dimensions by leaving their index strings blank.
from mpas_tools.viz.paraview_extractor import extract_vtk extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list='allOnCells', dimension_list=['nVertLevelsP1=', 'maxEdges='], topo_dim='nVertLevels', topo_cell_index='maxLevelCell', combine=True, include_mesh_vars=True, mesh_filename='init.nc', out_dir='vtk_files8', xtime='xtime_startDaily')
Together, this can be used to plot topography by using a Calculator Filter in ParaView, checking the “Coordinate Result” box, and entering the following:
coords*(1.0 + 100.0/mag(coords)*((1 - boundaryMask)*(-bottomDepth) + 10.0*boundaryMask))
The result is that the MPAS-Ocean topography is displayed with a vertical exaggeration of 100 and with a value equivalent to 10 m along boundary points of edge polygons (a “water-tight” surface).
Here is what that looks like for a 240-km (very coarse) ocean mesh:
The same approach can be used with
lonlat=True. In this case, the
Calculator Filter is a bit simpler:
coords + 1e-3*(1 - boundaryMask)*(-bottomDepth)*kHat + 1.0*boundaryMask*kHat
Here is the bottom temperature in such a plot:
Extracting a Region¶
Some simulations are focused on a small region, even though the entire globe is
included in the mesh. For such situations, we provide a way to extract a
subset of the data over a region before converting it to VTK format. The user
geometric_features package as an argument. The regions in this
feature collection are used to define a mask, and the MPAS data is culled to
lie within the mask before conversion to VTK proceeds.
The region should indicate the parts of the mesh to keep, not those to remove.
In this example, we extract sea surface temperature only in the Southern Ocean:
from mpas_tools.viz.paraview_extractor import extract_vtk from geometric_features import GeometricFeatures gf = GeometricFeatures() fc = gf.read(componentName='ocean', objectType='region', featureNames=['Southern Ocean']) extract_vtk( filename_pattern='analysis_members/mpaso.hist.am.timeSeriesStatsDaily.*.nc', variable_list='timeDaily_avg_activeTracers_temperature', dimension_list=['nVertLevels=0'], mesh_filename='init.nc', fc_region_mask=fc, out_dir='vtk_files9', xtime='xtime_startDaily')
The first of these adds a sphere that is just a bit smaller than the MPAS data
on the sphere so that continents are not holes in the data. The second can
be used to display the current time (extracted from the
xtime variable) in
a ParaVeiw animation.
MPAS Mesh to Triangles¶
A relatively new and under-explored functionality in MPAS-Tools is the
capability to extract a triangle mesh for visualization from an MPAS mesh.
This functionality comes from the function
mpas_tools.viz.mesh_to_triangles.mesh_to_triangles(). The function
takes an MPAS mesh as an
xarray.Dataset object as its only required input
and produces another
xarray.Dataset with the triangle mesh that connects
pairs of adjacent vertices to cell centers as its output. Thus, each hexagon
becomes 6 triangles, each pentagon becomes 5, and so on.
In addition to the points and connectivity data for defining these trinagles,
the output dataset,
dsTris, also includes the cell index that each triangle
is in and cell indices and weights for interpolating data defined at cell
centers to triangle nodes.
dsTris includes variables
the cell that each triangle is part of;
nodeCellWeights, the indices and weights used to interpolate from MPAS cell
centers to triangle nodes; Cartesian coordinates
latNode in radians.
guaranteed to be within 180 degrees of the cell center corresponding to
triCellIndices. Nodes always have a counterclockwise winding.
Here is an example workflow for using this function:
import xarray import numpy import matplotlib.pyplot as plt from matplotlib.tri import Triangulation from mpas_tools.viz import mesh_to_triangles dsMesh = xarray.open_dataset('mpaso.rst.0501-01-01_00000.nc') dsTris = mesh_to_triangles(dsMesh, periodicCopy=True) sst = dsMesh.temperature.isel(Time=0, nVertLevels=0).values sstTri = sst[dsTris.triCellIndices] sstNode = (sst[dsTris.nodeCellIndices]*dsTris.nodeCellWeights).sum('nInterp') nTriangles = dsTris.sizes['nTriangles'] tris = numpy.arange(3*nTriangles).reshape(nTriangles, 3) lonNode = numpy.rad2deg(dsTris.lonNode.values).ravel() latNode = numpy.rad2deg(dsTris.latNode.values).ravel() sstNode = sstNode.values.ravel() triangles = Triangulation(lonNode, latNode, tris) plt.figure(1) plt.tripcolor(triangles, sstNode, shading='gouraud') plt.xlim([0., 360.]) plt.ylim([-90., 90.]) plt.figure(2) plt.tripcolor(triangles, sstTri, shading='flat') plt.xlim([0., 360.]) plt.ylim([-90., 90.]) plt.show()
In this example,
mpaso.rst.0501-01-01_00000.nc is a restart file from a
simulation with an EC at the default 30 to 60 km resolution (see
Defining an Eddy-closure Mesh); the restart file contains both mesh information and a snapshot
of the 3D temperature field.
Here are the resulting plots (which look nearly identical at this resolution):
Some MPAS-Tools routines include plots of mesh resolution and related variables.
We have found it handy to use the
SciVisColor Colormaps for some of these plots.
Unfortunately, there is not a standard python package for adding these
matplotlib (as is the case for
cmocean, for example). To add the
SciVisColor colormaps, call the function
mpas_tools.viz.colormaps.register_sci_viz_colormaps(). No arguments
are required, as the XML files for defining the colormaps are included in the
In this example, we use the