diff --git a/conda_package/docs/api.rst b/conda_package/docs/api.rst index d32004b91..ec1fab7c0 100644 --- a/conda_package/docs/api.rst +++ b/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 diff --git a/conda_package/docs/mesh_conversion.rst b/conda_package/docs/mesh_conversion.rst index 37d9c1c7c..56d4ff47b 100644 --- a/conda_package/docs/mesh_conversion.rst +++ b/conda_package/docs/mesh_conversion.rst @@ -83,6 +83,8 @@ 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``. +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`. diff --git a/conda_package/mpas_tools/mesh/attrs.py b/conda_package/mpas_tools/mesh/attrs.py index f236384c3..bd5fcf8e1 100644 --- a/conda_package/mpas_tools/mesh/attrs.py +++ b/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')) diff --git a/conda_package/tests/test_cf_attrs.py b/conda_package/tests/test_cf_attrs.py index 88d80700f..a4865dfde 100644 --- a/conda_package/tests/test_cf_attrs.py +++ b/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 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 85035c63a..48895f024 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 @@ -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]); 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 a088a1e8c..00cc155d4 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 @@ -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]);