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. 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 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 48238f8b8..1113a27f5 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,17 @@ -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.attrs import MESH_VAR_ATTRS, cf_conventions +from mpas_tools.mesh.creation.util import circumcenter def triangle_to_netcdf(node, ele, output_name): @@ -30,8 +35,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 +45,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 +69,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 +79,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 +126,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,38 +143,51 @@ 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 + for name, var in grid.variables.items(): + var.setncatts(MESH_VAR_ATTRS.get(name, {})) + grid.Conventions = cf_conventions() + grid.sync() grid.close() 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..e560973ee 100755 --- a/conda_package/mpas_tools/planar_hex.py +++ b/conda_package/mpas_tools/planar_hex.py @@ -1,19 +1,32 @@ #!/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): +from mpas_tools.mesh.attrs import add_mesh_attrs + + +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 @@ -72,6 +85,8 @@ def make_planar_hex_mesh(nx, ny, dc, nonperiodic_x, # 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) @@ -85,8 +100,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 +113,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 +144,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 +159,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 +201,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 +210,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 +340,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 +386,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 +396,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 +480,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 +500,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__': diff --git a/conda_package/mpas_tools/scrip/from_mpas.py b/conda_package/mpas_tools/scrip/from_mpas.py index f24236f7c..1528d36ee 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,46 +76,52 @@ 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' 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[:] 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..ce22369c8 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.units = 'unitless' - grid_dims = fout.createVariable('grid_dims', 'i4', - ('grid_rank',)) + grid_imask = fout.createVariable('grid_imask', 'i4', ('grid_size',)) + grid_imask.units = '1' + 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 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 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..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 @@ -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); @@ -832,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]); @@ -902,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]); @@ -972,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]); @@ -1060,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]); @@ -1089,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]); @@ -1127,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]); @@ -1152,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]); @@ -1169,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]); @@ -1195,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; @@ -1312,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]); @@ -1396,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]); } @@ -1447,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]); @@ -1536,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 6f38e7424..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 @@ -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); @@ -2564,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]); @@ -2623,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]); @@ -2682,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]); @@ -2751,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]); @@ -2774,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]); @@ -2796,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]); @@ -2812,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]); @@ -2864,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]); @@ -2890,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]); @@ -2907,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]); @@ -2922,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]); @@ -2970,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]); @@ -2992,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]); @@ -3012,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]); @@ -3041,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]); @@ -3072,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]); @@ -3100,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]); @@ -3136,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]) ; @@ -3167,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]); @@ -3188,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]); @@ -3207,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]); 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__