Skip to content

Implements Grid.compute_face_node_angles() - #1672

Open
Sevans711 wants to merge 11 commits into
mainfrom
face_node_angles
Open

Implements Grid.compute_face_node_angles()#1672
Sevans711 wants to merge 11 commits into
mainfrom
face_node_angles

Conversation

@Sevans711

@Sevans711 Sevans711 commented Aug 12, 2026

Copy link
Copy Markdown
Collaborator

Closes #1566

Overview

Implements new function: Grid.compute_face_node_angles(), which computes the internal angles at each node, for each face in a grid. Intentionally excludes many options initially suggested in #1566, as the discussion there revealed:

  • There is no need to support non-spherical geometry or non-convex faces (see also Clarify in docs: UXarray supports only spherical 2D grids with convex faces #1626).
  • The implementation here is extremely fast (e.g., roughly 3 seconds for 3.75km MPAS grid, even when numba has only 4 threads), and the result's size is nontrivial (size=n_faces*n_max_face_nodes), so it doesn't really make sense to include a caching option right now.

Also adds corresponding tests, checking that the angle computations give correct values in a roughly 30,60,90 degree triangle and a hexagon grid. Also checks nan handling in a grid with both pentagons and hexagons (result should be nan wherever n_max_face_nodes dimension index is larger than n_nodes_per_face). For all three cases, checks that the generalized spherical excess formula holds, by comparing to Grid.face_areas results.

Within scope but possibly unexpected change included this PR: refactored _small_angle_of_2_vectors to improve efficiency when inputs are not numpy arrays, and rewrote docstring accordingly. (See also: #1648.) This was necessary to improve efficiency for compute_face_node_angles() without inlining all of the _small_angle_of_2_vectors logic. This is unlikely to affect efficiency elsewhere in the code, because other regions of the code seem to always construct tiny numpy arrays before calling _small_angle_of_2_vectors. Running ASV benchmarks here to check, just in case (but it would be acceptable if benchmarks show no improvements).

Also related (may help visualize face node angle results, if implemented): #1578

Expected Usage

import uxarray as ux

grid_path = "/path/to/grid.nc"
data_path = "/path/to/data.nc"

uxds = ux.open_dataset(grid_path, data_path)
grid = uxds.uxgrid   # or just do: grid = ux.open_grid(grid_path)

# internal angle [in radians] at each node, for each face, as an xr.DataArray
angles_rad = grid.compute_face_node_angles()

# internal angle [in degrees] at each node, for each face, as an xr.DataArray
angles_deg = grid.compute_face_node_angles(degrees=True)

# internal angle [in radians] at each node, for each face, as a ux.UxDataArray
angles_rad_uxarr = grid.compute_face_node_angles(as_uxarray=True)

PR Checklist

General

  • An issue is created and linked
  • Added appropriate labels (if your uxarray repo permissions allow it)
  • Filled out Overview and Expected Usage (if applicable) sections

Testing & Benchmarking

  • Adequate tests are created if there is new functionality
  • Tests are not too basic (such as simply calling a function and nothing else)
  • Tests cover all major paths in your new functions
  • If this PR could affect performance, ran ASV benchmarks and confirmed they show expected behavior (add a new benchmark if necessary)

Documentation

  • Docstrings have been added to all new functions
  • Docstrings have been updated with any function changes
  • User (public) functions have been added to docs/api.rst
  • Internal (private) function names start with an underscore (_)

AI Disclosure

AI Usage: discussions with Claude and ChatGPT, inline code suggestions from GitHub Copilot. Also (not in this PR, but contributed to decisions, see original issue for details) used Claude to make a numpy-only method for efficiency comparisons, and other people used Claude to iterate and run tests for more efficiency improvements.

  • I take responsibility for all AI-generated content in my PR.
  • I have tested all AI-generated content in my PR.

Sevans711 and others added 9 commits July 15, 2026 13:45
This provides an initial implementation which seems to be working. See #1566

TODO:
- support assume_convex=False
- type-hint possible UxDataArray return type?
- tests
- examples
removes "geometry" (always use spherical)
removes "assume_convex" (always assume convex)

Also, improves type-hinting for compute_face_node_angles().

Also, formats using pre-commit ruff formatting.
use tuples instead of tiny arrays in numba!

Cleans up docstring in grid.py (forgot to change in previous commit).
(it was making a type error in numba when trying to run the method)
compute_face_node_angles() seems to be really fast, so including a cache option in initial implementation isn't worthwhile.
@Sevans711 Sevans711 added new feature New feature or request run-benchmark Run ASV benchmark workflow labels Aug 12, 2026
@github-actions

github-actions Bot commented Aug 12, 2026

Copy link
Copy Markdown

ASV Benchmarking

Benchmark Comparison Results

Benchmarks that have improved:

Change Before [7adbd4c] After [f409fc7] Ratio Benchmark (Parameter)
- 518M 337M 0.65 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
- 636M 336M 0.53 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
- 3.54±0.07s 2.46±0s 0.69 mpas_ocean.ConstructFaceLatLon.time_welzl('120km')
- 222±0.9ms 159±0.7ms 0.72 mpas_ocean.ConstructFaceLatLon.time_welzl('480km')
- 442M 331M 0.75 mpas_ocean.FaceAreas.peakmem_compute_face_areas('480km')
- 466M 331M 0.71 mpas_ocean.Gradient.peakmem_gradient('480km')

Benchmarks that have stayed the same:

Change Before [7adbd4c] After [f409fc7] Ratio Benchmark (Parameter)
203±2ms 211±10ms 1.04 bench_connectivity.Connectivity.time_edge_face('120km')
12.2±0.07ms 12.4±0.5ms 1.02 bench_connectivity.Connectivity.time_edge_face('480km')
200±1ms 199±1ms 0.99 bench_connectivity.Connectivity.time_edge_node('120km')
11.1±0.09ms 11.2±0.2ms 1 bench_connectivity.Connectivity.time_edge_node('480km')
198±1ms 199±0.6ms 1.01 bench_connectivity.Connectivity.time_face_edge('120km')
11.3±0.1ms 11.6±0.4ms 1.02 bench_connectivity.Connectivity.time_face_edge('480km')
892±5ms 893±8ms 1 bench_connectivity.Connectivity.time_face_face('120km')
57.5±0.7ms 57.2±0.7ms 0.99 bench_connectivity.Connectivity.time_face_face('480km')
67.5±1μs 65.2±0.8μs 0.97 bench_connectivity.Connectivity.time_face_node('120km')
63.4±2μs 65.8±1μs 1.04 bench_connectivity.Connectivity.time_face_node('480km')
406±5μs 404±10μs 1 bench_connectivity.Connectivity.time_n_nodes_per_face('120km')
349±20μs 342±3μs 0.98 bench_connectivity.Connectivity.time_n_nodes_per_face('480km')
199±0.6ms 209±7ms 1.05 bench_connectivity.Connectivity.time_node_edge('120km')
11.5±0.1ms 11.4±0.05ms 0.99 bench_connectivity.Connectivity.time_node_edge('480km')
78.3±2ms 80.3±4ms 1.03 bench_connectivity.Connectivity.time_node_face('120km')
5.19±0.03ms 5.13±0.02ms 0.99 bench_connectivity.Connectivity.time_node_face('480km')
336M 335M 1 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
367M 367M 1 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
8.35±0.06ms 8.50±0.06ms 1.02 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
2.68±0.04ms 2.69±0.05ms 1.01 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
10.3±0.02ms 10.2±0.02ms 0.99 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
2.17±0.03ms 2.18±0.01ms 1 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
1.18±0.05μs 1.19±0.04μs 1.01 geometry_kernels.AccucrossKernels.time_accucross
2.71±0.02μs 2.71±0.05μs 1 geometry_kernels.AccucrossKernels.time_accucross_pair
471±20ns 471±10ns 1 geometry_kernels.EFTPrimitives.time_acc_sqrt_re
431±10ns 436±10ns 1.01 geometry_kernels.EFTPrimitives.time_diff_of_products
381±9ns 381±30ns 1 geometry_kernels.EFTPrimitives.time_two_prod
396±20ns 396±5ns 1 geometry_kernels.EFTPrimitives.time_two_sum
1.57±0.03μs 1.58±0.03μs 1.01 geometry_kernels.GCAConstLatIntersection.time_accux_constlat_kernel
1.14±0.03μs 1.18±0.03μs 1.04 geometry_kernels.GCAConstLatIntersection.time_gca_const_lat_intersection
1.94±0.04μs 1.97±0.02μs 1.01 geometry_kernels.GCAConstLatIntersection.time_try_gca_const_lat_intersection
1.69±0.03μs 1.77±0.06μs 1.05 geometry_kernels.GCAGCAIntersection.time_accux_gca_kernel
1.36±0.01μs 1.41±0.05μs 1.03 geometry_kernels.GCAGCAIntersection.time_gca_gca_intersection
2.19±0.03μs 2.25±0.04μs 1.03 geometry_kernels.GCAGCAIntersection.time_try_gca_gca_intersection
53.3±0.2μs 54.2±2μs 1.02 geometry_kernels.OrientPredicates.time_on_minor_arc
1.15±0.04μs 1.10±0.02μs 0.96 geometry_kernels.OrientPredicates.time_orient3d_on_sphere
2.71±0.1ms 2.60±0.01ms 0.96 geometry_samebody.SameBodyConstLat.time_accux_dispatch
1.18±0.01ms 1.17±0.01ms 0.99 geometry_samebody.SameBodyConstLat.time_accux_kernel
1.72±0.01ms 1.71±0.02ms 1 geometry_samebody.SameBodyConstLat.time_fp64_dispatch
146±0.4μs 147±0.6μs 1.01 geometry_samebody.SameBodyConstLat.time_fp64_kernel
32.4±0.09ms 32.2±0.01ms 0.99 geometry_samebody_gcagca.SameBodyGcaGca.time_accux_dispatch
10.2±0.01ms 10.2±0.02ms 1 geometry_samebody_gcagca.SameBodyGcaGca.time_accux_kernel
26.3±0ms 26.4±0.04ms 1 geometry_samebody_gcagca.SameBodyGcaGca.time_fp64_dispatch
4.82±0.02ms 4.87±0.01ms 1.01 geometry_samebody_gcagca.SameBodyGcaGca.time_fp64_kernel
804±4ms 815±4ms 1.01 import.Imports.timeraw_import_uxarray
2.70±0.01ms 2.75±0.03ms 1.02 mpas_ocean.CheckNorm.time_check_norm('120km')
2.23±0.03ms 2.23±0.02ms 1 mpas_ocean.CheckNorm.time_check_norm('480km')
863±9ms 844±4ms 0.98 mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')
54.6±0.7ms 55.2±0.5ms 1.01 mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')
665±8μs 669±9μs 1.01 mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km')
610±10μs 581±10μs 0.95 mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')
5.44±0.02ms 5.46±0.01ms 1 mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('120km')
3.94±0.03ms 3.98±0.02ms 1.01 mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('480km')
18.2±0.03ms 18.2±0.02ms 1 mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')
1.04±0.03ms 1.05±0.01ms 1.01 mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')
10.6±0.04ms 10.6±0.03ms 1 mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')
712±9μs 724±20μs 1.02 mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')
710±9ms 711±1ms 1 mpas_ocean.CrossSections.time_const_lat('120km', 1)
361±0.8ms 360±1ms 1 mpas_ocean.CrossSections.time_const_lat('120km', 2)
186±1ms 184±1ms 0.99 mpas_ocean.CrossSections.time_const_lat('120km', 4)
551±3ms 557±0.9ms 1.01 mpas_ocean.CrossSections.time_const_lat('480km', 1)
279±2ms 280±1ms 1 mpas_ocean.CrossSections.time_const_lat('480km', 2)
144±0.5ms 144±1ms 1 mpas_ocean.CrossSections.time_const_lat('480km', 4)
24.1±0.04ms 24.4±0.1ms 1.01 mpas_ocean.DualMesh.time_dual_mesh_construction('120km')
3.22±0.09ms 3.23±0.06ms 1 mpas_ocean.DualMesh.time_dual_mesh_construction('480km')
350M 350M 1 mpas_ocean.FaceAreas.peakmem_compute_face_areas('120km')
60.5±0.1ms 60.7±0.09ms 1 mpas_ocean.FaceAreas.time_compute_face_areas('120km')
6.92±0.1ms 6.84±0.1ms 0.99 mpas_ocean.FaceAreas.time_compute_face_areas('480km')
938±7ms 937±3ms 1 mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)
53.3±1ms 54.7±2ms 1.03 mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)
81.9±0.3ms 82.8±0.3ms 1.01 mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)
5.36±0.1ms 5.56±0.06ms 1.04 mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)
351M 351M 1 mpas_ocean.Gradient.peakmem_gradient('120km')
173±0.3ms 174±0.7ms 1 mpas_ocean.Gradient.time_gradient('120km')
12.6±0.04ms 12.3±0.07ms 0.98 mpas_ocean.Gradient.time_gradient('480km')
361±8μs 376±10μs 1.04 mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('120km')
197±10μs 199±7μs 1.01 mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('480km')
351M 351M 1 mpas_ocean.Integrate.peakmem_integrate('120km')
330M 332M 1.01 mpas_ocean.Integrate.peakmem_integrate('480km')
543±10μs 519±8μs 0.95 mpas_ocean.Integrate.time_integrate('120km')
453±20μs 463±7μs 1.02 mpas_ocean.Integrate.time_integrate('480km')
180±0.9ms 179±1ms 0.99 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude')
180±2ms 183±2ms 1.02 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include')
178±0.8ms 177±0.8ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')
13.0±0.05ms 13.4±0.06ms 1.03 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude')
13.3±0.06ms 13.5±0.1ms 1.02 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include')
13.5±0.3ms 13.5±0.1ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')
401±10μs 404±20μs 1.01 mpas_ocean.PointInPolygon.time_face_search_lonlat('120km')
387±9μs 394±20μs 1.02 mpas_ocean.PointInPolygon.time_face_search_lonlat('480km')
389±10μs 378±8μs 0.97 mpas_ocean.PointInPolygon.time_face_search_xyz('120km')
359±7μs 357±10μs 1 mpas_ocean.PointInPolygon.time_face_search_xyz('480km')
240±0.4ms 237±0.5ms 0.99 mpas_ocean.RemapDownsample.time_bilinear_remapping
284±0.7ms 286±3ms 1.01 mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping
15.5±0.03ms 15.4±0.03ms 1 mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping
1.42±0s 1.37±0s 0.97 mpas_ocean.RemapUpsample.time_bilinear_remapping
36.0±0.3ms 36.1±0.3ms 1 mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping
12.2±0.2ms 12.2±0.3ms 1 mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping
25.2±0.6ms 25.6±0.1ms 1.02 mpas_ocean.ZonalAverage.time_zonal_average('120km')
5.65±0.02ms 5.77±0.09ms 1.02 mpas_ocean.ZonalAverage.time_zonal_average('480km')
325M 326M 1 quad_hexagon.QuadHexagon.peakmem_open_dataset
324M 325M 1 quad_hexagon.QuadHexagon.peakmem_open_grid
6.98±0.09ms 6.95±0.07ms 1 quad_hexagon.QuadHexagon.time_open_dataset
5.94±0.08ms 5.87±0.05ms 0.99 quad_hexagon.QuadHexagon.time_open_grid

@Sevans711

Copy link
Copy Markdown
Collaborator Author

Looks like this PR is having a nice side effect: the optimization to _small_angle_of_2_vectors did in fact improve the mpas_ocean.ConstructFaceLatLon.time_welzl benchmarks. Following the code reveals that it does in fact get called there, so this is probably a real improvement, not just a rounding error.

(The chain of calls is deep, so I'm putting it here for reference: construct_face_centers(method='welzl') --> _populate_face_centerpoints --> _construct_face_centerpoints --> _smallest_enclosing_circle --> _welzl_recursive --> _circle_from_two_points or _circle_from_three_points --> _small_angle_of_2_vectors)

Ignoring the apparent peakmem improvements (see #1605).

Most importantly, none of the benchmarks show performance degradation, which helps demonstrate that the _small_angle_of_2_vectors are at best an improvement, at worst causing no change.

@erogluorhan erogluorhan left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks great to me! Thanks for all the invaluable discussion throughout this!

Comment thread uxarray/grid/utils.py
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

new feature New feature or request run-benchmark Run ASV benchmark workflow

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Compute face_node_angles from uxarray.Grid?

2 participants