Skip to content

Implements Grid.compute_face_node_angles() - #1672

Merged
Sevans711 merged 15 commits into
mainfrom
face_node_angles
Aug 18, 2026
Merged

Implements Grid.compute_face_node_angles()#1672
Sevans711 merged 15 commits into
mainfrom
face_node_angles

Conversation

@Sevans711

@Sevans711Sevans711 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

importuxarrayasuxgrid_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.DataArrayangles_rad=grid.compute_face_node_angles()
# internal angle [in degrees] at each node, for each face, as an xr.DataArrayangles_deg=grid.compute_face_node_angles(degrees=True)
# internal angle [in radians] at each node, for each face, as a ux.UxDataArrayangles_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.

Sevans711and 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.
@Sevans711Sevans711 added new feature New feature or request run-benchmark Run ASV benchmark workflow labels Aug 12, 2026
@github-actions

github-actionsBot commented Aug 12, 2026

Copy link
Copy Markdown

ASV Benchmarking

Benchmark Comparison Results

Benchmarks that have improved:

ChangeBefore [fda71bb]After [de36c56]RatioBenchmark (Parameter)
-3.45±0.01s2.49±0.01s0.72mpas_ocean.ConstructFaceLatLon.time_welzl('120km')
-224±0.7ms162±1ms0.72mpas_ocean.ConstructFaceLatLon.time_welzl('480km')

Benchmarks that have stayed the same:

