From 0f545615e972dabe3cd8b2f911bbe40d9c85f22a Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 1 Apr 2026 12:11:24 +0200 Subject: [PATCH 1/5] Fix computation of `angleEdge` in NetCDF-C mesh converter The old spherical `angleEdge` calculation in `MpasMeshConverter.x` inferred the edge angle from the latitude difference between edge vertices and then used a sign correction based on a nearby artificial "north pole" point. This approximation becomes inaccurate and unstable near the poles, where local east/north directions are poorly represented by latitude differences. This merge replaces that logic with a Cartesian geometric formulation: compute the edge normal from the vertex chord and edge position on the sphere, project it onto the local east/north basis at the edge, and compute `angleEdge` with `atan2`. This matches the downstream Polaris fix and gives robust spherical angles across the globe, including polar regions. --- .../mpas_mesh_converter.cpp | 83 +++++++++++-------- 1 file changed, 50 insertions(+), 33 deletions(-) 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 9ad3001a6..6f38e7424 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 @@ -97,6 +97,10 @@ int buildAreas(); int buildEdgesOnEdgeArrays(); int buildAngleEdge(); int buildMeshQualities(); +double clampForUnitInterval(const double value); +void buildLocalEastNorth(const pnt &location, pnt &east, pnt &north); +double sphericalAngleEdge(const pnt &edge_loc, const pnt &vertex_loc1, + const pnt &vertex_loc2); /*}}}*/ /* Output functions {{{*/ @@ -2182,10 +2186,10 @@ int buildAngleEdge(){/*{{{*/ int iEdge; int cell1, cell2; int vertex1, vertex2; - pnt np, x_axis, normal; + pnt x_axis, normal; pnt cell_loc1, cell_loc2; pnt vertex_loc1, vertex_loc2; - double angle, sign; + double angle; #ifdef _DEBUG cout << endl << endl << "Begin function: buildAngleEdge" << endl << endl; @@ -2225,52 +2229,66 @@ int buildAngleEdge(){/*{{{*/ cell_loc2.fixPeriodicity(cell_loc1, xPeriodicFix, yPeriodicFix); normal = cell_loc2 - cell_loc1; - angleEdge.at(iEdge) = acos( x_axis.dot(normal) / (x_axis.magnitude() * normal.magnitude())); + angleEdge.at(iEdge) = acos(clampForUnitInterval( + x_axis.dot(normal) / (x_axis.magnitude() * normal.magnitude()))); if (cell_loc2.y < cell_loc1.y) angleEdge.at(iEdge) = 2.0 * M_PI - angleEdge.at(iEdge); } else { - - np = pntFromLatLon(edges.at(iEdge).getLat()+0.05, edges.at(iEdge).getLon()); - np.normalize(); + angle = sphericalAngleEdge(edges.at(iEdge), vertex_loc1, vertex_loc2); #ifdef _DEBUG - cout << " NP: " << np << endl; + cout << " fangle: " << angle << endl; #endif - angle = (vertex_loc2.getLat() - vertex_loc1.getLat()) / dvEdge.at(iEdge); - angle = max( min(angle, 1.0), -1.0); - angle = acos(angle); + angleEdge.at(iEdge) = angle; + } + } -#ifdef _DEBUG - cout << " angle: " << angle << endl; -#endif + return 0; +}/*}}}*/ - sign = planeAngle(edges.at(iEdge), np, vertex_loc2, edges.at(iEdge)); - if(sign != 0.0){ - sign = sign / fabs(sign); - } else { - sign = 1.0; - } +double clampForUnitInterval(const double value){/*{{{*/ + return max(min(value, 1.0), -1.0); +}/*}}}*/ +void buildLocalEastNorth(const pnt &location, pnt &east, pnt &north){/*{{{*/ + /* + * buildLocalEastNorth constructs local eastward and northward unit vectors + * at a point on the unit sphere. + */ + const double lat = location.getLat(); + const double lon = location.getLon(); -#ifdef _DEBUG - cout << " sign : " << sign << endl; - cout << " a*s : " << angle * sign << endl; -#endif + east = pnt(-sin(lon), cos(lon), 0.0); + north = pnt(-sin(lat) * cos(lon), -sin(lat) * sin(lon), cos(lat)); +}/*}}}*/ - angle = angle * sign; - if(angle > M_PI) angle = angle - 2.0 * M_PI; - if(angle < -M_PI) angle = angle + 2.0 * M_PI; +double sphericalAngleEdge(const pnt &edge_loc, const pnt &vertex_loc1, + const pnt &vertex_loc2){/*{{{*/ + /* + * Compute the angle made by the positive normal direction relative to the + * local eastward direction using Cartesian geometry on the sphere. + */ + pnt edge_unit, dvertex, normal, east, north; + double east_component, north_component; -#ifdef _DEBUG - cout << " fangle: " << angle << endl; -#endif + edge_unit = edge_loc; + edge_unit.normalize(); - angleEdge.at(iEdge) = angle; - } + dvertex = vertex_loc2 - vertex_loc1; + normal = dvertex.cross(edge_unit); + + buildLocalEastNorth(edge_unit, east, north); + + east_component = east.dot(normal); + north_component = north.dot(normal); + + if(fabs(east_component) < DBL_EPSILON && fabs(north_component) < DBL_EPSILON){ + return 0.0; } - return 0; + return atan2(north_component, east_component); }/*}}}*/ + int buildMeshQualities(){/*{{{*/ /* * buildMeshQualities constructs fields describing the quality of the mesh, including: @@ -3270,4 +3288,3 @@ string gen_random(const int len) {/*{{{*/ return rand_str; }/*}}}*/ - From 02bf3ff77f501a1111c356336a15c9d091916b1e Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 1 Apr 2026 12:17:13 +0200 Subject: [PATCH 2/5] Add Python helpers to recompute spherical angleEdge Add `mpas_tools.mesh.spherical` with utilities to compute local east/north vectors, edge-normal vectors, and recomputed spherical `angleEdge` values from mesh geometry. These routines mirror the corrected spherical formulation used for the C++ mesh converter and provide a reusable Python reference implementation. --- conda_package/docs/api.rst | 9 ++ conda_package/docs/mesh_conversion.rst | 5 + conda_package/mpas_tools/mesh/spherical.py | 105 +++++++++++++++++++++ 3 files changed, 119 insertions(+) create mode 100644 conda_package/mpas_tools/mesh/spherical.py diff --git a/conda_package/docs/api.rst b/conda_package/docs/api.rst index bc3c95c48..0cd784d01 100644 --- a/conda_package/docs/api.rst +++ b/conda_package/docs/api.rst @@ -76,6 +76,15 @@ Mesh conversion compute_mpas_flood_fill_mask compute_lon_lat_region_masks +.. currentmodule:: mpas_tools.mesh.spherical + +.. autosummary:: + :toctree: generated/ + + recompute_angle_edge + calc_edge_normal_vector + calc_vector_east_north + .. currentmodule:: mpas_tools.merge_grids .. autosummary:: diff --git a/conda_package/docs/mesh_conversion.rst b/conda_package/docs/mesh_conversion.rst index 16f62b4b4..89c768373 100644 --- a/conda_package/docs/mesh_conversion.rst +++ b/conda_package/docs/mesh_conversion.rst @@ -83,6 +83,11 @@ 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. +For spherical meshes, the :py:mod:`mpas_tools.mesh.spherical` module provides +Python utilities for recomputing ``angleEdge`` and related local east/north +geometry directly from the mesh coordinates. This can be useful for +verification and diagnostics after conversion. + .. _cell_culler: Cell Culler diff --git a/conda_package/mpas_tools/mesh/spherical.py b/conda_package/mpas_tools/mesh/spherical.py new file mode 100644 index 000000000..3891f684d --- /dev/null +++ b/conda_package/mpas_tools/mesh/spherical.py @@ -0,0 +1,105 @@ +import numpy as np +import xarray as xr + +from mpas_tools.transects import lon_lat_to_cartesian + + +def recompute_angle_edge(ds_mesh): + """ + Recompute ``angleEdge`` from edge and vertex locations on the sphere. + + Parameters + ---------- + ds_mesh : xarray.Dataset + An MPAS spherical mesh dataset containing edge and vertex locations. + + Returns + ------- + angle_edge : xarray.DataArray + ``angleEdge`` recomputed from spherical geometry. + """ + normal_east_north = calc_edge_normal_vector(ds_mesh) + angle_edge = xr.zeros_like(ds_mesh.angleEdge) + angle_edge.values = np.atan2( + normal_east_north[:, 1], normal_east_north[:, 0] + ) + return angle_edge + + +def calc_edge_normal_vector(ds_mesh): + """ + Compute edge-normal vectors projected onto local east/north coordinates. + + Parameters + ---------- + ds_mesh : xarray.Dataset + An MPAS spherical mesh dataset containing edge and vertex locations. + + Returns + ------- + normal_east_north : numpy.ndarray + A ``(nEdges, 2)`` array of unit normal vectors in local east/north + coordinates. + """ + edge_cartesian = np.array( + lon_lat_to_cartesian( + ds_mesh.lonEdge, ds_mesh.latEdge, 1.0, degrees=False + ) + ) + + vertex_1 = ds_mesh.verticesOnEdge.isel(TWO=0).values - 1 + vertex_2 = ds_mesh.verticesOnEdge.isel(TWO=1).values - 1 + + lon_vertex_1 = ds_mesh.lonVertex.isel(nVertices=vertex_1) + lat_vertex_1 = ds_mesh.latVertex.isel(nVertices=vertex_1) + lon_vertex_2 = ds_mesh.lonVertex.isel(nVertices=vertex_2) + lat_vertex_2 = ds_mesh.latVertex.isel(nVertices=vertex_2) + + vertex_1_cartesian = np.array( + lon_lat_to_cartesian(lon_vertex_1, lat_vertex_1, 1.0, degrees=False) + ) + vertex_2_cartesian = np.array( + lon_lat_to_cartesian(lon_vertex_2, lat_vertex_2, 1.0, degrees=False) + ) + + dvertex_cartesian = vertex_2_cartesian - vertex_1_cartesian + normal_cartesian = np.cross(dvertex_cartesian, edge_cartesian, axis=0) + + edge_east, edge_north = calc_vector_east_north( + edge_cartesian[0, :], edge_cartesian[1, :], edge_cartesian[2, :] + ) + + normal_east_north = np.zeros((ds_mesh.sizes['nEdges'], 2)) + normal_east_north[:, 0] = np.sum(edge_east * normal_cartesian, axis=0) + normal_east_north[:, 1] = np.sum(edge_north * normal_cartesian, axis=0) + + norm = np.linalg.norm(normal_east_north, axis=1) + nonzero = norm > 0.0 + normal_east_north[nonzero, :] /= norm[nonzero, np.newaxis] + + return normal_east_north + + +def calc_vector_east_north(x, y, z): + """ + Compute local east and north unit vectors on the sphere. + + Parameters + ---------- + x, y, z : numpy.ndarray + Cartesian coordinates of points on the unit sphere. + + Returns + ------- + east, north : tuple of numpy.ndarray + Local east and north unit vectors, each with shape ``(3, nPoints)``. + """ + axis = np.array([0.0, 0.0, 1.0]) + xyz = np.stack((x, y, z), axis=1) + east = np.cross(axis, np.transpose(xyz), axis=0) + north = np.cross(np.transpose(xyz), east, axis=0) + + east /= np.linalg.norm(east, axis=0) + north /= np.linalg.norm(north, axis=0) + + return east, north From fcc8c325e51712a082debed3fe4560f512e65853 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 1 Apr 2026 12:18:06 +0200 Subject: [PATCH 3/5] Add regression test for spherical angleEdge in mesh conversion Add a conversion test that runs `MpasMeshConverter.x` on the low-resolution QU test mesh, recomputes `angleEdge` in Python from spherical geometry, and verifies the converted mesh agrees to tight tolerance. This gives us a lightweight regression test for the polar `angleEdge` fix without introducing a dedicated C++ unit-test harness. --- conda_package/tests/test_conversion.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/conda_package/tests/test_conversion.py b/conda_package/tests/test_conversion.py index 348da88cf..c972c9007 100755 --- a/conda_package/tests/test_conversion.py +++ b/conda_package/tests/test_conversion.py @@ -1,9 +1,11 @@ #!/usr/bin/env python import matplotlib +import numpy as np from mpas_tools.io import write_netcdf from mpas_tools.mesh.conversion import convert, cull, mask +from mpas_tools.mesh.spherical import recompute_angle_edge from .util import get_test_data_file @@ -30,5 +32,21 @@ def test_conversion(): write_netcdf(dsMask, 'antarctic_mask.nc') +def test_conversion_angle_edge(): + ds_mesh = xarray.open_dataset( + get_test_data_file('mesh.QU.1920km.151026.nc') + ) + ds_mesh = convert(dsIn=ds_mesh) + + angle_edge_python = recompute_angle_edge(ds_mesh) + angle_diff = np.angle( + np.exp(1j * (angle_edge_python.values - ds_mesh.angleEdge.values)) + ) + + assert np.all(np.isfinite(angle_diff)) + assert np.max(np.abs(angle_diff)) < 1.0e-10 + + if __name__ == '__main__': test_conversion() + test_conversion_angle_edge() From 9d2d7390065e84fd32d8ac33e3a34de433604a6e Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 3 Apr 2026 12:37:50 +0200 Subject: [PATCH 4/5] Fix URL formatting in the docs --- conda_package/docs/authors.rst | 4 ++-- conda_package/docs/building_docs.rst | 2 +- conda_package/docs/making_changes.rst | 12 ++++++------ conda_package/docs/releasing.rst | 18 +++++++++--------- conda_package/docs/visualization.rst | 2 +- 5 files changed, 19 insertions(+), 19 deletions(-) diff --git a/conda_package/docs/authors.rst b/conda_package/docs/authors.rst index 13e64be2a..3bb818ec6 100644 --- a/conda_package/docs/authors.rst +++ b/conda_package/docs/authors.rst @@ -23,5 +23,5 @@ Contributors * Phillip J. Wolfram * Tong Zhang -For a list of all the contributions: -https://github.com/MPAS-Dev/MPAS-Tools/graphs/contributors +For a list of all contributions, see the +`contributors graph `_. diff --git a/conda_package/docs/building_docs.rst b/conda_package/docs/building_docs.rst index 329e62e93..30894e25b 100644 --- a/conda_package/docs/building_docs.rst +++ b/conda_package/docs/building_docs.rst @@ -26,4 +26,4 @@ To preview the documentation locally, open the ``index.html`` file in the cd _build/html python -m http.server 8000 -Then, open http://0.0.0.0:8000/master/ in your browser. +Then, open ``_ in your browser. diff --git a/conda_package/docs/making_changes.rst b/conda_package/docs/making_changes.rst index 6b7c26e6e..3c36ebf8d 100644 --- a/conda_package/docs/making_changes.rst +++ b/conda_package/docs/making_changes.rst @@ -21,9 +21,9 @@ things like whitespace at the end of lines. The first time you set up the ``mpas_tools_dev`` environment, you will need to set up ``pre-commit``. This is done by running: -```bash -pre-commit install -``` +.. code-block:: bash + + pre-commit install You only need to do this once when you create the ``mpas_tools_dev`` environment. If you create a new version of ``mpas_tools_dev``, then you will @@ -43,9 +43,9 @@ PEP8 compliance, as well as sort, check and format imports, f-strings, and `mypy ` to check for consistent variable types. An example error might be: -```bash -example.py:77:1: E302 expected 2 blank lines, found 1 -``` +.. code-block:: bash + + example.py:77:1: E302 expected 2 blank lines, found 1 For this example, we would just add an additional blank line after line 77 and try the commit again to make sure we've resolved the issue. diff --git a/conda_package/docs/releasing.rst b/conda_package/docs/releasing.rst index 4366187cc..80a5be486 100644 --- a/conda_package/docs/releasing.rst +++ b/conda_package/docs/releasing.rst @@ -17,8 +17,8 @@ Version Bump and Dependency Updates - ``conda_package/mpas_tools/__init__.py`` - ``conda_package/recipe/recipe.yaml`` - - Make sure the version follows semantic versioning (see - https://semver.org/). + - Make sure the version follows + `semantic versioning `_. For release candidates, use versions like ``1.3.0rc1`` (no ``v`` prefix). 2. **Check and Update Dependencies** @@ -38,7 +38,7 @@ Version Bump and Dependency Updates - The dependencies in ``recipe.yaml`` are the ones that will be used for the released package on conda-forge. The dependencies in ``pyproject.toml`` are for PyPI and should be kept in sync as much as possible but are only - there as a sanity check when we run ```pip check``. The ``dev-spec.txt`` + there as a sanity check when we run ``pip check``. The ``dev-spec.txt`` file should include all dependencies needed for development and testing, and ``pixi.toml`` should remain equivalent for pixi users. @@ -100,8 +100,7 @@ Tagging and Publishing a Release Candidate - Update dependencies if needed - Commit, push to a new branch, and open a PR **against the ``dev`` branch** - of the feedstock: - https://github.com/conda-forge/mpas_tools-feedstock + of the `mpas_tools-feedstock `_. - Follow any instructions in the PR template and merge once approved @@ -112,7 +111,7 @@ Publishing a Stable Release - For stable releases, create a GitHub release page as follows: - - Go to https://github.com/MPAS-Dev/MPAS-Tools/releases + - Go to `the GitHub releases page `_ - Click "Draft a new release" @@ -127,8 +126,8 @@ Publishing a Stable Release 7. **Updating the conda-forge Feedstock for a Stable Release** - - Wait for the ``regro-cf-autotick-bot`` to open a PR at: - https://github.com/conda-forge/mpas_tools-feedstock + - Wait for the ``regro-cf-autotick-bot`` to open a PR at the + `mpas_tools-feedstock repository `_. - This may take several hours to a day. @@ -137,7 +136,8 @@ Publishing a Stable Release - Merge once CI checks pass **Note:** If you are impatient, you can accelerate this process by creating - a bot issue at: https://github.com/conda-forge/mpas_tools-feedstock/issues + a bot issue at the + `mpas_tools-feedstock issues page `_ with the subject ``@conda-forge-admin, please update version``. This will open a new PR with the version within a few minutes. diff --git a/conda_package/docs/visualization.rst b/conda_package/docs/visualization.rst index 4c42b762c..5c423a899 100644 --- a/conda_package/docs/visualization.rst +++ b/conda_package/docs/visualization.rst @@ -49,7 +49,7 @@ centers to triangle nodes. ``dsTris`` includes variables ``triCellIndices``, the cell that each triangle is part of; ``nodeCellIndices`` and ``nodeCellWeights``, the indices and weights used to interpolate from MPAS cell centers to triangle nodes; Cartesian coordinates ``xNode``, ``yNode``, and -``zNode``; and ``lonNode``` and ``latNode`` in radians. ``lonNode`` is +``zNode``; and ``lonNode`` and ``latNode`` in radians. ``lonNode`` is guaranteed to be within 180 degrees of the cell center corresponding to ``triCellIndices``. Nodes always have a counterclockwise winding. From 5f8f2725cde2ce3149ee52ce26cc2417e7ca8647 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 3 Apr 2026 13:46:36 +0200 Subject: [PATCH 5/5] Add a note about bumping the build number --- conda_package/docs/testing_changes.rst | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/conda_package/docs/testing_changes.rst b/conda_package/docs/testing_changes.rst index f4d706a34..fb05ee4ea 100644 --- a/conda_package/docs/testing_changes.rst +++ b/conda_package/docs/testing_changes.rst @@ -25,9 +25,9 @@ command: cd conda_package pixi install pixi shell - rattler-build build -m ci/linux_64_python3.14.____cpython.yaml -r recipe/ --output-dir ../output + rattler-build build -m ci/linux_64_python3.14.____cpython.yaml -r recipe/ --output-dir output -This writes package artifacts to ``output/`` in the repository root. +This writes package artifacts to ``output/`` under ``conda_package``. To install the locally built package into the pixi environment, add the local build output as a channel and then add ``mpas_tools`` from that channel: @@ -35,8 +35,20 @@ build output as a channel and then add ``mpas_tools`` from that channel: .. code-block:: bash cd conda_package - pixi workspace channel add "file://$PWD/../output" - pixi add --platform linux-64 "mpas_tools [channel='file://$PWD/../output']" + pixi workspace channel add "file://$PWD/output" + pixi add --platform linux-64 "mpas_tools [channel='file://$PWD/output']" + +.. important:: + + pixi, like other conda-family package managers, identifies a package by + its name, version and build string. If you rebuild a local package without + changing that identity, pixi may continue using an older cached artifact + even if the file in ``output/`` has changed. + + If you rebuild ``mpas_tools`` locally and need pixi to pick up the new + package contents reliably, bump the conda recipe build number in + ``conda_package/recipe/recipe.yaml`` before rebuilding. For Python-only + development, ``pixi run install-editable`` is often more convenient. .. warning:: @@ -75,7 +87,8 @@ Then run tools within the pixi shell (for example ``pytest``). A useful hybrid workflow is to install the latest release conda package first (to get compiled tools), then install your branch in editable mode on - top for Python development. + top for Python development. This also avoids the need to bump the conda + build number for every local Python-only rebuild. Legacy Method: Conda Editable Install *************************************