Mesh Conversion

Mesh Converter

The command-line tool MpasMeshConverter.x and the corresponding wrapper function mpas_tools.mesh.conversion.convert() are used to convert a dataset describing cell locations, vertex locations, and connectivity between cells and vertices into a valid MPAS mesh following the MPAS mesh specification.

Example call to the command-line tool:

$ planar_hex --nx 4 --ny 4 --dc 10e3 -o base_mesh.nc
$ MpasMeshConverter.x base_mesh.nc mesh.nc

This example uses the planar_hex tool to generate a small, doubly periodic MPAS mesh with 10-km cells, then uses the MPAS mesh converter to make sure the resulting mesh adheres to the mesh specification. MpasMeshConverter.x takes two arguments, the input and output mesh files, and will prompt the user for file names if these arguments are not supplied.

The same example in a python script can be accomplished with:

from mpas_tools.planar_hex import make_planar_hex_mesh
from mpas_tools.mesh.conversion import convert
from mpas_tools.io import write_netcdf

ds = make_planar_hex_mesh(nx=4, ny=4, dc=10e3, nonperiodic_x=False,
                          nonperiodic_y=False)
ds = convert(ds)
write_netcdf(ds, 'mesh.nc')

Regardless of which of these methods is used, the input mesh must define the following dimensions, variables and global attributes (not the dimension sizes are merely examples from the mesh generated in the previous examples):

