Only give lat/lon standard names to spherical meshes

conda_package/docs/api.rst

@@ -92,6 +92,7 @@ Mesh conversion
 
    add_mesh_attrs
    cf_conventions
+   get_mesh_var_attrs
 
 .. currentmodule:: mpas_tools.merge_grids
 

conda_package/docs/mesh_conversion.rst

@@ -83,6 +83,8 @@ The output mesh follows the `CF conventions <https://cfconventions.org/>`_.
 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``.
+The latitude and longitude variables only have a ``standard_name`` on
+spherical meshes, since they are not geographic coordinates on planar meshes.
 ``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`.

conda_package/mpas_tools/mesh/attrs.py

@@ -6,6 +6,11 @@
 udunits parses and CF uses (``m``, ``m2``, ``radians``, ``1`` for
 dimensionless).  Index and connectivity variables and masks have no units.
 
+``SPHERE_VAR_ATTRS`` holds the ``standard_name`` attributes of the latitude
+and longitude variables, which are only added on spherical meshes because the
+latitude and longitude of planar meshes are not geographic coordinates (they
+are all zero).
+
 ``MpasMeshConverter.x`` and ``MpasCellCuller.x`` write the same attributes,
 so the two need to be kept in sync.
 """
@@ -37,42 +42,21 @@
     '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',
-    },
+    'latCell': {'long_name': 'latitude of cell centers', 'units': 'radians'},
+    'lonCell': {'long_name': 'longitude of cell centers', 'units': 'radians'},
     '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',
-    },
+    'latEdge': {'long_name': 'latitude of edge midpoints', 'units': 'radians'},
     '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',
-    },
+    'latVertex': {'long_name': 'latitude of vertices', 'units': 'radians'},
+    'lonVertex': {'long_name': 'longitude of vertices', 'units': 'radians'},
     # geometry
     'areaCell': {
         'long_name': 'area of each cell in the primal mesh',
@@ -161,6 +145,40 @@
     },
 }
 
+SPHERE_VAR_ATTRS = {
+    'latCell': {'standard_name': 'latitude'},
+    'lonCell': {'standard_name': 'longitude'},
+    'latEdge': {'standard_name': 'latitude'},
+    'lonEdge': {'standard_name': 'longitude'},
+    'latVertex': {'standard_name': 'latitude'},
+    'lonVertex': {'standard_name': 'longitude'},
+}
+
+
+def get_mesh_var_attrs(on_a_sphere):
+    """
+    Get the CF metadata for the MPAS mesh variables on a spherical or planar
+    mesh
+
+    Parameters
+    ----------
+    on_a_sphere : bool
+        Whether the mesh is spherical, in which case latitude and longitude
+        variables get a ``standard_name``
+
+    Returns
+    -------
+    mesh_var_attrs : dict
+        The attributes of each mesh variable
+    """
+    mesh_var_attrs = {
+        name: dict(attrs) for name, attrs in MESH_VAR_ATTRS.items()
+    }
+    if on_a_sphere:
+        for name, attrs in SPHERE_VAR_ATTRS.items():
+            mesh_var_attrs[name].update(attrs)
+    return mesh_var_attrs
+
 
 def cf_conventions(conventions=None):
     """
@@ -194,7 +212,12 @@ def cf_conventions(conventions=None):
 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.
+    ``Conventions``.  Attributes that a variable already has are kept, except
+    that a ``latitude`` or ``longitude`` ``standard_name`` is removed from
+    the latitude and longitude variables of a planar mesh.
+
+    Latitude and longitude variables only get a ``standard_name`` if the
+    ``on_a_sphere`` attribute of the dataset is ``YES``.
 
     Parameters
     ----------
@@ -206,10 +229,15 @@ def add_mesh_attrs(ds):
     ds : xarray.Dataset
         The same dataset, for convenience
     """
-    for name, attrs in MESH_VAR_ATTRS.items():
+    on_a_sphere = str(ds.attrs.get('on_a_sphere', '')).strip() == 'YES'
+    for name, attrs in get_mesh_var_attrs(on_a_sphere).items():
         if name not in ds:
             continue
         var_attrs = ds[name].attrs
+        if not on_a_sphere and name in SPHERE_VAR_ATTRS:
+            # planar meshes made with older versions may have these
+            if var_attrs.get('standard_name') in ['latitude', 'longitude']:
+                del var_attrs['standard_name']
         for key, value in attrs.items():
             var_attrs.setdefault(key, value)
     ds.attrs['Conventions'] = cf_conventions(ds.attrs.get('Conventions'))

conda_package/tests/test_cf_attrs.py

@@ -1,6 +1,11 @@
 import xarray
 
-from mpas_tools.mesh.attrs import MESH_VAR_ATTRS, cf_conventions
+from mpas_tools.mesh.attrs import (
+    SPHERE_VAR_ATTRS,
+    add_mesh_attrs,
+    cf_conventions,
+    get_mesh_var_attrs,
+)
 from mpas_tools.mesh.conversion import convert, cull
 from mpas_tools.planar_hex import make_planar_hex_mesh
 
@@ -10,15 +15,19 @@
 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
+    ``get_mesh_var_attrs()`` and that the file follows CF and MPAS conventions
     """
+    on_a_sphere = ds.attrs['on_a_sphere'].strip() == 'YES'
+    mesh_var_attrs = get_mesh_var_attrs(on_a_sphere)
     for name in ds.data_vars:
-        if name not in MESH_VAR_ATTRS:
+        if name not in mesh_var_attrs:
             continue
-        for key, value in MESH_VAR_ATTRS[name].items():
+        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]:
+        if 'units' not in mesh_var_attrs[name]:
             assert 'units' not in ds[name].attrs, name
+        if 'standard_name' not in mesh_var_attrs[name]:
+            assert 'standard_name' not in ds[name].attrs, name
     assert ds.attrs['Conventions'].split()[0].startswith('CF-')
     assert 'MPAS' in ds.attrs['Conventions'].split()
 
@@ -40,6 +49,38 @@ def test_planar_hex_cf_attrs():
     assert ds.attrs['Conventions'] == 'CF-1.8 MPAS'
 
 
+def test_get_mesh_var_attrs():
+    for name, attrs in SPHERE_VAR_ATTRS.items():
+        assert 'standard_name' not in get_mesh_var_attrs(False)[name]
+        assert (
+            get_mesh_var_attrs(True)[name]['standard_name']
+            == attrs['standard_name']
+        )
+
+
+def test_add_mesh_attrs_planar_lat_lon():
+    ds = make_planar_hex_mesh(
+        nx=4, ny=4, dc=10e3, nonperiodic_x=True, nonperiodic_y=True
+    )
+    for name in SPHERE_VAR_ATTRS:
+        assert 'standard_name' not in ds[name].attrs, name
+    # planar meshes made with older versions have these standard names
+    for name, attrs in SPHERE_VAR_ATTRS.items():
+        ds[name].attrs.update(attrs)
+    add_mesh_attrs(ds)
+    _check_mesh_attrs(ds)
+
+
+def test_add_mesh_attrs_spherical_lat_lon():
+    ds = xarray.open_dataset(get_test_data_file('mesh.QU.1920km.151026.nc'))
+    for name in SPHERE_VAR_ATTRS:
+        ds[name].attrs = {}
+    add_mesh_attrs(ds)
+    _check_mesh_attrs(ds)
+    assert ds.latCell.attrs['standard_name'] == 'latitude'
+    assert ds.lonCell.attrs['standard_name'] == 'longitude'
+
+
 def test_cull_convert_planar_cf_attrs():
     ds = make_planar_hex_mesh(
         nx=10, ny=10, dc=1e3, nonperiodic_x=False, nonperiodic_y=True

mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp

@@ -843,10 +843,10 @@ int mapAndOutputGridCoordinates( const string inputFilename, const string output
 
     ncutil::def_var(outputFilename, "latCell",
         NC_DOUBLE, "latitude of cell centers",
-        {"nCells"}, "radians", "latitude");
+        {"nCells"}, "radians", spherical ? "latitude" : "");
     ncutil::def_var(outputFilename, "lonCell",
         NC_DOUBLE, "longitude of cell centers",
-        {"nCells"}, "radians", "longitude");
+        {"nCells"}, "radians", spherical ? "longitude" : "");
 
     ncutil::put_var(outputFilename, "latCell", &latNew[0]);
     ncutil::put_var(outputFilename, "lonCell", &lonNew[0]);
@@ -915,10 +915,10 @@ int mapAndOutputGridCoordinates( const string inputFilename, const string output
 
     ncutil::def_var(outputFilename, "latEdge",
         NC_DOUBLE, "latitude of edge midpoints",
-        {"nEdges"}, "radians", "latitude");
+        {"nEdges"}, "radians", spherical ? "latitude" : "");
     ncutil::def_var(outputFilename, "lonEdge",
         NC_DOUBLE, "longitude of edge midpoints",
-        {"nEdges"}, "radians", "longitude");
+        {"nEdges"}, "radians", spherical ? "longitude" : "");
 
     ncutil::put_var(outputFilename, "latEdge", &latNew[0]);
     ncutil::put_var(outputFilename, "lonEdge", &lonNew[0]);
@@ -987,10 +987,10 @@ int mapAndOutputGridCoordinates( const string inputFilename, const string output
 
     ncutil::def_var(outputFilename, "latVertex",
         NC_DOUBLE, "latitude of vertices",
-        {"nVertices"}, "radians", "latitude");
+        {"nVertices"}, "radians", spherical ? "latitude" : "");
     ncutil::def_var(outputFilename, "lonVertex",
         NC_DOUBLE, "longitude of vertices",
-        {"nVertices"}, "radians", "longitude");
+        {"nVertices"}, "radians", spherical ? "longitude" : "");
 
     ncutil::put_var(outputFilename, "latVertex", &latNew[0]);
     ncutil::put_var(outputFilename, "lonVertex", &lonNew[0]);

mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp

@@ -2575,10 +2575,10 @@ int outputGridCoordinates( const string outputFilename) {/*{{{*/
 
     ncutil::def_var(outputFilename, "latCell",
         NC_DOUBLE, "latitude of cell centers",
-        {"nCells"}, "radians", "latitude");
+        {"nCells"}, "radians", spherical ? "latitude" : "");
     ncutil::def_var(outputFilename, "lonCell",
         NC_DOUBLE, "longitude of cell centers",
-        {"nCells"}, "radians", "longitude");
+        {"nCells"}, "radians", spherical ? "longitude" : "");
 
     ncutil::put_var(outputFilename, "latCell", &lat[0]);
     ncutil::put_var(outputFilename, "lonCell", &lon[0]);
@@ -2636,10 +2636,10 @@ int outputGridCoordinates( const string outputFilename) {/*{{{*/
 
     ncutil::def_var(outputFilename, "latEdge",
         NC_DOUBLE, "latitude of edge midpoints",
-        {"nEdges"}, "radians", "latitude");
+        {"nEdges"}, "radians", spherical ? "latitude" : "");
     ncutil::def_var(outputFilename, "lonEdge",
         NC_DOUBLE, "longitude of edge midpoints",
-        {"nEdges"}, "radians", "longitude");
+        {"nEdges"}, "radians", spherical ? "longitude" : "");
 
     ncutil::put_var(outputFilename, "latEdge", &lat[0]);
     ncutil::put_var(outputFilename, "lonEdge", &lon[0]);
@@ -2697,10 +2697,10 @@ int outputGridCoordinates( const string outputFilename) {/*{{{*/
 
     ncutil::def_var(outputFilename, "latVertex",
         NC_DOUBLE, "latitude of vertices",
-        {"nVertices"}, "radians", "latitude");
+        {"nVertices"}, "radians", spherical ? "latitude" : "");
     ncutil::def_var(outputFilename, "lonVertex",
         NC_DOUBLE, "longitude of vertices",
-        {"nVertices"}, "radians", "longitude");
+        {"nVertices"}, "radians", spherical ? "longitude" : "");
 
     ncutil::put_var(outputFilename, "latVertex", &lat[0]);
     ncutil::put_var(outputFilename, "lonVertex", &lon[0]);