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.

Mask Creation with Python Multiprocessing

The Mask Creator is a serial code, and the algorithms it uses to find points in a region or cells, edges and vertices along a transect are not particularly efficient or sophisticated.

To provide better efficiency and to enable more sophisticated algorithms (now and in the future), a set of related python functions has been developed to provide much (but not all) of the functionality of the C++ Mask Creator described above.

Computing MPAS Region Masks

The function mpas_tools.mesh.mask.compute_mpas_region_masks() or the compute_mpas_region_masks command-line tool can be used to create region masks on cells, edges and/or vertices given an MPAS mesh xarray.Dataset dsMesh and a geometric_features.FeatureCollection fcMask containing regions. The resulting masks, in the variable regionCellMasks, are 1 where the center of the polygon corresponding to the cell, edge or vertex (see the MPAS Mesh Specification) are inside the given region and 0 where they are outside. This function is far more useful if the user provides a multiprocessing.Pool in the pool argument. pool should be created at the beginning of the calling code (when memory usage is small), possibly with mpas_tools.parallel.create_pool(), and terminated (multiprocessing.Pool.terminate()) before the code has finished. The same pool can be used in multiple calls to this and the other Python-based masking functions. If pool = None (the default), the masks are computed in serial, which will likely be frustratingly slow.

The chunkSize argument can be used to control how much work (how many cells, edges or vertices) each process computes on at one time. A very small chunkSize will incur a high overhead, while a very large chunkSize will lead to poor load balancing and infrequent progress updates (if showProgress = True). The default chunkSize of 1000 seems to perform well across a wide variety of mesh sizes and processor counts.

It is a good idea to provide a logger (see Logging) to get some output as the mask creation is progressing.

For efficiency, large shapes (e.g. the global coastline) are divided into smaller “tiles”. This subdivision is controlled with subdivisionThreshold, which should be set to a minimum size in degrees (latitude or longitude). If a shape is larger than this, it will be divided into tiles. The underlying algorithm first check bounding boxes of the resulting shapes against points before performing the more time-consuming step of determining if the point is inside the shape. The default value of 30 degrees performs much better than no subdivision for large shapes (again, such as the global coastline), but alternative values have not yet been explored.

The resulting variables are:

• regionCellMasks(nCells, nRegions) - a cell mask (1 for inside and 0 for outside the region) for each region

• regionEdgeMasks(nEdges, nRegions) - an edge mask for each region

• regionVertexMasks(nVertices, nRegions) - a vertex mask for each region

• regionNames(nRegions, string64) - the names of the regions

NetCDF fill values are used for invalid mask values, so nCellsInRegion, etc. are not produced.

The command-line tool takes the following arguments:

$ compute_mpas_region_masks --help
usage: compute_mpas_region_masks [-h] -m MESH_FILE_NAME -g GEOJSON_FILE_NAME
                                 -o MASK_FILE_NAME
                                 [-t MASK_TYPES [MASK_TYPES ...]]
                                 [-c CHUNK_SIZE] [--show_progress]
                                 [-s SUBDIVISION]
                                 [--process_count PROCESS_COUNT]
                                 [--multiprocessing_method MULTIPROCESSING_METHOD]

optional arguments:
  -h, --help            show this help message and exit
  -m MESH_FILE_NAME, --mesh_file_name MESH_FILE_NAME
                        An MPAS mesh file
  -g GEOJSON_FILE_NAME, --geojson_file_name GEOJSON_FILE_NAME
                        An Geojson file containing mask regions
  -o MASK_FILE_NAME, --mask_file_name MASK_FILE_NAME
                        An output MPAS region masks file
  -t MASK_TYPES [MASK_TYPES ...], --mask_types MASK_TYPES [MASK_TYPES ...]
                        Which type(s) of masks to make: cell, edge or vertex.
                        Default is cell and vertex.
  -c CHUNK_SIZE, --chunk_size CHUNK_SIZE
                        The number of cells, vertices or edges that are
                        processed in one operation
  --show_progress       Whether to show a progress bar
  -s SUBDIVISION, --subdivision SUBDIVISION
                        A threshold in degrees (lon or lat) above which the
                        mask region will be subdivided into smaller polygons
                        for faster intersection checking
  --process_count PROCESS_COUNT
                        The number of processes to use to compute masks. The
                        default is to use all available cores
  --multiprocessing_method MULTIPROCESSING_METHOD
                        The multiprocessing method use for python mask
                        creation ('fork', 'spawn' or 'forkserver')