netcdf mesh {
dimensions:
  nCells = 16 ;
  nVertices = 32 ;
  vertexDegree = 3 ;
variables:
  double xCell(nCells) ;
  double yCell(nCells) ;
  double zCell(nCells) ;
  double xVertex(nVertices) ;
  double yVertex(nVertices) ;
  double zVertex(nVertices) ;
  int cellsOnVertex(nVertices, vertexDegree) ;
  double meshDensity(nCells) ;

// global attributes:
      :on_a_sphere = "NO" ;
      :sphere_radius = 0. ;
      :is_periodic = "YES" ;

The variable meshDensity is required for historical reasons and is passed unchanged to the resulting mesh. It can contain all ones (or zeros), the resolution of the mesh in kilometers, or whatever field the user wishes.

The following global attributes are optional but will be passed on to the resulting mesh:

// global attributes:
      :x_period = 40000. ;
      :y_period = 34641.0161513775 ;
      :history = "Tue May 26 20:58:10 2020: /home/xylar/miniconda3/envs/mpas/bin/planar_hex --nx 4 --ny 4 --dc 10e3 -o base_mesh.nc" ;

The file_id global attribute is also optional and is preserved in the resulting mesh as parent_id. A new file_id (a random hash tag) is generated by the mesh converter.

The resulting dataset has all the dimensions and variables required for an MPAS mesh.

The mesh converter also generates a file called graph.info that is used to create graph partitions from tools like Metis. This file is not stored by default when the python cull function is used but can be specified with the graphInfoFileName argument. The python function also takes a logger that can be used to capture the output that would otherwise go to the screen via stdout and stderr.

Cell Culler

The command-line tool MpasCellCuller.x and the corresponding wrapper function mpas_tools.mesh.conversion.cull() are used to cull cells from a mesh based on the cullCell field in the input dataset and/or the provided masks. The contents of the cullCell field is merge with the mask(s) from a masking dataset and the inverse of the mask(s) from an inverse-masking dataset. Then, a preserve-masking dataset is used to determine where cells should not be culled.

Example call to the command-line tool, assuming you start with a spherical mesh called base_mesh.nc (not the doubly periodic planar mesh from the examples above):

$ merge_features -c natural_earth -b region -n "Land Coverage" -o land.geojson
$ MpasMaskCreator.x base_mesh.nc land.nc -f land.geojson
$ MpasCellCuller.x base_mesh.nc culled_mesh.nc -m land.nc

This example uses the merge_features tool from the geometric_features conda package to get a geojson file containing land coverage. Then, it uses the mask creator (see the next section) to create a mask on the MPAS mesh that is one inside this region and zero outside. Finally, it culls the base mesh to only those cells where the mask is zero (i.e. the mask indicates which cells are to be removed).

The same example in a python script can be accomplished with:

import xarray
from geometric_features import GeometricFeatures
from mpas_tools.mesh.conversion import mask, cull

gf = GeometricFeatures()

fcLandCoverage = gf.read(componentName='natural_earth', objectType='region',
                         featureNames=['Land Coverage'])

dsBaseMesh = xarray.open_dataset('base_mesh.nc')
dsLandMask = mask(dsBaseMesh, fcMask=fcLandCoverage)
dsCulledMesh = conversion.cull(dsBaseMesh, dsMask=dsLandMask)
write_netcdf(dsCulledMesh, 'culled_mesh.nc')

Here is the full usage of MpasCellCuller.x:

MpasCellCuller.x [input_name] [output_name] [[-m/-i/-p] masks_name] [-c]

    input_name:
        This argument specifies the input MPAS mesh.
    output_name:
        This argument specifies the output culled MPAS mesh.
        If not specified, it defaults to culled_mesh.nc, but
        it is required if additional arguments are specified.
    -m/-i/-p:
        These arguments control how a set of masks is used when
        culling a mesh.
        The -m argument applies a mask to cull based on (i.e.
        where the mask is 1, the mesh will be culled).
        The -i argument applies the inverse mask to cull based
        on (i.e. where the mask is 0, the mesh will be
        culled).
        The -p argument forces any marked cells to not be
        culled.
        If this argument is specified, the masks_name argument
        is required
    -c:
        Output the mapping from old to new mesh (cellMap) in
            cellMapForward.txt,
        and output the reverse mapping from new to old mesh in
            cellMapBackward.txt.

Mask Creator

The command-line tool MpasMaskCreator.x and the corresponding wrapper function mpas_tools.mesh.conversion.mask() are used to create a set of region masks either from mask features or from seed points to be used to flood fill a contiguous block of cells.

Examples usage of the mask creator can be found above under the Cell Culler.

Here is the full usage of MpasMaskCreator.x:

MpasMaskCreator.x in_file out_file [ [-f/-s] file.geojson ] [--positive_lon]
    in_file: This argument defines the input file that masks will be created for.
    out_file: This argument defines the file that masks will be written to.
    -s file.geojson: This argument pair defines a set of points (from the geojson point definition)
        that will be used as seed points in a flood fill algorithim. This is useful when trying to remove isolated cells from a mesh.
    -f file.geojson: This argument pair defines a set of geojson features (regions, transects, or points)
        that will be converted into masks / lists.
    --positive_lon: It is unlikely that you want this argument.  In rare cases when using a non-standard geojson
        file where the logitude ranges from 0 to 360 degrees (with the prime meridian at 0 degrees), use this flag.
        If this flag is not set, the logitude range is -180-180 with 0 degrees being the prime meridian, which is the
        case for standar geojson files including all features from the geometric_feature repo.
        The fact that longitudes in the input MPAS mesh range from 0 to 360 is not relevant to this flag,
        as latitude and longitude are recomputed internally from Cartesian coordinates.
        Whether this flag is passed in or not, any longitudes written are in the 0-360 range.

Merging and Splitting

In order to support running MPAS-Albany Land Ice (MALI) with both Greenland and Antarctica at the same time, tools have been added to support merging and splitting MPAS meshes.

Merging two meshes can be accomplished with mpas_tools.merge_grids.merge_grids():

from mpas_tools.translate import translate
from mpas_tools.merge_grids import merge_grids
from mpas_tools.planar_hex import make_planar_hex_mesh
from mpas_tools.io import write_netcdf


dsMesh1 = make_planar_hex_mesh(nx=10, ny=10, dc=1000., nonperiodic_x=True,
                               nonperiodic_y=True)

dsMesh2 = make_planar_hex_mesh(nx=10, ny=10, dc=1000., nonperiodic_x=True,
                               nonperiodic_y=True)

translate(dsMesh2, xOffset=20000., yOffset=0.)

write_netcdf(dsMesh1, 'mesh1.nc')
write_netcdf(dsMesh2, 'mesh2.nc')

merge_grids(infile1='mesh1.nc', infile2='mesh2.nc',
            outfile='merged_mesh.nc')

Typically, it will only make sense to merge non-periodic meshes in this way.

Later, perhaps during analysis or visualization, it can be useful to split apart the merged meshes. This can be done with mpas_tools.split_grids.split_grids()

from mpas_tools.translate import translate
from mpas_tools.split_grids import split_grids
from mpas_tools.planar_hex import make_planar_hex_mesh
from mpas_tools.io import write_netcdf


dsMesh1 = make_planar_hex_mesh(nx=10, ny=10, dc=1000., nonperiodic_x=True,
                               nonperiodic_y=True)

dsMesh2 = make_planar_hex_mesh(nx=10, ny=10, dc=1000., nonperiodic_x=True,
                               nonperiodic_y=True)

translate(dsMesh2, xOffset=20000., yOffset=0.)

write_netcdf(dsMesh1, 'mesh1.nc')
write_netcdf(dsMesh2, 'mesh2.nc')


split_grids(infile='merged_mesh.nc', outfile1='split_mesh1.nc',
            outfile='split_mesh2.nc')

Merging meshes can also be accomplished with the merge_grids command-line tool:

$ merge_grids --help

usage: merge_grids [-h] [-o FILENAME] FILENAME1 FILENAME2

Tool to merge 2 MPAS non-contiguous meshes together into a single file

positional arguments:
  FILENAME1    File name for first mesh to merge
  FILENAME2    File name for second mesh to merge

optional arguments:
  -h, --help   show this help message and exit
  -o FILENAME  The merged mesh file

Similarly, split_grids can be used to to split meshes:

$ split_grids --help

usage: split_grids [-h] [-1 FILENAME] [-2 FILENAME] [--nCells NCELLS]
                   [--nEdges NEDGES] [--nVertices NVERTICES]
                   [--maxEdges MAXEDGES1 MAXEDGES2]
                   MESHFILE

Tool to split 2 previously merged MPAS non-contiguous meshes into separate files.
Typical usage is:
    split_grids.py -1 outfile1.nc -2 outfile2.nc infile
The optional arguments for nCells, nEdges, nVertices, and maxEdges should
generally not be required as this information is saved in the combined mesh file
as global attributes by the merge_grids.py script.

positional arguments:
  MESHFILE              Mesh file to split

optional arguments:
  -h, --help            show this help message and exit
  -1 FILENAME, --outfile1 FILENAME
                        File name for first mesh output
                        (default: mesh1.nc)
  -2 FILENAME, --outfile2 FILENAME
                        File name for second mesh output
                        (default: mesh2.nc)
  --nCells NCELLS       The number of cells in the first mesh
                        (default: the value specified in MESHFILE global attribute merge_point)
  --nEdges NEDGES       The number of edges in the first mesh
                        (default: the value specified in MESHFILE global attribute merge_point)
  --nVertices NVERTICES
                        The number of vertices in the first mesh
                        (default: the value specified in MESHFILE global attribute merge_point)
  --maxEdges MAXEDGES1 MAXEDGES2
                        The number of maxEdges in each mesh
                        (default: the value specified in MESHFILE global attribute merge_point
                              OR: will use MESHFILE maxEdges dimension and assume same for both)

Translation

A planar mesh can be translated in x, y or both by calling mpas_tools.translate.translate():

from mpas_tools.translate import translate
from mpas_tools.planar_hex import make_planar_hex_mesh

dsMesh = make_planar_hex_mesh(nx=10, ny=20, dc=1000., nonperiodic_x=False,
                              nonperiodic_y=False)

translate(dsMesh, xOffset=1000., yOffset=2000.)

This creates a periodic, planar mesh and then translates it by 1 km in x and 2 km in y.

Note

All the functions in the mpas_tools.translate module modify the mesh inplace, rather than returning a new xarray.Dataset object. This is in contrast to typical xarray functions and methods.

A mesh can be translated so that its center is at x = 0., y = 0. with the function mpas_tools.translate.center():

from mpas_tools.translate import center
from mpas_tools.planar_hex import make_planar_hex_mesh

dsMesh = make_planar_hex_mesh(nx=10, ny=20, dc=1000., nonperiodic_x=False,
                              nonperiodic_y=False)

center(dsMesh)

A mesh can be translated so its center matches the center of another mesh by using mpas_tools.translate.center_on_mesh():

from mpas_tools.translate import center_on_mesh
from mpas_tools.planar_hex import make_planar_hex_mesh

dsMesh1 = make_planar_hex_mesh(nx=10, ny=20, dc=1000., nonperiodic_x=False,
                               nonperiodic_y=False)

dsMesh2 = make_planar_hex_mesh(nx=20, ny=40, dc=2000., nonperiodic_x=False,
                               nonperiodic_y=False)

center_on_mesh(dsMesh2, dsMesh1)

In this example, the coordinates of dsMesh2 are altered so its center matches that of dsMesh1.

The functionality of all three of these functions is also available via the translate_planar_grid command-line tool:

$ translate_planar_grid --help

== Gathering information.  (Invoke with --help for more details. All arguments are optional)
Usage: translate_planar_grid [options]

This script translates the coordinate system of the planar MPAS mesh specified
with the -f flag.  There are 3 possible methods to choose from: 1) shift the
origin to the center of the domain 2) arbirary shift in x and/or y 3) shift to
the center of the domain described in a separate file

