From 6e02082aa3ffdd51b7f02519e3f179715bf093b4 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 23 Oct 2026 20:40:38 +0200 Subject: [PATCH] Only give lat/lon standard names to spherical meshes The latitude and longitude of planar meshes are all zero and are not geographic coordinates, so CF-aware tools should not treat them as such. Move the latitude and longitude standard names out of MESH_VAR_ATTRS into SPHERE_VAR_ATTRS, add get_mesh_var_attrs() to combine them, and have add_mesh_attrs() only add them when on_a_sphere is YES. On a planar mesh, add_mesh_attrs() also removes these standard names, which meshes made with 2.2.0 have. MpasMeshConverter.x and MpasCellCuller.x likewise only write them for spherical meshes. Update the CF tests to check that planar meshes have no lat/lon standard names and spherical meshes do. Co-Authored-By: Claude Opus 5.5 (1M context) --- conda_package/docs/api.rst | 1 + conda_package/docs/mesh_conversion.rst | 2 + conda_package/mpas_tools/mesh/attrs.py | 84 ++++++++++++------- conda_package/tests/test_cf_attrs.py | 51 +++++++++-- .../mpas_cell_culler.cpp | 12 +-- .../mpas_mesh_converter.cpp | 12 +-- 6 files changed, 117 insertions(+), 45 deletions(-) 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]);