Computing Transect Masks

The function mpas_tools.mesh.mask.compute_mpas_transect_masks() and the compute_mpas_transect_masks command-line tool are similar to the function for computing region masks. The function takes a geometric_features.FeatureCollection fcMask that is made up of transects, rather than regions. One mask is produced for each feature in the collection, indicating where the transect intersects the cell, edge or vertex polygons (see the MPAS Mesh Specification).

The arguments logger, pool, chunkSize and showProgress are the same as for region-mask creation above.

The argument subdivisionResolution is a length in meters, above which segments of the transect are subdivided to provide a better representation of the spherical path in longitude/latitude space. The default value of 10 km is typically good enough to capture distortion at typical MPAS mesh resolutions.

The algorithm perform intersections in longitude/latitude space using the shapely library. Because shapely is designed for 2D shapes in a Cartesian plane, it is not designed for spherical coordinates. Care has been taken to handle periodicity at the dateline (antimeridian) but there may be issues with MPAS mesh polygons containing the north or south pole. If a user needs to handle a transect that is very close to the pole, it is likely worth contacting the developers to request modifications to the code to support this case.

The resulting variables are:

• transectCellMasks(nCells, nTransects) - a cell mask (1 if the transect intersects the cell and 0 if not) for each transect

• transectEdgeMasks(nEdges, nTransects) - an edge mask for each transect

• transectVertexMasks(nVertices, nTransects) - a vertex mask for each transect

• transectNames(nTransects, string64) - the names of the transects

We don’t currently provide cell, edge or vertex indices (e.g. transectCellGlobalIDs) for path along a transect. This is, in part, because the algorithm doesn’t keep track of the relative order of points along a transect. This could be updated in the future if there is sufficient demand.

The edge sign (transectEdgeMaskSigns) is computed only if addEdgeSign=True, since this takes extra time to compute and isn’t always needed.

Note

While the default subdivisionResolution is 10 km for mpas_tools.mesh.mask.compute_mpas_transect_masks(), the default behavior in the command-line tool compute_mpas_transect_masks is no subdivision because there is otherwise not a good way to specify at the command line that no subdivision is desired. Typically, users will want to request subdivision with something like -s 10e3

The command-line tool takes the following arguments:

$ compute_mpas_transect_masks --help
usage: compute_mpas_transect_masks [-h] -m MESH_FILE_NAME -g GEOJSON_FILE_NAME
                                   -o MASK_FILE_NAME
                                   [-t MASK_TYPES [MASK_TYPES ...]]
                                   [-c CHUNK_SIZE] [--show_progress]
                                   [-s SUBDIVISION]
                                   [--process_count PROCESS_COUNT]
                                   [--multiprocessing_method MULTIPROCESSING_METHOD]