ChangeBefore [fda71bb]After [de36c56]RatioBenchmark (Parameter)
200±0.2ms204±1ms1.02bench_connectivity.Connectivity.time_edge_face('120km')
12.2±0.1ms12.3±0.08ms1.01bench_connectivity.Connectivity.time_edge_face('480km')
200±1ms203±0.4ms1.01bench_connectivity.Connectivity.time_edge_node('120km')
11.3±0.03ms11.3±0.06ms1.00bench_connectivity.Connectivity.time_edge_node('480km')
199±1ms203±0.6ms1.02bench_connectivity.Connectivity.time_face_edge('120km')
11.4±0.1ms12.1±0.4ms1.06bench_connectivity.Connectivity.time_face_edge('480km')
898±6ms921±5ms1.03bench_connectivity.Connectivity.time_face_face('120km')
57.7±0.6ms58.0±4ms1.01bench_connectivity.Connectivity.time_face_face('480km')
70.6±6μs69.2±3μs0.98bench_connectivity.Connectivity.time_face_node('120km')
68.0±2μs67.5±1μs0.99bench_connectivity.Connectivity.time_face_node('480km')
426±5μs431±8μs1.01bench_connectivity.Connectivity.time_n_nodes_per_face('120km')
355±10μs360±8μs1.01bench_connectivity.Connectivity.time_n_nodes_per_face('480km')
200±0.6ms204±0.6ms1.02bench_connectivity.Connectivity.time_node_edge('120km')
11.6±0.09ms11.6±0.08ms1.00bench_connectivity.Connectivity.time_node_edge('480km')
82.4±0.9ms80.2±0.6ms0.97bench_connectivity.Connectivity.time_node_face('120km')
5.35±0.07ms5.29±0.06ms0.99bench_connectivity.Connectivity.time_node_face('480km')
8.52±0.09ms8.51±0.07ms1.00face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
2.76±0.1ms2.79±0.04ms1.01face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
10.5±10s10.3±10ms~0.00face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
2.26±0.01ms2.18±0.04ms0.97face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
57.3k57.3k1.00face_bounds.FaceBounds.track_nbytes_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
12.3k12.3k1.00face_bounds.FaceBounds.track_nbytes_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
123k123k1.00face_bounds.FaceBounds.track_nbytes_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
1281281.00face_bounds.FaceBounds.track_nbytes_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
1.27M1.27M1.00face_bounds.FaceBounds.track_nbytes_grid_with_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
50.1k50.1k1.00face_bounds.FaceBounds.track_nbytes_grid_with_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
1.48M1.48M1.00face_bounds.FaceBounds.track_nbytes_grid_with_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
7127121.00face_bounds.FaceBounds.track_nbytes_grid_with_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
1.98M1.98M1.00face_bounds.FaceBounds.track_peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
1.99M1.99M1.00face_bounds.FaceBounds.track_peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
2.15M2.15M1.00face_bounds.FaceBounds.track_peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
38.3k38.3k1.00face_bounds.FaceBounds.track_peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
336M336M1.00face_bounds.FaceBoundsColdStartRss.peakmem_open_and_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
365M367M1.01face_bounds.FaceBoundsColdStartRss.peakmem_open_and_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
338M338M1.00face_bounds.FaceBoundsColdStartRss.peakmem_open_and_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
337M337M1.00face_bounds.FaceBoundsColdStartRss.peakmem_open_and_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
1.20±0.06μs1.22±0.03μs1.02geometry_kernels.AccucrossKernels.time_accucross
2.79±0.04μs2.75±0.04μs0.99geometry_kernels.AccucrossKernels.time_accucross_pair
471±20ns451±10ns0.96geometry_kernels.EFTPrimitives.time_acc_sqrt_re
446±20ns441±30ns0.99geometry_kernels.EFTPrimitives.time_diff_of_products
390±20ns426±20ns1.09geometry_kernels.EFTPrimitives.time_two_prod
411±30ns400±20ns0.97geometry_kernels.EFTPrimitives.time_two_sum
1.68±0.06μs1.59±0.06μs0.95geometry_kernels.GCAConstLatIntersection.time_accux_constlat_kernel
1.13±0.02μs1.13±0.01μs1.00geometry_kernels.GCAConstLatIntersection.time_gca_const_lat_intersection
1.92±0.03μs1.94±0.02μs1.01geometry_kernels.GCAConstLatIntersection.time_try_gca_const_lat_intersection
1.72±0.02μs1.68±0.05μs0.98geometry_kernels.GCAGCAIntersection.time_accux_gca_kernel
1.44±0.04μs1.41±0.03μs0.98geometry_kernels.GCAGCAIntersection.time_gca_gca_intersection
2.19±0.02μs2.22±0.02μs1.02geometry_kernels.GCAGCAIntersection.time_try_gca_gca_intersection
52.9±0.8μs53.1±1μs1.00geometry_kernels.OrientPredicates.time_on_minor_arc
1.09±0.03μs1.12±0.04μs1.03geometry_kernels.OrientPredicates.time_orient3d_on_sphere
2.72±0.1ms2.63±0.02ms0.97geometry_samebody.SameBodyConstLat.time_accux_dispatch
1.17±0.01ms1.17±0.01ms1.00geometry_samebody.SameBodyConstLat.time_accux_kernel
1.73±0.01ms1.72±0.01ms1.00geometry_samebody.SameBodyConstLat.time_fp64_dispatch
147±0.7μs147±0.3μs1.00geometry_samebody.SameBodyConstLat.time_fp64_kernel
32.3±0ms32.5±0.2ms1.01geometry_samebody_gcagca.SameBodyGcaGca.time_accux_dispatch
10.3±0ms11.9±0.8ms~1.16geometry_samebody_gcagca.SameBodyGcaGca.time_accux_kernel
27.0±0.2ms26.7±0.1ms0.99geometry_samebody_gcagca.SameBodyGcaGca.time_fp64_dispatch
4.90±0ms4.96±0.07ms1.01geometry_samebody_gcagca.SameBodyGcaGca.time_fp64_kernel
833±3ms830±2ms1.00import.Imports.timeraw_import_uxarray
293M293M1.00import.Imports.track_peakmem_import_uxarray
2.76±0.02ms2.85±0.04ms1.03mpas_ocean.CheckNorm.time_check_norm('120km')
2.27±0.03ms2.31±0.01ms1.02mpas_ocean.CheckNorm.time_check_norm('480km')
858±10ms851±9ms0.99mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')
55.2±0.3ms56.0±0.5ms1.01mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')
678±10μs677±9μs1.00mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km')
587±10μs606±10μs1.03mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')
5.46±0.03ms5.47±0.04ms1.00mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('120km')
3.99±0.02ms3.96±0.04ms0.99mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('480km')
18.2±0.02ms18.2±0.01ms1.00mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')
1.05±0.01ms1.04±0.02ms0.99mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')
10.6±0.03ms10.6±0.03ms1.00mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')
730±7μs732±10μs1.00mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')
587±4ms591±2ms1.01mpas_ocean.CrossSections.time_const_lat('120km', 1)
297±0.7ms296±1ms1.00mpas_ocean.CrossSections.time_const_lat('120km', 2)
155±0.5ms156±1ms1.01mpas_ocean.CrossSections.time_const_lat('120km', 4)
534±0.7ms536±3ms1.00mpas_ocean.CrossSections.time_const_lat('480km', 1)
270±0.8ms270±1ms1.00mpas_ocean.CrossSections.time_const_lat('480km', 2)
139±0.3ms139±0.6ms1.00mpas_ocean.CrossSections.time_const_lat('480km', 4)
356M356M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('120km', 1)
356M356M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('120km', 2)
355M356M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('120km', 4)
341M339M0.99mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('480km', 1)
339M339M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('480km', 2)
339M339M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('480km', 4)
25.0±0.4ms24.5±0.07ms0.98mpas_ocean.DualMesh.time_dual_mesh_construction('120km')
3.21±0.07ms3.25±0.06ms1.01mpas_ocean.DualMesh.time_dual_mesh_construction('480km')
61.6±1ms61.7±0.6ms1.00mpas_ocean.FaceAreas.time_face_areas('120km')
4.90±5s7.88±5ms~0.00mpas_ocean.FaceAreas.time_face_areas('480km')
229k229k1.00mpas_ocean.FaceAreas.track_nbytes_face_areas('120km')
14.3k14.3k1.00mpas_ocean.FaceAreas.track_nbytes_face_areas('480km')
2.12M2.12M1.00mpas_ocean.FaceAreas.track_peakmem_face_areas('120km')
823k823k1.00mpas_ocean.FaceAreas.track_peakmem_face_areas('480km')
949±9ms943±10ms0.99mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)
54.5±0.8ms55.2±2ms1.01mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)
84.5±0.3ms84.6±0.5ms1.00mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)
5.76±0.05ms5.98±0.1ms1.04mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)
175±0.7ms175±0.2ms1.00mpas_ocean.Gradient.time_gradient('120km')
12.2±0.01ms12.2±0.06ms1.00mpas_ocean.Gradient.time_gradient('480km')
457k457k1.00mpas_ocean.Gradient.track_nbytes_gradient('120km')
28.7k28.7k1.00mpas_ocean.Gradient.track_nbytes_gradient('480km')
5.08M5.08M1.00mpas_ocean.Gradient.track_peakmem_gradient('120km')
328k328k1.00mpas_ocean.Gradient.track_peakmem_gradient('480km')
351M352M1.00mpas_ocean.GradientColdStartRss.peakmem_gradient('120km')
331M331M1.00mpas_ocean.GradientColdStartRss.peakmem_gradient('480km')
371±10μs370±10μs1.00mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('120km')
201±10μs199±3μs0.99mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('480km')
566±9μs579±6μs1.02mpas_ocean.Integrate.time_integrate('120km')
494±5μs500±10μs1.01mpas_ocean.Integrate.time_integrate('480km')
18.4M18.4M1.00mpas_ocean.Integrate.track_nbytes_integrate('120km')
1.2M1.2M1.00mpas_ocean.Integrate.track_nbytes_integrate('480km')
184±1ms183±1ms1.00mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude')
187±6ms184±1ms0.98mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include')
183±2ms182±1ms1.00mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')
13.7±0.09ms13.7±0.08ms1.00mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude')
14.0±0.08ms13.7±0.03ms0.98mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include')
13.6±0.09ms13.7±0.2ms1.01mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')
430±20μs433±10μs1.01mpas_ocean.PointInPolygon.time_face_search_lonlat('120km')
414±6μs438±10μs1.06mpas_ocean.PointInPolygon.time_face_search_lonlat('480km')
407±9μs402±10μs0.99mpas_ocean.PointInPolygon.time_face_search_xyz('120km')
408±10μs383±7μs0.94mpas_ocean.PointInPolygon.time_face_search_xyz('480km')
244±1ms239±0.7ms0.98mpas_ocean.RemapDownsample.time_bilinear_remapping
290±1ms281±0.8ms0.97mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping
15.6±0.1ms15.7±0.1ms1.00mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping
1.44±0.01s1.39±0s0.97mpas_ocean.RemapUpsample.time_bilinear_remapping
36.0±0.3ms37.0±0.4ms1.03mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping
12.5±0.2ms12.3±0.1ms0.99mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping
9.29±0.2ms9.12±0.04ms0.98mpas_ocean.ZonalAverage.time_zonal_average('120km')
4.82±0.06ms4.79±0.04ms1.00mpas_ocean.ZonalAverage.time_zonal_average('480km')
358M358M1.00mpas_ocean.ZonalAveragePeakMem.peakmem_zonal_average('120km')
341M341M1.00mpas_ocean.ZonalAveragePeakMem.peakmem_zonal_average('480km')
6.83±0.03ms6.84±0.01ms1.00quad_hexagon.QuadHexagon.time_open_dataset
5.86±0.01ms5.87±0.01ms1.00quad_hexagon.QuadHexagon.time_open_grid
4084081.00quad_hexagon.QuadHexagon.track_nbytes_open_dataset
3923921.00quad_hexagon.QuadHexagon.track_nbytes_open_grid
73.5k73.8k1.00quad_hexagon.QuadHexagon.track_peakmem_open_dataset
72.7k73k1.00quad_hexagon.QuadHexagon.track_peakmem_open_grid

@Sevans711

Copy link
Copy Markdown
CollaboratorAuthor

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.

@erogluorhanerogluorhan 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 threaduxarray/grid/utils.py
Comment threadtest/grid/geometry/test_angles.py Outdated
Comment threadtest/grid/geometry/test_angles.py Outdated
Comment threadtest/grid/geometry/test_angles.py Outdated

@cmdupuis3cmdupuis3 left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

LGTM, tests pass and the tests and numba parallelization look reasonable. The vector arithmetic will probably need to sync up with other incoming PRs in the future though

@Sevans711
Sevans711 merged commit ab39109 into mainAug 18, 2026
13 checks passed
@Sevans711
Sevans711 deleted the face_node_angles branch August 18, 2026 18:38
Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Labels

new featureNew feature or requestrun-benchmarkRun ASV benchmark workflow

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Compute face_node_angles from uxarray.Grid?

3 participants

@Sevans711@erogluorhan@cmdupuis3