Options:
  -h, --help            show this help message and exit
  -f FILENAME, --file=FILENAME
                        MPAS planar grid file name. [default: grid.nc]
  -d FILENAME, --datafile=FILENAME
                        data file name to which to match the domain center of.
                        Uses xCell,yCell or, if those fields do not exist,
                        will secondly try x1,y1 fields.
  -x SHIFT_VALUE        user-specified shift in the x-direction. [default:
                        0.0]
  -y SHIFT_VALUE        user-specified shift in the y-direction. [default:
                        0.0]
  -c                    shift so origin is at center of domain [default:
                        False]

Converting Between Mesh Formats

MSH to MPAS NetCDF

jigsawpy produces meshes in .msh format that need to be converted to NetCDF files for use by MPAS components. A utility function mpas_tools.mesh.creation.jigsaw_to_netcdf.jigsaw_to_netcdf() or the command-line utility jigsaw_to_netcdf are used for this purpose.

In addition to the input .msh and output .nc files, the user must specify whether this is a spherical or planar mesh and, if it is spherical, provide the radius of the Earth in meters.

Triangle to MPAS NetCDF

Meshes in Triangle format can be converted to MPAS NetCDF format using mpas_tools.mesh.creation.triangle_to_netcdf.triangle_to_netcdf() or the triangle_to_netcdf command-line tool.

The user supplies the names of input .node and .ele files and the name of an output MPAS mesh file.

MPAS NetCDF to Triangle

MPAS meshes in NetCDF format can be converted to Triangle format using mpas_tools.mesh.creation.mpas_to_triangle.mpas_to_triangle() or the mpas_to_triangle command-line tool.

The user supplies the name of an input MPAS mesh file and the output prefix for the resulting Triangle .node and .ele files.

MPAS NetCDF to SCRIP

The function mpas_tools.scrip.from_mpas.scrip_from_mpas() can be used to convert an MPAS mesh file in NetCDF format to SCRIP format. SCRIP files are typically used to create mapping files used to interpolate between meshes.

A command-line tools is also available for this purpose:

$ scrip_from_mpas --help
== Gathering information.  (Invoke with --help for more details. All arguments are optional)
Usage: scrip_from_mpas [options]

This script takes an MPAS grid file and generates a SCRIP grid file.

Options:
  -h, --help            show this help message and exit
  -m FILENAME, --mpas=FILENAME
                        MPAS grid file name used as input. [default: grid.nc]
  -s FILENAME, --scrip=FILENAME
                        SCRIP grid file to output. [default: scrip.nc]
  -l, --landice         If flag is on, landice masks will be computed and
                        used.