optional arguments:
  -h, --help            show this help message and exit
  -m MESH_FILE_NAME, --mesh_file_name MESH_FILE_NAME
                        An MPAS mesh file
  -g GEOJSON_FILE_NAME, --geojson_file_name GEOJSON_FILE_NAME
                        An Geojson file containing transects
  -o MASK_FILE_NAME, --mask_file_name MASK_FILE_NAME
                        An output MPAS transect masks file
  -t MASK_TYPES [MASK_TYPES ...], --mask_types MASK_TYPES [MASK_TYPES ...]
                        Which type(s) of masks to make: cell, edge or vertex.
                        Default is cell, edge and vertex.
  -c CHUNK_SIZE, --chunk_size CHUNK_SIZE
                        The number of cells, vertices or edges that are
                        processed in one operation
  --show_progress       Whether to show a progress bar
  -s SUBDIVISION, --subdivision SUBDIVISION
                        The maximum resolution (in meters) of segments in a
                        transect. If a transect is too coarse, it will be
                        subdivided. Default is no subdivision.
  --process_count PROCESS_COUNT
                        The number of processes to use to compute masks. The
                        default is to use all available cores
  --multiprocessing_method MULTIPROCESSING_METHOD
                        The multiprocessing method use for python mask
                        creation ('fork', 'spawn' or 'forkserver')
  --add_edge_sign       Whether to add the transectEdgeMaskSigns variable

Computing a Flood-fill Mask

The function mpas_tools.mesh.mask.compute_mpas_flood_fill_mask() and the command-line tool compute_mpas_flood_fill_mask fill in a mask, starting with the cell centers closest to the seed points given in geometric_features.FeatureCollection fcSeed. This algorithm runs in serial, and will be more efficient the more seed points are provided and the more widely scattered over the mesh they are.

An optional daGrow argument to the function (not currently available from the command-line tool) provides a mask into which the flood fill is allowed to grow. The default is all ones.

The resulting dataset contains a single variable:

• cellSeedMask(nCells) - a cell mask that is 1 where the flood fill (following cellsOnCell) propagated starting from the seed points and 0 elsewhere

The command-line tool takes the following arguments:

$ compute_mpas_flood_fill_mask --help
usage: compute_mpas_flood_fill_mask [-h] -m MESH_FILE_NAME -g
                                    GEOJSON_FILE_NAME -o MASK_FILE_NAME

optional arguments:
  -h, --help            show this help message and exit
  -m MESH_FILE_NAME, --mesh_file_name MESH_FILE_NAME
                        An MPAS mesh file
  -g GEOJSON_FILE_NAME, --geojson_file_name GEOJSON_FILE_NAME
                        An Geojson file containing points at which to start
                        the flood fill
  -o MASK_FILE_NAME, --mask_file_name MASK_FILE_NAME
                        An output MPAS region masks file

Computing Lon/Lat Region Masks

The function mpas_tools.mesh.mask.compute_lon_lat_region_masks() or the compute_lon_lat_region_masks command-line tool compute region masks on a longitude/latitude grid but are otherwise functionally very similar to the corresponding tools for compute MPAS region masks. The major difference is that 1D arrays of longitude and latitude are provided instead of an MPAS mesh dataset. There is no argument equivalent to the mask type for MPAS meshes. Instead, mask values are given at each point on the 2D longitude/latitude grid. All other arguments serve the same purpose as for the MPAS region mask creation described above.

The command-line tool takes the following arguments:

$ compute_lon_lat_region_masks --help
usage: compute_lon_lat_region_masks [-h] -i GRID_FILE_NAME [--lon LON]
                                    [--lat LAT] -g GEOJSON_FILE_NAME -o
                                    MASK_FILE_NAME [-c CHUNK_SIZE]
                                    [--show_progress] [-s SUBDIVISION]
                                    [--process_count PROCESS_COUNT]
                                    [--multiprocessing_method MULTIPROCESSING_METHOD]

