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 (followingcellsOnCell) 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.