From 8467d43c637e14b1498362ce65b796eea85925ba Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:53:47 -0500 Subject: [PATCH 01/10] Apply pre-commit formatting to planar_hex and triangle_to_netcdf These modules predate the ruff and flynt hooks. Reformat them before changing them so the functional changes are easy to review. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../mesh/creation/triangle_to_netcdf.py | 91 +++--- conda_package/mpas_tools/planar_hex.py | 269 ++++++++++++------ 2 files changed, 228 insertions(+), 132 deletions(-) diff --git a/conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py b/conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py index 48238f8b8..86ba339e3 100644 --- a/conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py +++ b/conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py @@ -1,12 +1,16 @@ -from __future__ import absolute_import, division, print_function, \ - unicode_literals +from __future__ import ( + absolute_import, + division, + print_function, + unicode_literals, +) -import numpy as np +import argparse +import numpy as np from netCDF4 import Dataset as NetCDFFile -from mpas_tools.mesh.creation.util import circumcenter -import argparse +from mpas_tools.mesh.creation.util import circumcenter def triangle_to_netcdf(node, ele, output_name): @@ -30,8 +34,8 @@ def triangle_to_netcdf(node, ele, output_name): # Get nCells cell_info = open(node, 'r') nCells = -1 # There is one header line - for block in iter(lambda: cell_info.readline(), ""): - if block.startswith("#"): + for block in iter(lambda: cell_info.readline(), ''): + if block.startswith('#'): continue # skip comment lines nCells = nCells + 1 cell_info.close() @@ -40,15 +44,17 @@ def triangle_to_netcdf(node, ele, output_name): cov_info = open(ele, 'r') vertexDegree = 3 # always triangles with Triangle! nVertices = -1 # There is one header line - for block in iter(lambda: cov_info.readline(), ""): - if block.startswith("#"): + for block in iter(lambda: cov_info.readline(), ''): + if block.startswith('#'): continue # skip comment lines nVertices = nVertices + 1 cov_info.close() if vertexDegree != 3: - ValueError("This script can only compute vertices with triangular " - "dual meshes currently.") + ValueError( + 'This script can only compute vertices with triangular ' + 'dual meshes currently.' + ) grid.createDimension('nCells', nCells) grid.createDimension('nVertices', nVertices) @@ -62,9 +68,9 @@ def triangle_to_netcdf(node, ele, output_name): cell_info = open(node, 'r') cell_info.readline() # read header i = 0 - for block in iter(lambda: cell_info.readline(), ""): + for block in iter(lambda: cell_info.readline(), ''): block_arr = block.split() - if block_arr[0] == "#": + if block_arr[0] == '#': continue # skip comment lines xCell_full[i] = float(block_arr[1]) yCell_full[i] = float(block_arr[2]) @@ -72,18 +78,17 @@ def triangle_to_netcdf(node, ele, output_name): i = i + 1 cell_info.close() - grid.on_a_sphere = "NO" + grid.on_a_sphere = 'NO' grid.sphere_radius = 0.0 - cellsOnVertex_full = np.zeros( - (nVertices, vertexDegree), dtype=np.int32) + cellsOnVertex_full = np.zeros((nVertices, vertexDegree), dtype=np.int32) cov_info = open(ele, 'r') cov_info.readline() # read header iVertex = 0 - for block in iter(lambda: cov_info.readline(), ""): + for block in iter(lambda: cov_info.readline(), ''): block_arr = block.split() - if block_arr[0] == "#": + if block_arr[0] == '#': continue # skip comment lines cellsOnVertex_full[iVertex, :] = int(-1) # skip the first column, which is the triangle number, and then @@ -120,8 +125,7 @@ def triangle_to_netcdf(node, ele, output_name): yVertex_full[iVertex] = pv.y zVertex_full[iVertex] = pv.z - meshDensity_full = grid.createVariable( - 'meshDensity', 'f8', ('nCells',)) + meshDensity_full = grid.createVariable('meshDensity', 'f8', ('nCells',)) meshDensity_full[0:nCells] = 1.0 @@ -138,7 +142,13 @@ def triangle_to_netcdf(node, ele, output_name): var = grid.createVariable('zVertex', 'f8', ('nVertices',)) var[:] = zVertex_full var = grid.createVariable( - 'cellsOnVertex', 'i4', ('nVertices', 'vertexDegree',)) + 'cellsOnVertex', + 'i4', + ( + 'nVertices', + 'vertexDegree', + ), + ) var[:] = cellsOnVertex_full grid.sync() @@ -147,29 +157,32 @@ def triangle_to_netcdf(node, ele, output_name): def main(): parser = argparse.ArgumentParser( - description=__doc__, - formatter_class=argparse.RawTextHelpFormatter) + description=__doc__, formatter_class=argparse.RawTextHelpFormatter + ) parser.add_argument( - "-n", - "--node", - dest="node", + '-n', + '--node', + dest='node', required=True, - help="input .node file generated by Triangle.", - metavar="FILE") + help='input .node file generated by Triangle.', + metavar='FILE', + ) parser.add_argument( - "-e", - "--ele", - dest="ele", + '-e', + '--ele', + dest='ele', required=True, - help="input .ele file generated by Triangle.", - metavar="FILE") + help='input .ele file generated by Triangle.', + metavar='FILE', + ) parser.add_argument( - "-o", - "--output", - dest="output", - default="grid.nc", - help="output file name.", - metavar="FILE") + '-o', + '--output', + dest='output', + default='grid.nc', + help='output file name.', + metavar='FILE', + ) options = parser.parse_args() triangle_to_netcdf(options.node, options.ele, options.output) diff --git a/conda_package/mpas_tools/planar_hex.py b/conda_package/mpas_tools/planar_hex.py index 3a32ebd0b..335377f81 100755 --- a/conda_package/mpas_tools/planar_hex.py +++ b/conda_package/mpas_tools/planar_hex.py @@ -1,19 +1,31 @@ #!/usr/bin/env python -from __future__ import absolute_import, division, print_function, \ - unicode_literals +from __future__ import ( + absolute_import, + division, + print_function, + unicode_literals, +) + +import argparse import numpy import xarray -import argparse from mpas_tools.io import write_netcdf -def make_planar_hex_mesh(nx, ny, dc, nonperiodic_x, - nonperiodic_y, outFileName=None, - compareWithFileName=None, - format=None, engine=None): +def make_planar_hex_mesh( + nx, + ny, + dc, + nonperiodic_x, + nonperiodic_y, + outFileName=None, + compareWithFileName=None, + format=None, + engine=None, +): """ Builds an MPAS periodic, planar hexagonal mesh with the requested dimensions, optionally saving it to a file, and returns it as an @@ -85,8 +97,10 @@ def make_planar_hex_mesh(nx, ny, dc, nonperiodic_x, def initial_setup(nx, ny, dc, nonperiodic_x, nonperiodic_y): """Setup the dimensions and add placeholders for some index variables""" if ny % 2 != 0: - raise ValueError('ny must be divisible by 2 for the grid\'s ' - 'periodicity to work properly.') + raise ValueError( + "ny must be divisible by 2 for the grid's " + 'periodicity to work properly.' + ) mesh = xarray.Dataset() @@ -96,19 +110,19 @@ def initial_setup(nx, ny, dc, nonperiodic_x, nonperiodic_y): mesh.attrs['is_periodic'] = 'YES' if nonperiodic_x: - mesh.attrs['x_period'] = 0. + mesh.attrs['x_period'] = 0.0 else: mesh.attrs['x_period'] = nx * dc if nonperiodic_y: - mesh.attrs['y_period'] = 0. + mesh.attrs['y_period'] = 0.0 else: - mesh.attrs['y_period'] = ny * dc * numpy.sqrt(3.) / 2. + mesh.attrs['y_period'] = ny * dc * numpy.sqrt(3.0) / 2.0 mesh.attrs['dc'] = dc mesh.attrs['nx'] = nx mesh.attrs['ny'] = ny mesh.attrs['on_a_sphere'] = 'NO' - mesh.attrs['sphere_radius'] = 0. + mesh.attrs['sphere_radius'] = 0.0 if nonperiodic_x: nx = nx + 2 @@ -127,8 +141,9 @@ def initial_setup(nx, ny, dc, nonperiodic_x, nonperiodic_y): indexToVertexID = numpy.arange(nVertices, dtype='i4') cellIdx = indexToCellID.reshape(ny, nx) - cellCol, cellRow = numpy.meshgrid(numpy.arange(nx, dtype='i4'), - numpy.arange(ny, dtype='i4')) + cellCol, cellRow = numpy.meshgrid( + numpy.arange(nx, dtype='i4'), numpy.arange(ny, dtype='i4') + ) mesh['cellIdx'] = (('ny', 'nx'), cellIdx) mesh['cellRow'] = (('nCells',), cellRow.ravel()) @@ -141,25 +156,38 @@ def initial_setup(nx, ny, dc, nonperiodic_x, nonperiodic_y): mesh['cullCell'] = (('nCells',), numpy.zeros(nCells, 'i4')) mesh['nEdgesOnCell'] = (('nCells',), 6 * numpy.ones((nCells,), 'i4')) - mesh['cellsOnCell'] = (('nCells', 'maxEdges'), - numpy.zeros((nCells, maxEdges), 'i4')) - mesh['edgesOnCell'] = (('nCells', 'maxEdges'), - numpy.zeros((nCells, maxEdges), 'i4')) - mesh['verticesOnCell'] = (('nCells', 'maxEdges'), - numpy.zeros((nCells, maxEdges), 'i4')) + mesh['cellsOnCell'] = ( + ('nCells', 'maxEdges'), + numpy.zeros((nCells, maxEdges), 'i4'), + ) + mesh['edgesOnCell'] = ( + ('nCells', 'maxEdges'), + numpy.zeros((nCells, maxEdges), 'i4'), + ) + mesh['verticesOnCell'] = ( + ('nCells', 'maxEdges'), + numpy.zeros((nCells, maxEdges), 'i4'), + ) mesh['nEdgesOnEdge'] = (('nEdges',), 10 * numpy.ones((nEdges,), 'i4')) - mesh['cellsOnEdge'] = (('nEdges', 'TWO'), - numpy.zeros((nEdges, 2), 'i4')) - mesh['edgesOnEdge'] = (('nEdges', 'maxEdges2'), - -1 * numpy.ones((nEdges, 2 * maxEdges), 'i4')) - mesh['verticesOnEdge'] = (('nEdges', 'TWO'), - numpy.zeros((nEdges, 2), 'i4')) - - mesh['cellsOnVertex'] = (('nVertices', 'vertexDegree'), - numpy.zeros((nVertices, vertexDegree), 'i4')) - mesh['edgesOnVertex'] = (('nVertices', 'vertexDegree'), - numpy.zeros((nVertices, vertexDegree), 'i4')) + mesh['cellsOnEdge'] = (('nEdges', 'TWO'), numpy.zeros((nEdges, 2), 'i4')) + mesh['edgesOnEdge'] = ( + ('nEdges', 'maxEdges2'), + -1 * numpy.ones((nEdges, 2 * maxEdges), 'i4'), + ) + mesh['verticesOnEdge'] = ( + ('nEdges', 'TWO'), + numpy.zeros((nEdges, 2), 'i4'), + ) + + mesh['cellsOnVertex'] = ( + ('nVertices', 'vertexDegree'), + numpy.zeros((nVertices, vertexDegree), 'i4'), + ) + mesh['edgesOnVertex'] = ( + ('nVertices', 'vertexDegree'), + numpy.zeros((nVertices, vertexDegree), 'i4'), + ) return mesh @@ -170,7 +198,7 @@ def mark_cull_cell_nonperiodic_y(mesh): nCells = mesh.sizes['nCells'] nx = mesh.sizes['nx'] cullCell[0:nx] = 1 - cullCell[nCells - nx:nCells + 1] = 1 + cullCell[nCells - nx : nCells + 1] = 1 def mark_cull_cell_nonperiodic_x(mesh): @@ -179,7 +207,7 @@ def mark_cull_cell_nonperiodic_x(mesh): nCells = mesh.sizes['nCells'] nx = mesh.sizes['nx'] cullCell[::nx] = 1 - cullCell[nx - 1:nCells + 1:nx] = 1 + cullCell[nx - 1 : nCells + 1 : nx] = 1 def compute_indices_on_cell(mesh): @@ -309,14 +337,19 @@ def compute_weights_on_edge(mesh): nEdges = mesh.sizes['nEdges'] maxEdges2 = mesh.sizes['maxEdges2'] - mesh['weightsOnEdge'] = (('nEdges', 'maxEdges2'), - numpy.zeros((nEdges, maxEdges2), 'f8')) + mesh['weightsOnEdge'] = ( + ('nEdges', 'maxEdges2'), + numpy.zeros((nEdges, maxEdges2), 'f8'), + ) weightsOnEdge = mesh.weightsOnEdge - weights = (1. / numpy.sqrt(3.)) * numpy.array( - [[1. / 3., 1. / 6., 0., 1. / 6., 1. / 3.], - [1. / 3., -1. / 6., 0., 1. / 6., -1. / 3.], - [-1. / 3., -1. / 6., 0., -1. / 6., -1. / 3.]]) + weights = (1.0 / numpy.sqrt(3.0)) * numpy.array( + [ + [1.0 / 3.0, 1.0 / 6.0, 0.0, 1.0 / 6.0, 1.0 / 3.0], + [1.0 / 3.0, -1.0 / 6.0, 0.0, 1.0 / 6.0, -1.0 / 3.0], + [-1.0 / 3.0, -1.0 / 6.0, 0.0, -1.0 / 6.0, -1.0 / 3.0], + ] + ) for i in range(3): for j in range(5): weightsOnEdge[edgesOnCell[:, i + 3], j] = weights[i, j] @@ -350,7 +383,7 @@ def compute_coordinates(mesh): mask = numpy.mod(cellRow, 2) == 0 mesh['xCell'] = (dc * (cellCol + 0.5)).where(mask, dc * (cellCol + 1)) - mesh['yCell'] = dc * (cellRow + 1) * numpy.sqrt(3.) / 2. + mesh['yCell'] = dc * (cellRow + 1) * numpy.sqrt(3.0) / 2.0 mesh['zCell'] = (('nCells',), numpy.zeros((nCells,), 'f8')) mesh['xEdge'] = (('nEdges',), numpy.zeros((nEdges,), 'f8')) @@ -360,54 +393,77 @@ def compute_coordinates(mesh): mesh.xEdge[edgesOnCell[:, 0]] = mesh.xCell - 0.5 * dc mesh.yEdge[edgesOnCell[:, 0]] = mesh.yCell - mesh.xEdge[edgesOnCell[:, 1]] = mesh.xCell - \ - 0.5 * dc * numpy.cos(numpy.pi / 3.) - mesh.yEdge[edgesOnCell[:, 1]] = mesh.yCell - \ - 0.5 * dc * numpy.sin(numpy.pi / 3.) + mesh.xEdge[edgesOnCell[:, 1]] = mesh.xCell - 0.5 * dc * numpy.cos( + numpy.pi / 3.0 + ) + mesh.yEdge[edgesOnCell[:, 1]] = mesh.yCell - 0.5 * dc * numpy.sin( + numpy.pi / 3.0 + ) - mesh.xEdge[edgesOnCell[:, 2]] = mesh.xCell + \ - 0.5 * dc * numpy.cos(numpy.pi / 3.) - mesh.yEdge[edgesOnCell[:, 2]] = mesh.yCell - \ - 0.5 * dc * numpy.sin(numpy.pi / 3.) + mesh.xEdge[edgesOnCell[:, 2]] = mesh.xCell + 0.5 * dc * numpy.cos( + numpy.pi / 3.0 + ) + mesh.yEdge[edgesOnCell[:, 2]] = mesh.yCell - 0.5 * dc * numpy.sin( + numpy.pi / 3.0 + ) mesh['xVertex'] = (('nVertices',), numpy.zeros((nVertices,), 'f8')) mesh['yVertex'] = (('nVertices',), numpy.zeros((nVertices,), 'f8')) mesh['zVertex'] = (('nVertices',), numpy.zeros((nVertices,), 'f8')) mesh.xVertex[verticesOnCell[:, 0]] = mesh.xCell - 0.5 * dc - mesh.yVertex[verticesOnCell[:, 0]] = mesh.yCell + dc * numpy.sqrt(3.) / 6. + mesh.yVertex[verticesOnCell[:, 0]] = ( + mesh.yCell + dc * numpy.sqrt(3.0) / 6.0 + ) mesh.xVertex[verticesOnCell[:, 1]] = mesh.xCell - 0.5 * dc - mesh.yVertex[verticesOnCell[:, 1]] = mesh.yCell - dc * numpy.sqrt(3.) / 6. + mesh.yVertex[verticesOnCell[:, 1]] = ( + mesh.yCell - dc * numpy.sqrt(3.0) / 6.0 + ) mesh['angleEdge'] = (('nEdges',), numpy.zeros((nEdges,), 'f8')) - mesh.angleEdge[edgesOnCell[:, 1]] = numpy.pi / 3. - mesh.angleEdge[edgesOnCell[:, 2]] = 2. * numpy.pi / 3. + mesh.angleEdge[edgesOnCell[:, 1]] = numpy.pi / 3.0 + mesh.angleEdge[edgesOnCell[:, 2]] = 2.0 * numpy.pi / 3.0 mesh['dcEdge'] = (('nEdges',), dc * numpy.ones((nEdges,), 'f8')) - mesh['dvEdge'] = mesh.dcEdge * numpy.sqrt(3.) / 3. - - mesh['areaCell'] = \ - (('nCells',), dc**2 * numpy.sqrt(3.) / 2. * numpy.ones((nCells,), 'f8')) - - mesh['areaTriangle'] = \ - (('nVertices',), dc**2 * numpy.sqrt(3.) / - 4. * numpy.ones((nVertices,), 'f8')) - - mesh['kiteAreasOnVertex'] = \ - (('nVertices', 'vertexDegree'), - dc**2 * numpy.sqrt(3.) / 12. * numpy.ones((nVertices, vertexDegree), - 'f8')) + mesh['dvEdge'] = mesh.dcEdge * numpy.sqrt(3.0) / 3.0 + + mesh['areaCell'] = ( + ('nCells',), + dc**2 * numpy.sqrt(3.0) / 2.0 * numpy.ones((nCells,), 'f8'), + ) + + mesh['areaTriangle'] = ( + ('nVertices',), + dc**2 * numpy.sqrt(3.0) / 4.0 * numpy.ones((nVertices,), 'f8'), + ) + + mesh['kiteAreasOnVertex'] = ( + ('nVertices', 'vertexDegree'), + dc**2 + * numpy.sqrt(3.0) + / 12.0 + * numpy.ones((nVertices, vertexDegree), 'f8'), + ) mesh['meshDensity'] = (('nCells',), numpy.ones((nCells,), 'f8')) def add_one_to_indices(mesh): """Needed to adhere to Fortran indexing""" - indexVars = ['indexToCellID', 'indexToEdgeID', 'indexToVertexID', - 'cellsOnCell', 'edgesOnCell', 'verticesOnCell', - 'cellsOnEdge', 'edgesOnEdge', 'verticesOnEdge', - 'cellsOnVertex', 'edgesOnVertex'] + indexVars = [ + 'indexToCellID', + 'indexToEdgeID', + 'indexToVertexID', + 'cellsOnCell', + 'edgesOnCell', + 'verticesOnCell', + 'cellsOnEdge', + 'edgesOnEdge', + 'verticesOnEdge', + 'cellsOnVertex', + 'edgesOnVertex', + ] for var in indexVars: mesh[var] = mesh[var] + 1 @@ -421,19 +477,19 @@ def make_diff(mesh, refMeshFileName, diffFileName): diff[variable] = mesh[variable] - refMesh[variable] print(diff[variable].name, float(numpy.abs(diff[variable]).max())) else: - print('mesh has extra variable {}'.format(mesh[variable].name)) + print(f'mesh has extra variable {mesh[variable].name}') for variable in refMesh.data_vars: if variable not in mesh: - print('mesh mising variable {}'.format(refMesh[variable].name)) + print(f'mesh mising variable {refMesh[variable].name}') for attr in refMesh.attrs: if attr not in mesh.attrs: - print('mesh mising attribute {}'.format(attr)) + print(f'mesh mising attribute {attr}') for attr in mesh.attrs: if attr not in refMesh.attrs: - print('mesh has extra attribute {}'.format(attr)) + print(f'mesh has extra attribute {attr}') write_netcdf(diff, diffFileName) @@ -441,28 +497,55 @@ def make_diff(mesh, refMeshFileName, diffFileName): def main(): parser = argparse.ArgumentParser( - description=__doc__, formatter_class=argparse.RawTextHelpFormatter) - parser.add_argument('--nx', dest='nx', type=int, required=True, - help='Cells in x direction') - parser.add_argument('--ny', dest='ny', type=int, required=True, - help='Cells in y direction') - parser.add_argument('--dc', dest='dc', type=float, required=True, - help='Distance between cell centers in meters') - parser.add_argument('--npx', '--nonperiodic_x', dest='nonperiodic_x', - action="store_true", - help='non-periodic in x direction') - parser.add_argument('--npy', '--nonperiodic_y', dest='nonperiodic_y', - action="store_true", - help='non-periodic in y direction') - parser.add_argument('-o', '--outFileName', dest='outFileName', type=str, - required=False, default='grid.nc', - help='The name of the output file') + description=__doc__, formatter_class=argparse.RawTextHelpFormatter + ) + parser.add_argument( + '--nx', dest='nx', type=int, required=True, help='Cells in x direction' + ) + parser.add_argument( + '--ny', dest='ny', type=int, required=True, help='Cells in y direction' + ) + parser.add_argument( + '--dc', + dest='dc', + type=float, + required=True, + help='Distance between cell centers in meters', + ) + parser.add_argument( + '--npx', + '--nonperiodic_x', + dest='nonperiodic_x', + action='store_true', + help='non-periodic in x direction', + ) + parser.add_argument( + '--npy', + '--nonperiodic_y', + dest='nonperiodic_y', + action='store_true', + help='non-periodic in y direction', + ) + parser.add_argument( + '-o', + '--outFileName', + dest='outFileName', + type=str, + required=False, + default='grid.nc', + help='The name of the output file', + ) args = parser.parse_args() - make_planar_hex_mesh(args.nx, args.ny, args.dc, - args.nonperiodic_x, args.nonperiodic_y, - args.outFileName) + make_planar_hex_mesh( + args.nx, + args.ny, + args.dc, + args.nonperiodic_x, + args.nonperiodic_y, + args.outFileName, + ) if __name__ == '__main__': From 7a2c31dc18704c9f8bbb8179e15497fcb491117f Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:54:38 -0500 Subject: [PATCH 02/10] Add CF metadata for MPAS mesh variables Add mpas_tools.mesh.attrs with a table of long_name, units and standard_name for the MPAS mesh variables, cf_conventions() to add CF-1.8 (unless another CF version is present) and MPAS to a Conventions attribute, and add_mesh_attrs() to apply both to a mesh dataset. The table is adapted from the one in E3SM-Project/polaris#784, with standard names for latitude, longitude and cell area. Co-Authored-By: Claude Opus 5.5 (1M context) --- conda_package/mpas_tools/mesh/attrs.py | 216 +++++++++++++++++++++++++ 1 file changed, 216 insertions(+) create mode 100644 conda_package/mpas_tools/mesh/attrs.py diff --git a/conda_package/mpas_tools/mesh/attrs.py b/conda_package/mpas_tools/mesh/attrs.py new file mode 100644 index 000000000..f236384c3 --- /dev/null +++ b/conda_package/mpas_tools/mesh/attrs.py @@ -0,0 +1,216 @@ +""" +CF metadata for MPAS mesh variables and files + +``MESH_VAR_ATTRS`` holds the ``long_name``, ``units`` and ``standard_name`` +attributes of the MPAS mesh variables. Units are in the plain form that +udunits parses and CF uses (``m``, ``m2``, ``radians``, ``1`` for +dimensionless). Index and connectivity variables and masks have no units. + +``MpasMeshConverter.x`` and ``MpasCellCuller.x`` write the same attributes, +so the two need to be kept in sync. +""" + +CF_VERSION = 'CF-1.8' + +MESH_VAR_ATTRS = { + # global indices + 'indexToCellID': {'long_name': 'global index of each cell'}, + 'indexToEdgeID': {'long_name': 'global index of each edge'}, + 'indexToVertexID': {'long_name': 'global index of each vertex'}, + # connectivity + 'cellsOnCell': {'long_name': 'cells that neighbor each cell'}, + 'nEdgesOnCell': {'long_name': 'number of edges that border each cell'}, + 'edgesOnCell': {'long_name': 'edges that border each cell'}, + 'verticesOnCell': {'long_name': 'vertices that border each cell'}, + 'cellsOnEdge': {'long_name': 'cells that straddle each edge'}, + 'edgesOnEdge': { + 'long_name': 'edges that border the cells that straddle each edge' + }, + 'nEdgesOnEdge': { + 'long_name': 'number of edges that border the cells that straddle ' + 'each edge' + }, + 'verticesOnEdge': {'long_name': 'vertices that straddle each edge'}, + 'cellsOnVertex': {'long_name': 'cells that share each vertex'}, + 'edgesOnVertex': {'long_name': 'edges that share each vertex'}, + # coordinates + 'xCell': {'long_name': 'x coordinate of cell centers', 'units': 'm'}, + 'yCell': {'long_name': 'y coordinate of cell centers', 'units': 'm'}, + 'zCell': {'long_name': 'z coordinate of cell centers', 'units': 'm'}, + 'latCell': { + 'long_name': 'latitude of cell centers', + 'units': 'radians', + 'standard_name': 'latitude', + }, + 'lonCell': { + 'long_name': 'longitude of cell centers', + 'units': 'radians', + 'standard_name': 'longitude', + }, + 'xEdge': {'long_name': 'x coordinate of edge midpoints', 'units': 'm'}, + 'yEdge': {'long_name': 'y coordinate of edge midpoints', 'units': 'm'}, + 'zEdge': {'long_name': 'z coordinate of edge midpoints', 'units': 'm'}, + 'latEdge': { + 'long_name': 'latitude of edge midpoints', + 'units': 'radians', + 'standard_name': 'latitude', + }, + 'lonEdge': { + 'long_name': 'longitude of edge midpoints', + 'units': 'radians', + 'standard_name': 'longitude', + }, + 'xVertex': {'long_name': 'x coordinate of vertices', 'units': 'm'}, + 'yVertex': {'long_name': 'y coordinate of vertices', 'units': 'm'}, + 'zVertex': {'long_name': 'z coordinate of vertices', 'units': 'm'}, + 'latVertex': { + 'long_name': 'latitude of vertices', + 'units': 'radians', + 'standard_name': 'latitude', + }, + 'lonVertex': { + 'long_name': 'longitude of vertices', + 'units': 'radians', + 'standard_name': 'longitude', + }, + # geometry + 'areaCell': { + 'long_name': 'area of each cell in the primal mesh', + 'units': 'm2', + 'standard_name': 'cell_area', + }, + 'areaTriangle': { + 'long_name': 'area of each triangle in the dual mesh', + 'units': 'm2', + }, + 'kiteAreasOnVertex': { + 'long_name': 'area of the part of each dual cell in each cell on the ' + 'vertex', + 'units': 'm2', + }, + 'dvEdge': { + 'long_name': 'distance between the vertices at the ends of each edge', + 'units': 'm', + }, + 'dcEdge': { + 'long_name': 'distance between the centers of the cells on each edge', + 'units': 'm', + }, + 'angleEdge': { + 'long_name': 'angle between the normal of each edge and local east', + 'units': 'radians', + }, + 'weightsOnEdge': { + 'long_name': 'weights for reconstructing tangential velocity from ' + 'edges on edge', + 'units': '1', + }, + 'meshDensity': { + 'long_name': 'value of the density function used to generate the mesh', + 'units': '1', + }, + # Coriolis + 'fCell': { + 'long_name': 'Coriolis parameter at cell centers', + 'units': 'radians s-1', + }, + 'fEdge': { + 'long_name': 'Coriolis parameter at edges', + 'units': 'radians s-1', + }, + 'fVertex': { + 'long_name': 'Coriolis parameter at vertices', + 'units': 'radians s-1', + }, + # masks + 'cullCell': { + 'long_name': 'mask of cells to be removed by the cell culler', + }, + 'boundaryVertex': { + 'long_name': 'mask of vertices with at least one inactive neighboring ' + 'cell', + }, + 'boundaryEdge': { + 'long_name': 'mask of edges with only one active neighboring cell', + }, + 'boundaryCell': { + 'long_name': 'mask of cells with at least one inactive neighboring ' + 'cell', + }, + # mesh quality from the mesh converter + 'cellQuality': { + 'long_name': 'ratio of the shortest to the longest edge of each cell', + 'units': '1', + }, + 'triangleQuality': { + 'long_name': 'ratio of the shortest to the longest edge of each dual ' + 'triangle', + 'units': '1', + }, + 'triangleAngleQuality': { + 'long_name': 'ratio of the smallest to the largest angle of each dual ' + 'triangle', + 'units': '1', + }, + 'obtuseTriangle': { + 'long_name': 'mask of dual triangles with an obtuse angle', + }, + 'gridSpacing': { + 'long_name': 'mean distance from each cell center to its neighbors', + 'units': 'm', + }, +} + + +def cf_conventions(conventions=None): + """ + Get a ``Conventions`` attribute for an MPAS mesh file that includes CF + + Entries in ``conventions`` (including any CF version) are kept, and + ``CF-1.8`` and ``MPAS`` are added if they are missing. + + Parameters + ---------- + conventions : str, optional + The existing ``Conventions`` attribute, with entries separated by + blanks or commas + + Returns + ------- + conventions : str + The ``Conventions`` attribute with a CF entry first + """ + if conventions is None: + conventions = '' + entries = conventions.replace(',', ' ').split() + cf = [entry for entry in entries if entry.startswith('CF-')] + other = [entry for entry in entries if not entry.startswith('CF-')] + if 'MPAS' not in other: + other.append('MPAS') + cf_entry = cf[0] if cf else CF_VERSION + return ' '.join([cf_entry] + other) + + +def add_mesh_attrs(ds): + """ + Add CF metadata to the mesh variables in a dataset and add CF to its + ``Conventions``. Attributes that a variable already has are kept. + + Parameters + ---------- + ds : xarray.Dataset + An MPAS mesh dataset, modified in place + + Returns + ------- + ds : xarray.Dataset + The same dataset, for convenience + """ + for name, attrs in MESH_VAR_ATTRS.items(): + if name not in ds: + continue + var_attrs = ds[name].attrs + for key, value in attrs.items(): + var_attrs.setdefault(key, value) + ds.attrs['Conventions'] = cf_conventions(ds.attrs.get('Conventions')) + return ds From d76892fb612fc19a715793540315cdeb5be5f09b Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:54:57 -0500 Subject: [PATCH 03/10] Add CF helpers to the mesh conversion tools ncutil::def_var() takes optional units and standard_name attributes, and cf_conventions() builds a Conventions attribute with a CF entry from an input one, matching mpas_tools.mesh.attrs.cf_conventions(). Co-Authored-By: Claude Opus 5.5 (1M context) --- .../netcdf_utils.h | 22 +++++++++- .../string_utils.h | 41 +++++++++++++++++++ 2 files changed, 62 insertions(+), 1 deletion(-) diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/netcdf_utils.h b/mesh_tools/mesh_conversion_tools_netcdf_c/netcdf_utils.h index add28c81b..91c8fbae6 100755 --- a/mesh_tools/mesh_conversion_tools_netcdf_c/netcdf_utils.h +++ b/mesh_tools/mesh_conversion_tools_netcdf_c/netcdf_utils.h @@ -286,7 +286,9 @@ std::string const&_name, // name of variable nc_type _type, // NetCDF data-type std::string const&_long, // NetCDF long_name - std::initializer_list _dims // dim. name list + std::initializer_list _dims, // dim. name list + std::string const&_units = "", // CF units; none if empty + std::string const&_std_name = "" // CF standard_name; none if empty ) { int _retv, _ncid, _vtag, _dtag[256]; @@ -333,6 +335,24 @@ _name + ": " + std::to_string(_retv)); } + if (!_units.empty() && (_retv = nc_put_att_text(_ncid, _vtag, + "units", _units.size(), _units.c_str()))) + { + nc_close(_ncid) ; + throw std::invalid_argument( + "Error putting variable " + + _name + ": " + std::to_string(_retv)); + } + + if (!_std_name.empty() && (_retv = nc_put_att_text(_ncid, _vtag, + "standard_name", _std_name.size(), _std_name.c_str()))) + { + nc_close(_ncid) ; + throw std::invalid_argument( + "Error putting variable " + + _name + ": " + std::to_string(_retv)); + } + if ((_retv = nc_close(_ncid))) throw std::invalid_argument( "Error handling " + diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/string_utils.h b/mesh_tools/mesh_conversion_tools_netcdf_c/string_utils.h index 39592936b..6088f8dd0 100755 --- a/mesh_tools/mesh_conversion_tools_netcdf_c/string_utils.h +++ b/mesh_tools/mesh_conversion_tools_netcdf_c/string_utils.h @@ -1,5 +1,6 @@ # include +# include # pragma once @@ -101,4 +102,44 @@ _fext = std::string(_pos4, _pos5); } + /* + -------------------------------------------------------- + * CF-CONVENTIONS: the Conventions attribute for output, + * keeping the entries (including any CF version) from + * the input and adding "CF-1.8" and "MPAS" if missing + -------------------------------------------------------- + */ + + inline std::string cf_conventions ( + std::string const& _in // input Conventions, may be empty + ) + { + std::string _list = _in, _item, _cf, _rest; + bool _mpas = false; + + // entries are separated by blanks or commas; also drop the + // trailing null from reading the attribute + for (auto &_char : _list) + { + if (_char == ',' || _char == '\0') _char = ' '; + } + + std::istringstream _stream(_list); + while (_stream >> _item) + { + if (_item.compare(0, 3, "CF-") == 0) + { + if (_cf.empty()) _cf = _item; + continue; + } + if (_item == "MPAS") _mpas = true; + _rest += " " + _item; + } + + if (_cf.empty()) _cf = "CF-1.8"; + if (!_mpas) _rest += " MPAS"; + + return _cf + _rest; + } + # endif //__STRING_UTILS__ From 894029861e5a681429f77a8801dcd7e61e6f14e3 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:55:35 -0500 Subject: [PATCH 04/10] Write CF Conventions from the mesh converter and culler MpasMeshConverter.x and MpasCellCuller.x wrote Conventions = "MPAS", which CF checkers report as an error. They now write "CF-1.8 MPAS", keeping the entries (including any CF version) of the input's Conventions. Also fix the debug message for reading mesh_spec in the culler. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../mpas_cell_culler.cpp | 14 ++++++++++++-- .../mpas_mesh_converter.cpp | 12 +++++++++++- 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp index 3ac3859a1..4138bbeca 100755 --- a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp +++ b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp @@ -26,6 +26,7 @@ double sphere_radius, xPeriod, yPeriod; string in_history = ""; string in_file_id = ""; string in_parent_id = ""; +string in_conventions = ""; string in_mesh_spec = "1.0"; bool outputMap = false; @@ -357,6 +358,14 @@ int readGridInput(const string inputFilename){/*{{{*/ // allow errors for optional attr. not found } try { +#ifdef _DEBUG + cout << " Reading Conventions" << endl; +#endif + ncutil::get_str(inputFilename, "Conventions", in_conventions); + } catch (...) { + // allow errors for optional attr. not found + } + try { #ifdef _DEBUG cout << " Reading parent_id" << endl; #endif @@ -366,7 +375,7 @@ int readGridInput(const string inputFilename){/*{{{*/ } try { #ifdef _DEBUG - cout << " Reading parent_id" << endl; + cout << " Reading mesh_spec" << endl; #endif ncutil::get_str(inputFilename, "mesh_spec", in_mesh_spec); } catch (...) { @@ -771,7 +780,8 @@ int outputGridAttributes( const string inputFilename, const string outputFilenam ncutil::put_str(outputFilename, "history", history_str); ncutil::put_str(outputFilename, "mesh_spec", in_mesh_spec); - ncutil::put_str(outputFilename, "Conventions", "MPAS"); + ncutil::put_str(outputFilename, "Conventions", + cf_conventions(in_conventions)); ncutil::put_str(outputFilename, "source", "MpasCellCuller.x"); ncutil::put_str(outputFilename, "file_id", id_str); diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp index 6f38e7424..f6e84bf17 100755 --- a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp +++ b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp @@ -40,6 +40,7 @@ double xPeriodicFix, yPeriodicFix; string in_history = ""; string in_file_id = ""; string in_parent_id = ""; +string in_conventions = ""; // Connectivity and location information {{{ @@ -357,6 +358,14 @@ int readGridInput(const string inputFilename){/*{{{*/ // allow errors for optional attr. not found } try { +#ifdef _DEBUG + cout << " Reading Conventions" << endl; +#endif + ncutil::get_str(inputFilename, "Conventions", in_conventions); + } catch (...) { + // allow errors for optional attr. not found + } + try { #ifdef _DEBUG cout << " Reading parent_id" << endl; #endif @@ -2512,7 +2521,8 @@ int outputGridAttributes( const string outputFilename, const string inputFilenam ncutil::put_str(outputFilename, "history", history_str); ncutil::put_str(outputFilename, "mesh_spec", mesh_spec_str); - ncutil::put_str(outputFilename, "Conventions", "MPAS"); + ncutil::put_str(outputFilename, "Conventions", + cf_conventions(in_conventions)); ncutil::put_str(outputFilename, "source", "MpasMeshConverter.x"); ncutil::put_str(outputFilename, "file_id", id_str); From 6b7c1b896eddc3e39a30bdd2f628e06e3de83585 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:55:43 -0500 Subject: [PATCH 05/10] Add CF attributes to mesh converter and culler variables Every variable that MpasMeshConverter.x and MpasCellCuller.x write now has the long_name, units and standard_name from mpas_tools.mesh.attrs. Several long names were wrong or unclear (e.g. cellsOnVertex was "vertices adj. to each vertex", both triangle qualities were "quality of mesh dual cells" and the culler's dvEdge was "length of arc between centres"), and none of the variables had units. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../mpas_cell_culler.cpp | 89 +++++++++------ .../mpas_mesh_converter.cpp | 107 +++++++++++------- 2 files changed, 115 insertions(+), 81 deletions(-) diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp index 4138bbeca..85035c63a 100755 --- a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp +++ b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp @@ -842,26 +842,28 @@ int mapAndOutputGridCoordinates( const string inputFilename, const string output } ncutil::def_var(outputFilename, "latCell", - NC_DOUBLE, "latitudes of cell centres", {"nCells"}); + NC_DOUBLE, "latitude of cell centers", + {"nCells"}, "radians", "latitude"); ncutil::def_var(outputFilename, "lonCell", - NC_DOUBLE, "longitudes of cell centres", {"nCells"}); + NC_DOUBLE, "longitude of cell centers", + {"nCells"}, "radians", "longitude"); ncutil::put_var(outputFilename, "latCell", &latNew[0]); ncutil::put_var(outputFilename, "lonCell", &lonNew[0]); ncutil::def_var(outputFilename, "xCell", - NC_DOUBLE, "x-coordinates of cell centres", {"nCells"}); + NC_DOUBLE, "x coordinate of cell centers", {"nCells"}, "m"); ncutil::def_var(outputFilename, "yCell", - NC_DOUBLE, "y-coordinates of cell centres", {"nCells"}); + NC_DOUBLE, "y coordinate of cell centers", {"nCells"}, "m"); ncutil::def_var(outputFilename, "zCell", - NC_DOUBLE, "z-coordinates of cell centres", {"nCells"}); + NC_DOUBLE, "z coordinate of cell centers", {"nCells"}, "m"); ncutil::put_var(outputFilename, "xCell", &xNew[0]); ncutil::put_var(outputFilename, "yCell", &yNew[0]); ncutil::put_var(outputFilename, "zCell", &zNew[0]); ncutil::def_var(outputFilename, "indexToCellID", - NC_INT, "index to cell ID mapping", {"nCells"}); + NC_INT, "global index of each cell", {"nCells"}); ncutil::put_var(outputFilename, "indexToCellID", &idxToNew[0]); @@ -912,26 +914,28 @@ int mapAndOutputGridCoordinates( const string inputFilename, const string output } ncutil::def_var(outputFilename, "latEdge", - NC_DOUBLE, "latitudes of edge centres", {"nEdges"}); + NC_DOUBLE, "latitude of edge midpoints", + {"nEdges"}, "radians", "latitude"); ncutil::def_var(outputFilename, "lonEdge", - NC_DOUBLE, "longitudes of edge centres", {"nEdges"}); + NC_DOUBLE, "longitude of edge midpoints", + {"nEdges"}, "radians", "longitude"); ncutil::put_var(outputFilename, "latEdge", &latNew[0]); ncutil::put_var(outputFilename, "lonEdge", &lonNew[0]); ncutil::def_var(outputFilename, "xEdge", - NC_DOUBLE, "x-coordinates of edge centres", {"nEdges"}); + NC_DOUBLE, "x coordinate of edge midpoints", {"nEdges"}, "m"); ncutil::def_var(outputFilename, "yEdge", - NC_DOUBLE, "y-coordinates of edge centres", {"nEdges"}); + NC_DOUBLE, "y coordinate of edge midpoints", {"nEdges"}, "m"); ncutil::def_var(outputFilename, "zEdge", - NC_DOUBLE, "z-coordinates of edge centres", {"nEdges"}); + NC_DOUBLE, "z coordinate of edge midpoints", {"nEdges"}, "m"); ncutil::put_var(outputFilename, "xEdge", &xNew[0]); ncutil::put_var(outputFilename, "yEdge", &yNew[0]); ncutil::put_var(outputFilename, "zEdge", &zNew[0]); ncutil::def_var(outputFilename, "indexToEdgeID", - NC_INT, "index to edge ID mapping", {"nEdges"}); + NC_INT, "global index of each edge", {"nEdges"}); ncutil::put_var(outputFilename, "indexToEdgeID", &idxToNew[0]); @@ -982,26 +986,28 @@ int mapAndOutputGridCoordinates( const string inputFilename, const string output } ncutil::def_var(outputFilename, "latVertex", - NC_DOUBLE, "latitudes of vertices", {"nVertices"}); + NC_DOUBLE, "latitude of vertices", + {"nVertices"}, "radians", "latitude"); ncutil::def_var(outputFilename, "lonVertex", - NC_DOUBLE, "longitudes of vertices", {"nVertices"}); + NC_DOUBLE, "longitude of vertices", + {"nVertices"}, "radians", "longitude"); ncutil::put_var(outputFilename, "latVertex", &latNew[0]); ncutil::put_var(outputFilename, "lonVertex", &lonNew[0]); ncutil::def_var(outputFilename, "xVertex", - NC_DOUBLE, "x-coordinates of vertices", {"nVertices"}); + NC_DOUBLE, "x coordinate of vertices", {"nVertices"}, "m"); ncutil::def_var(outputFilename, "yVertex", - NC_DOUBLE, "y-coordinates of vertices", {"nVertices"}); + NC_DOUBLE, "y coordinate of vertices", {"nVertices"}, "m"); ncutil::def_var(outputFilename, "zVertex", - NC_DOUBLE, "z-coordinates of vertices", {"nVertices"}); + NC_DOUBLE, "z coordinate of vertices", {"nVertices"}, "m"); ncutil::put_var(outputFilename, "xVertex", &xNew[0]); ncutil::put_var(outputFilename, "yVertex", &yNew[0]); ncutil::put_var(outputFilename, "zVertex", &zNew[0]); ncutil::def_var(outputFilename, "indexToVertexID", - NC_INT, "index to vertex ID mapping", {"nVertices"}); + NC_INT, "global index of each vertex", {"nVertices"}); ncutil::put_var(outputFilename, "indexToVertexID", &idxToNew[0]); @@ -1070,7 +1076,7 @@ int mapAndOutputCellFields( const string inputFilename, const string outputPath, // Write nEdgesOncell to output file ncutil::def_var(outputFilename, "nEdgesOnCell", - NC_INT, "number of edges on each cell", {"nCells"}); + NC_INT, "number of edges that border each cell", {"nCells"}); ncutil::put_var(outputFilename, "nEdgesOnCell", &nEdgesOnCellNew[0]); @@ -1099,7 +1105,7 @@ int mapAndOutputCellFields( const string inputFilename, const string outputPath, } ncutil::def_var(outputFilename, "edgesOnCell", - NC_INT, "edges on each cell", {"nCells", "maxEdges"}); + NC_INT, "edges that border each cell", {"nCells", "maxEdges"}); ncutil::put_var(outputFilename, "edgesOnCell", &tmp_arr_new[0]); @@ -1137,7 +1143,7 @@ int mapAndOutputCellFields( const string inputFilename, const string outputPath, graph.close(); ncutil::def_var(outputFilename, "cellsOnCell", - NC_INT, "cells adj. to each cell", {"nCells", "maxEdges"}); + NC_INT, "cells that neighbor each cell", {"nCells", "maxEdges"}); ncutil::put_var(outputFilename, "cellsOnCell", &tmp_arr_new[0]); @@ -1162,7 +1168,7 @@ int mapAndOutputCellFields( const string inputFilename, const string outputPath, } ncutil::def_var(outputFilename, "verticesOnCell", - NC_INT, "vertices on each cell", {"nCells", "maxEdges"}); + NC_INT, "vertices that border each cell", {"nCells", "maxEdges"}); ncutil::put_var(outputFilename, "verticesOnCell", &tmp_arr_new[0]); @@ -1179,7 +1185,8 @@ int mapAndOutputCellFields( const string inputFilename, const string outputPath, } ncutil::def_var(outputFilename, "areaCell", - NC_DOUBLE, "surface area of each cell", {"nCells"}); + NC_DOUBLE, "area of each cell in the primal mesh", + {"nCells"}, "m2", "cell_area"); ncutil::put_var(outputFilename, "areaCell", &areaCellNew[0]); @@ -1205,7 +1212,8 @@ int mapAndOutputCellFields( const string inputFilename, const string outputPath, } ncutil::def_var(outputFilename, "meshDensity", - NC_DOUBLE, "mesh density distribution", {"nCells"}); + NC_DOUBLE, "value of the density function used to generate the mesh", + {"nCells"}, "1"); ncutil::put_var(outputFilename, "meshDensity", &meshDensityNew[0]); delete[] meshDensityNew; @@ -1322,9 +1330,9 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen } ncutil::def_var(outputFilename, "verticesOnEdge", - NC_INT, "vertices on each edge", {"nEdges", "TWO"}); + NC_INT, "vertices that straddle each edge", {"nEdges", "TWO"}); ncutil::def_var(outputFilename, "cellsOnEdge", - NC_INT, "cells adj. to each edge", {"nEdges", "TWO"}); + NC_INT, "cells that straddle each edge", {"nEdges", "TWO"}); ncutil::put_var(outputFilename, "verticesOnEdge", &verticesOnEdgeNew[0]); ncutil::put_var(outputFilename, "cellsOnEdge", &cellsOnEdgeNew[0]); @@ -1406,16 +1414,19 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen } ncutil::def_var(outputFilename, "nEdgesOnEdge", - NC_INT, "number of edges adj. to each edge", {"nEdges"}); + NC_INT, "number of edges that border the cells that straddle each edge", + {"nEdges"}); ncutil::def_var(outputFilename, "edgesOnEdge", - NC_INT, "edges adj. to each edge", {"nEdges", "maxEdges2"}); + NC_INT, "edges that border the cells that straddle each edge", + {"nEdges", "maxEdges2"}); ncutil::put_var(outputFilename, "nEdgesOnEdge", &nEdgesOnEdgeNew[0]); ncutil::put_var(outputFilename, "edgesOnEdge", &edgesOnEdgeNew[0]); if(hasWeightsOnEdge) { ncutil::def_var(outputFilename, "weightsOnEdge", - NC_DOUBLE, "tangential flux reconstruction weights", {"nEdges", "maxEdges2"}); + NC_DOUBLE, "weights for reconstructing tangential velocity from edges on edge", + {"nEdges", "maxEdges2"}, "1"); ncutil::put_var(outputFilename, "weightsOnEdge", &weightsOnEdgeNew[0]); } @@ -1457,11 +1468,14 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen } ncutil::def_var(outputFilename, "dvEdge", - NC_DOUBLE, "length of arc between centres", {"nEdges"}); + NC_DOUBLE, "distance between the vertices at the ends of each edge", + {"nEdges"}, "m"); ncutil::def_var(outputFilename, "dcEdge", - NC_DOUBLE, "length of arc between centres", {"nEdges"}); + NC_DOUBLE, "distance between the centers of the cells on each edge", + {"nEdges"}, "m"); ncutil::def_var(outputFilename, "angleEdge", - NC_DOUBLE, "angle to edges", {"nEdges"}) ; + NC_DOUBLE, "angle between the normal of each edge and local east", + {"nEdges"}, "radians"); ncutil::put_var(outputFilename, "dvEdge", &dvEdgeNew[0]); ncutil::put_var(outputFilename, "dcEdge", &dcEdgeNew[0]); @@ -1546,18 +1560,19 @@ int mapAndOutputVertexFields( const string inputFilename, const string outputFil } ncutil::def_var(outputFilename, "edgesOnVertex", - NC_INT, "edges adj. to each vertex", {"nVertices", "vertexDegree"}); + NC_INT, "edges that share each vertex", {"nVertices", "vertexDegree"}); ncutil::def_var(outputFilename, "cellsOnVertex", - NC_INT, "cells adj. to each vertex", {"nVertices", "vertexDegree"}); + NC_INT, "cells that share each vertex", {"nVertices", "vertexDegree"}); ncutil::put_var(outputFilename, "edgesOnVertex", &edgesOnVertexNew [0]); ncutil::put_var(outputFilename, "cellsOnVertex", &cellsOnVertexNew [0]); ncutil::def_var(outputFilename, "areaTriangle", - NC_DOUBLE, "surface area of dual cells", {"nVertices"}); + NC_DOUBLE, "area of each triangle in the dual mesh", + {"nVertices"}, "m2"); ncutil::def_var(outputFilename, "kiteAreasOnVertex", - NC_DOUBLE, - "surface areas of overlap between cells and dual cells", {"nVertices", "vertexDegree"}); + NC_DOUBLE, "area of the part of each dual cell in each cell on the vertex", + {"nVertices", "vertexDegree"}, "m2"); ncutil::put_var(outputFilename, "areaTriangle", &areaTriangleNew [0]); ncutil::put_var(outputFilename, "kiteAreasOnVertex", &kiteAreasOnVertexNew [0]); diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp index f6e84bf17..1e93a4753 100755 --- a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp +++ b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp @@ -2574,26 +2574,28 @@ int outputGridCoordinates( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "latCell", - NC_DOUBLE, "latitudes of cell centres", {"nCells"}); + NC_DOUBLE, "latitude of cell centers", + {"nCells"}, "radians", "latitude"); ncutil::def_var(outputFilename, "lonCell", - NC_DOUBLE, "longitudes of cell centres", {"nCells"}); + NC_DOUBLE, "longitude of cell centers", + {"nCells"}, "radians", "longitude"); ncutil::put_var(outputFilename, "latCell", &lat[0]); ncutil::put_var(outputFilename, "lonCell", &lon[0]); ncutil::def_var(outputFilename, "xCell", - NC_DOUBLE, "x-coordinates of cell centres", {"nCells"}); + NC_DOUBLE, "x coordinate of cell centers", {"nCells"}, "m"); ncutil::def_var(outputFilename, "yCell", - NC_DOUBLE, "y-coordinates of cell centres", {"nCells"}); + NC_DOUBLE, "y coordinate of cell centers", {"nCells"}, "m"); ncutil::def_var(outputFilename, "zCell", - NC_DOUBLE, "z-coordinates of cell centres", {"nCells"}); + NC_DOUBLE, "z coordinate of cell centers", {"nCells"}, "m"); ncutil::put_var(outputFilename, "xCell", &x[0]); ncutil::put_var(outputFilename, "yCell", &y[0]); ncutil::put_var(outputFilename, "zCell", &z[0]); ncutil::def_var(outputFilename, "indexToCellID", - NC_INT, "index to cell ID mapping", {"nCells"}); + NC_INT, "global index of each cell", {"nCells"}); ncutil::put_var(outputFilename, "indexToCellID", &idxTo[0]); @@ -2633,26 +2635,28 @@ int outputGridCoordinates( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "latEdge", - NC_DOUBLE, "latitudes of edge centres", {"nEdges"}); + NC_DOUBLE, "latitude of edge midpoints", + {"nEdges"}, "radians", "latitude"); ncutil::def_var(outputFilename, "lonEdge", - NC_DOUBLE, "longitudes of edge centres", {"nEdges"}); + NC_DOUBLE, "longitude of edge midpoints", + {"nEdges"}, "radians", "longitude"); ncutil::put_var(outputFilename, "latEdge", &lat[0]); ncutil::put_var(outputFilename, "lonEdge", &lon[0]); ncutil::def_var(outputFilename, "xEdge", - NC_DOUBLE, "x-coordinates of edge centres", {"nEdges"}); + NC_DOUBLE, "x coordinate of edge midpoints", {"nEdges"}, "m"); ncutil::def_var(outputFilename, "yEdge", - NC_DOUBLE, "y-coordinates of edge centres", {"nEdges"}); + NC_DOUBLE, "y coordinate of edge midpoints", {"nEdges"}, "m"); ncutil::def_var(outputFilename, "zEdge", - NC_DOUBLE, "z-coordinates of edge centres", {"nEdges"}); + NC_DOUBLE, "z coordinate of edge midpoints", {"nEdges"}, "m"); ncutil::put_var(outputFilename, "xEdge", &x[0]); ncutil::put_var(outputFilename, "yEdge", &y[0]); ncutil::put_var(outputFilename, "zEdge", &z[0]); ncutil::def_var(outputFilename, "indexToEdgeID", - NC_INT, "index to edge ID mapping", {"nEdges"}); + NC_INT, "global index of each edge", {"nEdges"}); ncutil::put_var(outputFilename, "indexToEdgeID", &idxTo[0]); @@ -2692,26 +2696,28 @@ int outputGridCoordinates( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "latVertex", - NC_DOUBLE, "latitudes of vertices", {"nVertices"}); + NC_DOUBLE, "latitude of vertices", + {"nVertices"}, "radians", "latitude"); ncutil::def_var(outputFilename, "lonVertex", - NC_DOUBLE, "longitudes of vertices", {"nVertices"}); + NC_DOUBLE, "longitude of vertices", + {"nVertices"}, "radians", "longitude"); ncutil::put_var(outputFilename, "latVertex", &lat[0]); ncutil::put_var(outputFilename, "lonVertex", &lon[0]); ncutil::def_var(outputFilename, "xVertex", - NC_DOUBLE, "x-coordinates of vertices", {"nVertices"}); + NC_DOUBLE, "x coordinate of vertices", {"nVertices"}, "m"); ncutil::def_var(outputFilename, "yVertex", - NC_DOUBLE, "y-coordinates of vertices", {"nVertices"}); + NC_DOUBLE, "y coordinate of vertices", {"nVertices"}, "m"); ncutil::def_var(outputFilename, "zVertex", - NC_DOUBLE, "z-coordinates of vertices", {"nVertices"}); + NC_DOUBLE, "z coordinate of vertices", {"nVertices"}, "m"); ncutil::put_var(outputFilename, "xVertex", &x[0]); ncutil::put_var(outputFilename, "yVertex", &y[0]); ncutil::put_var(outputFilename, "zVertex", &z[0]); ncutil::def_var(outputFilename, "indexToVertexID", - NC_INT, "index to vertex ID mapping", {"nVertices"}); + NC_INT, "global index of each vertex", {"nVertices"}); ncutil::put_var(outputFilename, "indexToVertexID", &idxTo[0]); @@ -2761,7 +2767,7 @@ int outputCellConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "cellsOnCell", - NC_INT, "cells adj. to each cell", {"nCells", "maxEdges"}); + NC_INT, "cells that neighbor each cell", {"nCells", "maxEdges"}); ncutil::put_var(outputFilename, "cellsOnCell", &tmp_arr[0]); @@ -2784,7 +2790,7 @@ int outputCellConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "edgesOnCell", - NC_INT, "edges on each cell", {"nCells", "maxEdges"}); + NC_INT, "edges that border each cell", {"nCells", "maxEdges"}); ncutil::put_var(outputFilename, "edgesOnCell", &tmp_arr[0]); @@ -2806,7 +2812,7 @@ int outputCellConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "verticesOnCell", - NC_INT, "vertices on each cell", {"nCells", "maxEdges"}); + NC_INT, "vertices that border each cell", {"nCells", "maxEdges"}); ncutil::put_var(outputFilename, "verticesOnCell", &tmp_arr[0]); @@ -2822,7 +2828,7 @@ int outputCellConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "nEdgesOnCell", - NC_INT, "number of edges on each cell", {"nCells"}); + NC_INT, "number of edges that border each cell", {"nCells"}); ncutil::put_var(outputFilename, "nEdgesOnCell", &tmp_arr[0]); @@ -2874,7 +2880,8 @@ int outputEdgeConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "edgesOnEdge", - NC_INT, "edges adj. to each edge", {"nEdges", "maxEdges2"}); + NC_INT, "edges that border the cells that straddle each edge", + {"nEdges", "maxEdges2"}); ncutil::put_var(outputFilename, "edgesOnEdge", &tmp_arr[0]); @@ -2900,7 +2907,7 @@ int outputEdgeConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "cellsOnEdge", - NC_INT, "cells adj. to each edge", {"nEdges", "TWO"}); + NC_INT, "cells that straddle each edge", {"nEdges", "TWO"}); ncutil::put_var(outputFilename, "cellsOnEdge", &tmp_arr[0]); @@ -2917,7 +2924,7 @@ int outputEdgeConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "verticesOnEdge", - NC_INT, "vertices on each edge", {"nEdges", "TWO"}); + NC_INT, "vertices that straddle each edge", {"nEdges", "TWO"}); ncutil::put_var(outputFilename, "verticesOnEdge", &tmp_arr[0]); @@ -2932,7 +2939,8 @@ int outputEdgeConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "nEdgesOnEdge", - NC_INT, "number of edges on each edge", {"nEdges"}); + NC_INT, "number of edges that border the cells that straddle each edge", + {"nEdges"}); ncutil::put_var(outputFilename, "nEdgesOnEdge", &tmp_arr[0]); @@ -2980,7 +2988,7 @@ int outputVertexConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "cellsOnVertex", - NC_INT, "vertices adj. to each vertex", {"nVertices", "vertexDegree"}); + NC_INT, "cells that share each vertex", {"nVertices", "vertexDegree"}); ncutil::put_var(outputFilename, "cellsOnVertex", &tmp_arr[0]); @@ -3002,7 +3010,7 @@ int outputVertexConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "edgesOnVertex", - NC_INT, "edges adj. to each vertex", {"nVertices", "vertexDegree"}); + NC_INT, "edges that share each vertex", {"nVertices", "vertexDegree"}); ncutil::put_var(outputFilename, "edgesOnVertex", &tmp_arr[0]); @@ -3022,7 +3030,8 @@ int outputVertexConnectivity( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "boundaryVertex", - NC_INT, "non-zero for each vertex on mesh boundary", {"nVertices"}); + NC_INT, "mask of vertices with at least one inactive neighboring cell", + {"nVertices"}); ncutil::put_var(outputFilename, "boundaryVertex", &tmp_arr[0]); @@ -3051,7 +3060,8 @@ int outputCellParameters( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "areaCell", - NC_DOUBLE, "surface areas of cells", {"nCells"}); + NC_DOUBLE, "area of each cell in the primal mesh", + {"nCells"}, "m2", "cell_area"); ncutil::put_var(outputFilename, "areaCell", &areaCell[0]); @@ -3082,7 +3092,8 @@ int outputVertexParameters( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "areaTriangle", - NC_DOUBLE, "surface areas of dual cells", {"nVertices"}); + NC_DOUBLE, "area of each triangle in the dual mesh", + {"nVertices"}, "m2"); ncutil::put_var(outputFilename, "areaTriangle", &areaTriangle[0]); @@ -3110,8 +3121,8 @@ int outputVertexParameters( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "kiteAreasOnVertex", - NC_DOUBLE, - "surface areas of overlap between cells and dual cells", {"nVertices", "vertexDegree"}); + NC_DOUBLE, "area of the part of each dual cell in each cell on the vertex", + {"nVertices", "vertexDegree"}, "m2"); ncutil::put_var(outputFilename, "kiteAreasOnVertex", &tmp_arr[0]); @@ -3146,14 +3157,17 @@ int outputEdgeParameters( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "angleEdge", - NC_DOUBLE, "angle to edges", {"nEdges"}) ; + NC_DOUBLE, "angle between the normal of each edge and local east", + {"nEdges"}, "radians"); ncutil::put_var(outputFilename, "angleEdge", &angleEdge[0]); ncutil::def_var(outputFilename, "dcEdge", - NC_DOUBLE, "length of arc between centres", {"nEdges"}); + NC_DOUBLE, "distance between the centers of the cells on each edge", + {"nEdges"}, "m"); ncutil::def_var(outputFilename, "dvEdge", - NC_DOUBLE, "length of arc between vertices", {"nEdges"}); + NC_DOUBLE, "distance between the vertices at the ends of each edge", + {"nEdges"}, "m"); ncutil::put_var(outputFilename, "dcEdge", &dcEdge[0]) ; ncutil::put_var(outputFilename, "dvEdge", &dvEdge[0]) ; @@ -3177,7 +3191,8 @@ int outputEdgeParameters( const string outputFilename) {/*{{{*/ } ncutil::def_var(outputFilename, "weightsOnEdge", - NC_DOUBLE, "tangential flux reconstruction weights", {"nEdges", "maxEdges2"}); + NC_DOUBLE, "weights for reconstructing tangential velocity from edges on edge", + {"nEdges", "maxEdges2"}, "1"); ncutil::put_var(outputFilename, "weightsOnEdge", &tmp_arr[0]); @@ -3198,7 +3213,8 @@ int outputMeshDensity( const string outputFilename) {/*{{{*/ * *************************************************************************/ ncutil::def_var(outputFilename, "meshDensity", - NC_DOUBLE, "mesh density distribution", {"nCells"}); + NC_DOUBLE, "value of the density function used to generate the mesh", + {"nCells"}, "1"); ncutil::put_var(outputFilename, "meshDensity", &meshDensity[0]); @@ -3217,29 +3233,32 @@ int outputMeshQualities( const string outputFilename) {/*{{{*/ * *************************************************************************/ ncutil::def_var(outputFilename, "cellQuality", - NC_DOUBLE, "quality of mesh cells", {"nCells"}); + NC_DOUBLE, "ratio of the shortest to the longest edge of each cell", + {"nCells"}, "1"); ncutil::put_var(outputFilename, "cellQuality", &cellQuality[0]); ncutil::def_var(outputFilename, "gridSpacing", - NC_DOUBLE, "grid spacing distribution", {"nCells"}); + NC_DOUBLE, "mean distance from each cell center to its neighbors", + {"nCells"}, "m"); ncutil::put_var(outputFilename, "gridSpacing", &cellQuality[0]); ncutil::def_var(outputFilename, "triangleQuality", - NC_DOUBLE, "quality of mesh dual cells", {"nVertices"}); + NC_DOUBLE, "ratio of the shortest to the longest edge of each dual triangle", + {"nVertices"}, "1"); ncutil::put_var(outputFilename, "triangleQuality", &triangleQuality[0]); ncutil::def_var(outputFilename, "triangleAngleQuality", - NC_DOUBLE, "quality of mesh dual cells", {"nVertices"}); + NC_DOUBLE, "ratio of the smallest to the largest angle of each dual triangle", + {"nVertices"}, "1"); ncutil::put_var(outputFilename, "triangleAngleQuality", &triangleAngleQuality [0]); ncutil::def_var(outputFilename, "obtuseTriangle", - NC_INT, - "non-zero for any dual cell containing obtuse angles", {"nVertices"}); + NC_INT, "mask of dual triangles with an obtuse angle", {"nVertices"}); ncutil::put_var(outputFilename, "obtuseTriangle", &obtuseTriangle[0]); From 8bf8202b719bf1fbb75ecb2615c56e54fba91002 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:56:03 -0500 Subject: [PATCH 06/10] Add CF metadata to meshes from Python mesh creation tools make_planar_hex_mesh(), jigsaw_to_netcdf() and triangle_to_netcdf() wrote no units or long_name on the mesh variables and no Conventions. They now add the metadata from mpas_tools.mesh.attrs. Co-Authored-By: Claude Opus 5.5 (1M context) --- conda_package/mpas_tools/mesh/creation/jigsaw_to_netcdf.py | 3 +++ conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py | 5 +++++ conda_package/mpas_tools/planar_hex.py | 3 +++ 3 files changed, 11 insertions(+) diff --git a/conda_package/mpas_tools/mesh/creation/jigsaw_to_netcdf.py b/conda_package/mpas_tools/mesh/creation/jigsaw_to_netcdf.py index 75664e72d..894712eb0 100644 --- a/conda_package/mpas_tools/mesh/creation/jigsaw_to_netcdf.py +++ b/conda_package/mpas_tools/mesh/creation/jigsaw_to_netcdf.py @@ -4,6 +4,7 @@ import xarray as xr from mpas_tools.io import write_netcdf +from mpas_tools.mesh.attrs import add_mesh_attrs from mpas_tools.mesh.creation.open_msh import readmsh from mpas_tools.mesh.creation.util import circumcenter @@ -108,6 +109,8 @@ def jigsaw_to_netcdf(msh_filename, output_name, on_sphere, sphere_radius=None): attrs=attrs, ) + add_mesh_attrs(ds) + # Write to NetCDF using write_netcdf write_netcdf(ds, output_name) diff --git a/conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py b/conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py index 86ba339e3..1113a27f5 100644 --- a/conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py +++ b/conda_package/mpas_tools/mesh/creation/triangle_to_netcdf.py @@ -10,6 +10,7 @@ import numpy as np from netCDF4 import Dataset as NetCDFFile +from mpas_tools.mesh.attrs import MESH_VAR_ATTRS, cf_conventions from mpas_tools.mesh.creation.util import circumcenter @@ -151,6 +152,10 @@ def triangle_to_netcdf(node, ele, output_name): ) var[:] = cellsOnVertex_full + for name, var in grid.variables.items(): + var.setncatts(MESH_VAR_ATTRS.get(name, {})) + grid.Conventions = cf_conventions() + grid.sync() grid.close() diff --git a/conda_package/mpas_tools/planar_hex.py b/conda_package/mpas_tools/planar_hex.py index 335377f81..e560973ee 100755 --- a/conda_package/mpas_tools/planar_hex.py +++ b/conda_package/mpas_tools/planar_hex.py @@ -13,6 +13,7 @@ import xarray from mpas_tools.io import write_netcdf +from mpas_tools.mesh.attrs import add_mesh_attrs def make_planar_hex_mesh( @@ -84,6 +85,8 @@ def make_planar_hex_mesh( # the hex mesh mesh = mesh.drop_vars(['cellIdx', 'cellRow', 'cellCol']) + add_mesh_attrs(mesh) + if outFileName is not None: write_netcdf(mesh, outFileName, format=format, engine=engine) From 19ddb82b6b8bb2ca3329dd68a3a83dc2922c8b17 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:56:59 -0500 Subject: [PATCH 07/10] Apply pre-commit formatting to the SCRIP modules Reformat before changing them so the functional change is easy to review. Co-Authored-By: Claude Opus 5.5 (1M context) --- conda_package/mpas_tools/scrip/from_mpas.py | 166 +++++++++++------- conda_package/mpas_tools/scrip/from_planar.py | 146 +++++++++------ 2 files changed, 196 insertions(+), 116 deletions(-) diff --git a/conda_package/mpas_tools/scrip/from_mpas.py b/conda_package/mpas_tools/scrip/from_mpas.py index f24236f7c..c1ba0c544 100755 --- a/conda_package/mpas_tools/scrip/from_mpas.py +++ b/conda_package/mpas_tools/scrip/from_mpas.py @@ -1,10 +1,11 @@ # Create a SCRIP file from an MPAS mesh. # See for details: http://www.earthsystemmodeling.org/esmf_releases/public/ESMF_5_2_0rp1/ESMF_refdoc/node3.html#SECTION03024000000000000000 +from optparse import OptionParser + import netCDF4 import numpy as np -from optparse import OptionParser from mpas_tools.cime.constants import constants @@ -24,9 +25,9 @@ def scrip_from_mpas(mpasFile, scripFile, useLandIceMask=False): Whether to use the landIceMask field for masking """ if useLandIceMask: - print(" -- Landice Masks are enabled") + print(' -- Landice Masks are enabled') else: - print(" -- Landice Masks are disabled") + print(' -- Landice Masks are disabled') # make a space in stdout before further output print('') @@ -49,20 +50,24 @@ def scrip_from_mpas(mpasFile, scripFile, useLandIceMask=False): # check the longitude convention to use positive values [0 2pi] if np.any(np.logical_or(lonCell < 0, lonCell > 2.0 * np.pi)): - raise ValueError("lonCell is not in the desired range (0, 2pi)") + raise ValueError('lonCell is not in the desired range (0, 2pi)') if np.any(np.logical_or(lonVertex < 0, lonVertex > 2.0 * np.pi)): - raise ValueError("lonVertex is not in the desired range (0, 2pi)") + raise ValueError('lonVertex is not in the desired range (0, 2pi)') if sphereRadius <= 0: sphereRadius = constants['SHR_CONST_REARTH'] - print(f" -- WARNING: sphereRadius<0 so setting sphereRadius = " - f"{constants['SHR_CONST_REARTH']}") - - if on_a_sphere == "NO": - print(" -- WARNING: 'on_a_sphere' attribute is 'NO', which means that " - "there may be some disagreement regarding area between the " - "planar (source) and spherical (target) mesh") + print( + f' -- WARNING: sphereRadius<0 so setting sphereRadius = ' + f'{constants["SHR_CONST_REARTH"]}' + ) + + if on_a_sphere == 'NO': + print( + " -- WARNING: 'on_a_sphere' attribute is 'NO', which means that " + 'there may be some disagreement regarding area between the ' + 'planar (source) and spherical (target) mesh' + ) if useLandIceMask: landIceMask = fin.variables['landIceMask'][:] @@ -71,22 +76,26 @@ def scrip_from_mpas(mpasFile, scripFile, useLandIceMask=False): # Write to output file # Dimensions - fout.createDimension("grid_size", nCells) - fout.createDimension("grid_corners", maxVertices) - fout.createDimension("grid_rank", 1) + fout.createDimension('grid_size', nCells) + fout.createDimension('grid_corners', maxVertices) + fout.createDimension('grid_rank', 1) # Variables - grid_center_lat = fout.createVariable('grid_center_lat', 'f8', - ('grid_size',)) + grid_center_lat = fout.createVariable( + 'grid_center_lat', 'f8', ('grid_size',) + ) grid_center_lat.units = 'radians' - grid_center_lon = fout.createVariable('grid_center_lon', 'f8', - ('grid_size',)) + grid_center_lon = fout.createVariable( + 'grid_center_lon', 'f8', ('grid_size',) + ) grid_center_lon.units = 'radians' - grid_corner_lat = fout.createVariable('grid_corner_lat', 'f8', - ('grid_size', 'grid_corners')) + grid_corner_lat = fout.createVariable( + 'grid_corner_lat', 'f8', ('grid_size', 'grid_corners') + ) grid_corner_lat.units = 'radians' - grid_corner_lon = fout.createVariable('grid_corner_lon', 'f8', - ('grid_size', 'grid_corners')) + grid_corner_lon = fout.createVariable( + 'grid_corner_lon', 'f8', ('grid_size', 'grid_corners') + ) grid_corner_lon.units = 'radians' grid_area = fout.createVariable('grid_area', 'f8', ('grid_size',)) grid_area.units = 'radian^2' @@ -97,20 +106,22 @@ def scrip_from_mpas(mpasFile, scripFile, useLandIceMask=False): grid_center_lat[:] = latCell[:] grid_center_lon[:] = lonCell[:] # SCRIP uses square radians - grid_area[:] = areaCell[:]/(sphereRadius**2) + grid_area[:] = areaCell[:] / (sphereRadius**2) grid_dims[:] = nCells # grid corners: grid_corner_lon_local = np.zeros((nCells, maxVertices)) grid_corner_lat_local = np.zeros((nCells, maxVertices)) cellIndices = np.arange(nCells) - lastValidVertex = verticesOnCell[cellIndices, nEdgesOnCell-1] + lastValidVertex = verticesOnCell[cellIndices, nEdgesOnCell - 1] for iVertex in range(maxVertices): mask = iVertex < nEdgesOnCell - grid_corner_lat_local[mask, iVertex] = \ - latVertex[verticesOnCell[mask, iVertex]] - grid_corner_lon_local[mask, iVertex] = \ - lonVertex[verticesOnCell[mask, iVertex]] + grid_corner_lat_local[mask, iVertex] = latVertex[ + verticesOnCell[mask, iVertex] + ] + grid_corner_lon_local[mask, iVertex] = lonVertex[ + verticesOnCell[mask, iVertex] + ] mask = iVertex >= nEdgesOnCell grid_corner_lat_local[mask, iVertex] = latVertex[lastValidVertex[mask]] @@ -126,50 +137,85 @@ def scrip_from_mpas(mpasFile, scripFile, useLandIceMask=False): grid_corner_lat[:] = grid_corner_lat_local[:] grid_corner_lon[:] = grid_corner_lon_local[:] - print("Input latCell min/max values (radians): {}, {}".format( - latCell[:].min(), latCell[:].max())) - print("Input lonCell min/max values (radians): {}, {}".format( - lonCell[:].min(), lonCell[:].max())) - print("Calculated grid_center_lat min/max values (radians): {}, {}".format( - grid_center_lat[:].min(), grid_center_lat[:].max())) - print("Calculated grid_center_lon min/max values (radians): {}, {}".format( - grid_center_lon[:].min(), grid_center_lon[:].max())) - print("Calculated grid_area min/max values (sq radians): {}, {}".format( - grid_area[:].min(), grid_area[:].max())) + print( + 'Input latCell min/max values (radians): {}, {}'.format( + latCell[:].min(), latCell[:].max() + ) + ) + print( + 'Input lonCell min/max values (radians): {}, {}'.format( + lonCell[:].min(), lonCell[:].max() + ) + ) + print( + 'Calculated grid_center_lat min/max values (radians): {}, {}'.format( + grid_center_lat[:].min(), grid_center_lat[:].max() + ) + ) + print( + 'Calculated grid_center_lon min/max values (radians): {}, {}'.format( + grid_center_lon[:].min(), grid_center_lon[:].max() + ) + ) + print( + 'Calculated grid_area min/max values (sq radians): {}, {}'.format( + grid_area[:].min(), grid_area[:].max() + ) + ) fin.close() fout.close() - print("Creation of SCRIP file is complete.") + print('Creation of SCRIP file is complete.') def main(): - print("== Gathering information. (Invoke with --help for more details. " - "All arguments are optional)") + print( + '== Gathering information. (Invoke with --help for more details. ' + 'All arguments are optional)' + ) parser = OptionParser() - parser.description = "This script takes an MPAS grid file and generates " \ - "a SCRIP grid file." - parser.add_option("-m", "--mpas", dest="mpasFile", - help="MPAS grid file name used as input.", - default="grid.nc", metavar="FILENAME") - parser.add_option("-s", "--scrip", dest="scripFile", - help="SCRIP grid file to output.", default="scrip.nc", - metavar="FILENAME") - parser.add_option("-l", "--landice", dest="landiceMasks", - help="If flag is on, landice masks will be computed " - "and used.", - action="store_true") + parser.description = ( + 'This script takes an MPAS grid file and generates a SCRIP grid file.' + ) + parser.add_option( + '-m', + '--mpas', + dest='mpasFile', + help='MPAS grid file name used as input.', + default='grid.nc', + metavar='FILENAME', + ) + parser.add_option( + '-s', + '--scrip', + dest='scripFile', + help='SCRIP grid file to output.', + default='scrip.nc', + metavar='FILENAME', + ) + parser.add_option( + '-l', + '--landice', + dest='landiceMasks', + help='If flag is on, landice masks will be computed and used.', + action='store_true', + ) for option in parser.option_list: - if option.default != ("NO", "DEFAULT"): - option.help += (" " if option.help else "") + "[default: %default]" + if option.default != ('NO', 'DEFAULT'): + option.help += (' ' if option.help else '') + '[default: %default]' options, args = parser.parse_args() if not options.mpasFile: - raise ValueError('MPAS input grid file is required. Specify with -m ' - 'command line argument.') + raise ValueError( + 'MPAS input grid file is required. Specify with -m ' + 'command line argument.' + ) if not options.scripFile: - raise ValueError('SCRIP output grid file is required. Specify with ' - '-s command line argument.') + raise ValueError( + 'SCRIP output grid file is required. Specify with ' + '-s command line argument.' + ) if not options.landiceMasks: options.landiceMasks = False diff --git a/conda_package/mpas_tools/scrip/from_planar.py b/conda_package/mpas_tools/scrip/from_planar.py index 5d62a4717..825682d87 100644 --- a/conda_package/mpas_tools/scrip/from_planar.py +++ b/conda_package/mpas_tools/scrip/from_planar.py @@ -1,11 +1,12 @@ # Create a SCRIP file from a planar rectanfular mesh. # See for details: http://www.earthsystemmodeling.org/esmf_releases/public/ESMF_5_2_0rp1/ESMF_refdoc/node3.html#SECTION03024000000000000000 -import netCDF4 -import numpy as np from optparse import OptionParser + import matplotlib.pyplot as plt -from pyproj import Transformer, CRS +import netCDF4 +import numpy as np +from pyproj import CRS, Transformer from mpas_tools.landice.projections import projections @@ -15,49 +16,77 @@ def main(): Create a SCRIP file from a planar rectanfular mesh """ - print("== Gathering information. (Invoke with --help for more details. " - "All arguments are optional)") + print( + '== Gathering information. (Invoke with --help for more details. ' + 'All arguments are optional)' + ) parser = OptionParser() - parser.description = \ - "This script takes an MPAS grid file and generates a SCRIP grid file." + parser.description = ( + 'This script takes an MPAS grid file and generates a SCRIP grid file.' + ) parser.add_option( - "-i", "--input", dest="inputFile", - help="input grid file name used as input.", default="input.nc", - metavar="FILENAME") + '-i', + '--input', + dest='inputFile', + help='input grid file name used as input.', + default='input.nc', + metavar='FILENAME', + ) parser.add_option( - "-s", "--scrip", dest="scripFile", - help="SCRIP grid file to output.", default="scrip.nc", - metavar="FILENAME") + '-s', + '--scrip', + dest='scripFile', + help='SCRIP grid file to output.', + default='scrip.nc', + metavar='FILENAME', + ) parser.add_option( - "-p", "--proj", dest="projection", - help=f"projection used by the input data file. Valid options are: " - f"{projections.keys()}", - metavar="PROJ") + '-p', + '--proj', + dest='projection', + help=f'projection used by the input data file. Valid options are: ' + f'{projections.keys()}', + metavar='PROJ', + ) parser.add_option( - "-r", "--rank", dest="gridRank", - help="desired rank of the output SCRIP grid data") + '-r', + '--rank', + dest='gridRank', + help='desired rank of the output SCRIP grid data', + ) parser.add_option( - "--plot", dest="plot", action="store_true", - help="if this flag is used, destination grid points are plotted") + '--plot', + dest='plot', + action='store_true', + help='if this flag is used, destination grid points are plotted', + ) for option in parser.option_list: - if option.default != ("NO", "DEFAULT"): - option.help += (" " if option.help else "") + "[default: %default]" + if option.default != ('NO', 'DEFAULT'): + option.help += (' ' if option.help else '') + '[default: %default]' options, args = parser.parse_args() if not options.inputFile: - raise ValueError('Data input grid file is required. Specify with -c ' - 'command line argument.') + raise ValueError( + 'Data input grid file is required. Specify with -c ' + 'command line argument.' + ) if not options.scripFile: - raise ValueError('SCRIP output grid file is required. Specify with ' - '-s command line argument.') + raise ValueError( + 'SCRIP output grid file is required. Specify with ' + '-s command line argument.' + ) if not options.projection: - raise ValueError(f'data projection required with -p or --proj command ' - f'line argument. Valid options are: ' - f'{projections.keys()}') + raise ValueError( + f'data projection required with -p or --proj command ' + f'line argument. Valid options are: ' + f'{projections.keys()}' + ) if not options.gridRank: - raise ValueError('desired rank of SCRIP output grid data is required. ' - 'Valid options are 1 (for unstructured grid) or 2') + raise ValueError( + 'desired rank of SCRIP output grid data is required. ' + 'Valid options are 1 (for unstructured grid) or 2' + ) # make a space in stdout before further output print('') @@ -77,47 +106,51 @@ def main(): # Write to output file # Dimensions - fout.createDimension("grid_size", nx * ny) - fout.createDimension("grid_corners", 4) + fout.createDimension('grid_size', nx * ny) + fout.createDimension('grid_corners', 4) if int(options.gridRank) == 1: print('grid rank is 1') - fout.createDimension("grid_rank", 1) + fout.createDimension('grid_rank', 1) elif int(options.gridRank) == 2: print('grid rank is 2') - fout.createDimension("grid_rank", 2) + fout.createDimension('grid_rank', 2) else: - raise ValueError(f'grid rank value is invalid: valid options are ' - f'1 or 2 but {options.gridRank} was given.') + raise ValueError( + f'grid rank value is invalid: valid options are ' + f'1 or 2 but {options.gridRank} was given.' + ) # Variables - grid_center_lat = fout.createVariable('grid_center_lat', 'f8', - ('grid_size',)) + grid_center_lat = fout.createVariable( + 'grid_center_lat', 'f8', ('grid_size',) + ) grid_center_lat.units = 'degrees' - grid_center_lon = fout.createVariable('grid_center_lon', 'f8', - ('grid_size',)) + grid_center_lon = fout.createVariable( + 'grid_center_lon', 'f8', ('grid_size',) + ) grid_center_lon.units = 'degrees' - grid_corner_lat = fout.createVariable('grid_corner_lat', 'f8', - ('grid_size', 'grid_corners')) + grid_corner_lat = fout.createVariable( + 'grid_corner_lat', 'f8', ('grid_size', 'grid_corners') + ) grid_corner_lat.units = 'degrees' - grid_corner_lon = fout.createVariable('grid_corner_lon', 'f8', - ('grid_size', 'grid_corners')) + grid_corner_lon = fout.createVariable( + 'grid_corner_lon', 'f8', ('grid_size', 'grid_corners') + ) grid_corner_lon.units = 'degrees' - grid_imask = fout.createVariable('grid_imask', 'i4', - ('grid_size',)) + grid_imask = fout.createVariable('grid_imask', 'i4', ('grid_size',)) grid_imask.units = 'unitless' - grid_dims = fout.createVariable('grid_dims', 'i4', - ('grid_rank',)) + grid_dims = fout.createVariable('grid_dims', 'i4', ('grid_rank',)) # Create matrices of x,y print('Building matrix version of x, y locations.') xmatrix, ymatrix = np.meshgrid(x, y) # get a copy of x that is on the staggered grid and includes both bounding # edges - xc = np.append(x[:] - dx/2.0, x[-1] + dx / 2.0) + xc = np.append(x[:] - dx / 2.0, x[-1] + dx / 2.0) # get a copy of y that is on the staggered grid and includes both bounding # edges - yc = np.append(y[:] - dy/2.0, y[-1] + dy / 2.0) + yc = np.append(y[:] - dy / 2.0, y[-1] + dy / 2.0) xcmatrix, ycmatrix = np.meshgrid(xc, yc) # Unproject to lat/long for grid centers and grid corners @@ -135,8 +168,9 @@ def main(): t = Transformer.from_crs(crs_in, crs_out) # transform the original grid into the lat-lon grid - grid_center_lon[:], grid_center_lat[:] = t.transform(xmatrix_flat, - ymatrix_flat) + grid_center_lon[:], grid_center_lat[:] = t.transform( + xmatrix_flat, ymatrix_flat + ) # Now fill in the corners in the right locations stag_lon, stag_lat = t.transform(xcmatrix, ycmatrix) @@ -168,12 +202,12 @@ def main(): # set the grid dimension based on the grid rank if int(options.gridRank) == 1: - grid_dims[:] = (nx * ny) + grid_dims[:] = nx * ny elif int(options.gridRank) == 2: grid_dims[:] = [nx, ny] if options.plot: - print("plotting is on") + print('plotting is on') # plot some stuff # plot a single point i = -1 From cecf939c62918414ef4bab87286feb1bdde63fc6 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:57:06 -0500 Subject: [PATCH 08/10] Use CF units for grid_imask in SCRIP files "unitless" is not a valid udunits string; CF uses "1" for dimensionless quantities. Co-Authored-By: Claude Opus 5.5 (1M context) --- conda_package/mpas_tools/scrip/from_mpas.py | 2 +- conda_package/mpas_tools/scrip/from_planar.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/conda_package/mpas_tools/scrip/from_mpas.py b/conda_package/mpas_tools/scrip/from_mpas.py index c1ba0c544..1528d36ee 100755 --- a/conda_package/mpas_tools/scrip/from_mpas.py +++ b/conda_package/mpas_tools/scrip/from_mpas.py @@ -100,7 +100,7 @@ def scrip_from_mpas(mpasFile, scripFile, useLandIceMask=False): grid_area = fout.createVariable('grid_area', 'f8', ('grid_size',)) grid_area.units = 'radian^2' grid_imask = fout.createVariable('grid_imask', 'i4', ('grid_size',)) - grid_imask.units = 'unitless' + grid_imask.units = '1' grid_dims = fout.createVariable('grid_dims', 'i4', ('grid_rank',)) grid_center_lat[:] = latCell[:] diff --git a/conda_package/mpas_tools/scrip/from_planar.py b/conda_package/mpas_tools/scrip/from_planar.py index 825682d87..ce22369c8 100644 --- a/conda_package/mpas_tools/scrip/from_planar.py +++ b/conda_package/mpas_tools/scrip/from_planar.py @@ -139,7 +139,7 @@ def main(): ) grid_corner_lon.units = 'degrees' grid_imask = fout.createVariable('grid_imask', 'i4', ('grid_size',)) - grid_imask.units = 'unitless' + grid_imask.units = '1' grid_dims = fout.createVariable('grid_dims', 'i4', ('grid_rank',)) # Create matrices of x,y From 2b6f400fe836e40ae74ee2c777270940e6677ce8 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:57:18 -0500 Subject: [PATCH 09/10] Add tests for CF metadata on meshes Check cf_conventions() and that planar hex meshes, and meshes from MpasCellCuller.x and MpasMeshConverter.x, have the attributes in mpas_tools.mesh.attrs and a CF Conventions entry, including keeping an existing CF version. Co-Authored-By: Claude Opus 5.5 (1M context) --- conda_package/tests/test_cf_attrs.py | 62 ++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 conda_package/tests/test_cf_attrs.py diff --git a/conda_package/tests/test_cf_attrs.py b/conda_package/tests/test_cf_attrs.py new file mode 100644 index 000000000..88d80700f --- /dev/null +++ b/conda_package/tests/test_cf_attrs.py @@ -0,0 +1,62 @@ +import xarray + +from mpas_tools.mesh.attrs import MESH_VAR_ATTRS, cf_conventions +from mpas_tools.mesh.conversion import convert, cull +from mpas_tools.planar_hex import make_planar_hex_mesh + +from .util import get_test_data_file + + +def _check_mesh_attrs(ds): + """ + Check that every mesh variable in the dataset has the attributes from + ``MESH_VAR_ATTRS`` and that the file follows CF and MPAS conventions + """ + for name in ds.data_vars: + if name not in MESH_VAR_ATTRS: + continue + for key, value in MESH_VAR_ATTRS[name].items(): + assert ds[name].attrs.get(key) == value, (name, key) + if 'units' not in MESH_VAR_ATTRS[name]: + assert 'units' not in ds[name].attrs, name + assert ds.attrs['Conventions'].split()[0].startswith('CF-') + assert 'MPAS' in ds.attrs['Conventions'].split() + + +def test_cf_conventions(): + assert cf_conventions() == 'CF-1.8 MPAS' + assert cf_conventions('') == 'CF-1.8 MPAS' + assert cf_conventions('MPAS') == 'CF-1.8 MPAS' + assert cf_conventions('CF-1.8 MPAS') == 'CF-1.8 MPAS' + assert cf_conventions('MPAS CF-1.10') == 'CF-1.10 MPAS' + assert cf_conventions('CF-1.6, ACDD-1.3') == 'CF-1.6 ACDD-1.3 MPAS' + + +def test_planar_hex_cf_attrs(): + ds = make_planar_hex_mesh( + nx=10, ny=10, dc=1e3, nonperiodic_x=False, nonperiodic_y=True + ) + _check_mesh_attrs(ds) + assert ds.attrs['Conventions'] == 'CF-1.8 MPAS' + + +def test_cull_convert_planar_cf_attrs(): + ds = make_planar_hex_mesh( + nx=10, ny=10, dc=1e3, nonperiodic_x=False, nonperiodic_y=True + ) + ds_culled = cull(ds) + _check_mesh_attrs(ds_culled) + ds_converted = convert(ds_culled) + _check_mesh_attrs(ds_converted) + assert ds_converted.attrs['Conventions'] == 'CF-1.8 MPAS' + + +def test_convert_spherical_cf_attrs(): + ds = xarray.open_dataset(get_test_data_file('mesh.QU.1920km.151026.nc')) + # an existing CF version should be kept + ds.attrs['Conventions'] = 'MPAS CF-1.10' + ds_converted = convert(ds) + _check_mesh_attrs(ds_converted) + assert ds_converted.attrs['Conventions'] == 'CF-1.10 MPAS' + for name in ['cellQuality', 'gridSpacing', 'triangleQuality']: + assert name in ds_converted From 69c236544a5a1de245dd1f3079d594cbf5caf1fe Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 06:57:24 -0500 Subject: [PATCH 10/10] Document CF metadata on MPAS meshes Co-Authored-By: Claude Opus 5.5 (1M context) --- conda_package/docs/api.rst | 8 ++++++++ conda_package/docs/mesh_conversion.rst | 10 ++++++++++ 2 files changed, 18 insertions(+) diff --git a/conda_package/docs/api.rst b/conda_package/docs/api.rst index e9dcac9b5..d32004b91 100644 --- a/conda_package/docs/api.rst +++ b/conda_package/docs/api.rst @@ -85,6 +85,14 @@ Mesh conversion calc_edge_normal_vector calc_vector_east_north +.. currentmodule:: mpas_tools.mesh.attrs + +.. autosummary:: + :toctree: generated/ + + add_mesh_attrs + cf_conventions + .. currentmodule:: mpas_tools.merge_grids .. autosummary:: diff --git a/conda_package/docs/mesh_conversion.rst b/conda_package/docs/mesh_conversion.rst index 89c768373..37d9c1c7c 100644 --- a/conda_package/docs/mesh_conversion.rst +++ b/conda_package/docs/mesh_conversion.rst @@ -79,6 +79,16 @@ Optional global attributes (passed through): If present, the ``file_id`` attribute is preserved as ``parent_id`` in the output mesh, and a new ``file_id`` is generated. +The output mesh follows the `CF conventions `_. +Its ``Conventions`` attribute is ``CF-1.8 MPAS``, keeping any CF version +already in the input's ``Conventions``, and each mesh variable has a +``long_name`` and, where they apply, ``units`` and ``standard_name``. +``MpasCellCuller.x``, :py:func:`mpas_tools.planar_hex.make_planar_hex_mesh` +and :py:func:`mpas_tools.mesh.creation.jigsaw_to_netcdf.jigsaw_to_netcdf` +write the same metadata, which is listed in :py:mod:`mpas_tools.mesh.attrs`. +Use :py:func:`mpas_tools.mesh.attrs.add_mesh_attrs` to add it to mesh +variables that other tools create. + The converter also generates a ``graph.info`` file for graph partitioning tools (e.g., Metis). In Python, this file is only written if the ``graphInfoFileName`` argument is provided.