optional arguments:
  -h, --help            show this help message and exit
  -i GRID_FILE_NAME, --grid_file_name GRID_FILE_NAME
                        An input lon/lat grid file
  --lon LON             The name of the longitude coordinate
  --lat LAT             The name of the latitude coordinate
  -g GEOJSON_FILE_NAME, --geojson_file_name GEOJSON_FILE_NAME
                        An Geojson file containing mask regions
  -o MASK_FILE_NAME, --mask_file_name MASK_FILE_NAME
                        An output MPAS region masks file
  -c CHUNK_SIZE, --chunk_size CHUNK_SIZE
                        The number of grid points that are processed in one
                        operation
  --show_progress       Whether to show a progress bar
  -s SUBDIVISION, --subdivision SUBDIVISION
                        A threshold in degrees (lon or lat) above which the
                        mask region will be subdivided into smaller polygons
                        for faster intersection checking
  --process_count PROCESS_COUNT
                        The number of processes to use to compute masks. The
                        default is to use all available cores
  --multiprocessing_method MULTIPROCESSING_METHOD
                        The multiprocessing method use for python mask
                        creation ('fork', 'spawn' or 'forkserver')

Culling MPAS Datasets

The tools described in Cell Culler can be used to create a culled horizontal MPAS mesh. Once a culled MPAS mesh has been created, an MPAS dataset on the unculled mesh can be cropped to the culled mesh using the the mpas_tools.mesh.cull.cull_dataset() or mpas_tools.mesh.cull.write_culled_dataset() functions. These functions take a dataset (or filename) to crop as well as datasets (or filenames) for the unculled and culled horizontal MPAS meshes. They return (or write out) the culled version of the data set. Fields that exist in the culled horizonal mesh are copied from the culled mesh, rather than cropped from the dataset. This because we wish to keep the cropped horizontal mesh exactly as it was produced by the culling tool, which may not correspond to a cropped version of the field from the original mesh. For example, fields are reindexed during culling and coordinates are recomputed.

It may be useful to compute and store the maps from cells, edges and vertices on the culled mesh back to the unculled mesh for reuse. This can be accomplished by calling the mpas_tools.mesh.cull.map_culled_to_base() or mpas_tools.mesh.cull.write_map_culled_to_base() functions.

An example workflow that culls out ice-shelf cavities from an MPAS-Ocean initial condition might look like the following. In this case the file culled_mesh.nc is a mesh where land (and the grounded portion of the ice sheet) has been removed but where ice-shelf cavities are still present. It serves as the “base” mesh for the purposes of this example. culled_mesh_no_isc.nc is created (if it doesn’t already exist) with the ice-shelf cavities removed as well, so it is the “culled” mesh in this example. We store the mapping betwen the two horizontal meshes in no_isc_to_culled_map.nc in case we want to resue it later. The initial condition is read from initial_state.nc and the culled version is written to initial_state_no_isc.nc:

import os

import xarray as xr

from mpas_tools.io import write_netcdf
from mpas_tools.mesh.conversion import cull
from mpas_tools.mesh.cull import write_map_culled_to_base, write_culled_dataset
from mpas_tools.logging import LoggingContext


in_filename = 'initial_state.nc'
out_filename = 'initial_state_no_isc.nc'
base_mesh_filename = 'culled_mesh.nc'
culled_mesh_filename = 'culled_mesh_no_isc.nc'
map_filename = 'no_isc_to_culled_map.nc'

if not os.path.exists(culled_mesh_filename):
    ds_culled_mesh = xr.open_dataset(base_mesh_filename)
    ds_init = xr.open_dataset(in_filename)
    ds_culled_mesh['cullCell'] = ds_init.landIceMask
    ds_culled_mesh_no_isc = cull(ds_culled_mesh)
    write_netcdf(ds_culled_mesh_no_isc, culled_mesh_filename)

if not os.path.exists(map_filename):
    write_map_culled_to_base(base_mesh_filename=base_mesh_filename,
                             culled_mesh_filename=culled_mesh_filename,
                             out_filename=map_filename)

with LoggingContext('test') as logger:
    write_culled_dataset(in_filename=in_filename, out_filename=out_filename,
                         base_mesh_filename=base_mesh_filename,
                         culled_mesh_filename=culled_mesh_filename,
                         map_culled_to_base_filename=map_filename,
                         logger=logger)

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.