Skip to content

Clean up and vectorize the face_areas API#1577

Open
rajeeja wants to merge 9 commits into
mainfrom
rajeeja/area_cleanup
Open

Clean up and vectorize the face_areas API#1577
rajeeja wants to merge 9 commits into
mainfrom
rajeeja/area_cleanup

Conversation

@rajeeja

@rajeeja rajeeja commented Jul 15, 2026

Copy link
Copy Markdown
Contributor

Closes #1571.

  1. Un-deprecate compute_face_areas() as the public quadrature-capable entry point — returns areas by default, with return_jacobian and as_uxarray flags (the latter returns a UxDataArray paired with the grid).
  2. Rename the internal worker to _compute_face_areas_and_jacobian and privatize the low-level area.py helpers.
  3. Document the quadrature_rule/order/latitude_adjusted_area options; note calculate_total_face_area is equivalent to compute_face_areas().sum().
  4. Vectorize get_all_face_area_from_coords with prange and hoist the per-face quadrature setup out of the loop — ~6x faster with identical areas; adds a FaceAreas asv benchmark.
  5. Plot face areas in area_calc.ipynb and the geometric/theoretical ratio in healpix.ipynb using the new as_uxarray flag.

Un-deprecate compute_face_areas() as the public quadrature-capable entry
point (areas by default, return_jacobian and as_dataarray flags), rename
the internal worker to _compute_face_areas_and_jacobian, privatize the
low-level area.py helpers, and document the quadrature options. Vectorize
get_all_face_area_from_coords with prange and hoist the per-face
quadrature setup out of the loop for a 6.3x speedup with identical areas.
Closes #1571.
@review-notebook-app

Copy link
Copy Markdown

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@Sevans711 Sevans711 added documentation Improvements or additions to documentation redesign Content relating to the ongoing redesign improvement Improvements on existing features or infrastructure run-benchmark Run ASV benchmark workflow labels Jul 15, 2026
@Sevans711
Sevans711 self-requested a review July 15, 2026 23:16
@Sevans711

Copy link
Copy Markdown
Collaborator

Thank you for these changes! I added run-benchmarks label but then realized I'm not sure if there are any benchmarks which actually rely on face_areas. If there's no visible impact, then would you be able to add at least one face_areas benchmark? If the impact is visible already, I don't feel strongly about needing a face_areas specific benchmark.

I'll come back to this soon (probably tomorrow), and if the benchmark impact is visible then I will review the rest!

@github-actions

github-actions Bot commented Jul 15, 2026

Copy link
Copy Markdown

ASV Benchmarking

Benchmark Comparison Results

Benchmarks that have improved:

Change Before [a5715ef] After [f111ad7] Ratio Benchmark (Parameter)
- 580M 391M 0.67 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
- 698M 392M 0.56 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
- 487M 385M 0.79 mpas_ocean.FaceAreas.peakmem_compute_face_areas('480km')
- 156±1ms 57.4±0.1ms 0.37 mpas_ocean.FaceAreas.time_compute_face_areas('120km')
- 12.1±0.03ms 5.96±0.02ms 0.49 mpas_ocean.FaceAreas.time_compute_face_areas('480km')
- 489M 384M 0.79 mpas_ocean.Gradient.peakmem_gradient('480km')

Benchmarks that have stayed the same:

Change Before [a5715ef] After [f111ad7] Ratio Benchmark (Parameter)
10.2±0.1μs 9.75±0.07μs 0.96 bench_connectivity.Connectivity.time_edge_face('120km')
10.2±0.04μs 10.1±0.06μs 0.99 bench_connectivity.Connectivity.time_edge_face('480km')
10.2±0.2μs 9.87±0.05μs 0.97 bench_connectivity.Connectivity.time_edge_node('120km')
10.3±0.2μs 10.1±0.04μs 0.99 bench_connectivity.Connectivity.time_edge_node('480km')
10.4±0.3μs 10.3±0.5μs 0.99 bench_connectivity.Connectivity.time_face_edge('120km')
10.2±0.07μs 10.2±0.04μs 1.01 bench_connectivity.Connectivity.time_face_edge('480km')
9.99±0.09μs 10.2±0.4μs 1.02 bench_connectivity.Connectivity.time_face_face('120km')
10.2±0.1μs 10.1±0.06μs 0.99 bench_connectivity.Connectivity.time_face_face('480km')
20.4±0.3μs 20.4±0.1μs 1 bench_connectivity.Connectivity.time_face_node('120km')
21.3±0.2μs 21.0±0.7μs 0.99 bench_connectivity.Connectivity.time_face_node('480km')
9.91±0.05μs 9.97±0.1μs 1.01 bench_connectivity.Connectivity.time_node_edge('120km')
10.4±0.07μs 10.2±0.1μs 0.99 bench_connectivity.Connectivity.time_node_edge('480km')
9.93±0.08μs 9.82±0.05μs 0.99 bench_connectivity.Connectivity.time_node_face('120km')
10.2±0.08μs 10.3±0.1μs 1 bench_connectivity.Connectivity.time_node_face('480km')
389M 389M 1 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
421M 421M 1 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
12.8±0.03ms 13.1±0.3ms 1.02 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
3.43±0.07ms 3.45±0.09ms 1.01 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
17.3±0.09ms 17.3±0.09ms 1 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
2.01±0.04ms 1.99±0.01ms 0.99 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
989±5ms 989±5ms 1 import.Imports.timeraw_import_uxarray
901±6ns 889±6ns 0.99 mpas_ocean.CheckNorm.time_check_norm('120km')
849±10ns 855±5ns 1.01 mpas_ocean.CheckNorm.time_check_norm('480km')
714±6ms 716±5ms 1 mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')
44.8±0.1ms 45.3±0.3ms 1.01 mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')
629±5μs 630±6μs 1 mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km')
557±6μs 555±10μs 1 mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')
5.03±0.03ms 5.10±0.08ms 1.01 mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('120km')
3.55±0.03ms 3.55±0.03ms 1 mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('480km')
3.33±0.01s 3.35±0.01s 1.01 mpas_ocean.ConstructFaceLatLon.time_welzl('120km')
212±0.8ms 213±1ms 1 mpas_ocean.ConstructFaceLatLon.time_welzl('480km')
19.7±0.03ms 19.7±0.04ms 1 mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')
1.12±0.02ms 1.11±0.01ms 0.99 mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')
10.6±0.02ms 10.6±0.02ms 1 mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')
731±20μs 759±20μs 1.04 mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')
652±1ms 655±1ms 1 mpas_ocean.CrossSections.time_const_lat('120km', 1)
330±2ms 331±2ms 1 mpas_ocean.CrossSections.time_const_lat('120km', 2)
171±0.4ms 170±0.7ms 1 mpas_ocean.CrossSections.time_const_lat('120km', 4)
486±3ms 485±4ms 1 mpas_ocean.CrossSections.time_const_lat('480km', 1)
246±0.6ms 248±0.5ms 1.01 mpas_ocean.CrossSections.time_const_lat('480km', 2)
128±0.6ms 127±0.2ms 1 mpas_ocean.CrossSections.time_const_lat('480km', 4)
22.0±0.08ms 22.6±0.2ms 1.02 mpas_ocean.DualMesh.time_dual_mesh_construction('120km')
2.57±0.02ms 2.59±0.03ms 1.01 mpas_ocean.DualMesh.time_dual_mesh_construction('480km')
399M 400M 1 mpas_ocean.FaceAreas.peakmem_compute_face_areas('120km')
844±2ms 836±4ms 0.99 mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)
50.3±1ms 51.3±2ms 1.02 mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)
73.2±0.5ms 73.4±0.3ms 1 mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)
5.53±0.1ms 5.57±0.07ms 1.01 mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)
404M 404M 1 mpas_ocean.Gradient.peakmem_gradient('120km')
165±0.4ms 166±1ms 1 mpas_ocean.Gradient.time_gradient('120km')
11.6±0.3ms 11.5±0.04ms 0.98 mpas_ocean.Gradient.time_gradient('480km')
189±1μs 189±0.7μs 1 mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('120km')
93.5±5μs 87.6±0.7μs 0.94 mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('480km')
351M 351M 1 mpas_ocean.Integrate.peakmem_integrate('120km')
331M 330M 1 mpas_ocean.Integrate.peakmem_integrate('480km')
181±1μs 184±2μs 1.02 mpas_ocean.Integrate.time_integrate('120km')
162±2μs 156±1μs 0.96 mpas_ocean.Integrate.time_integrate('480km')
187±2ms 186±0.7ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude')
185±0.8ms 186±2ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include')
186±0.9ms 186±1ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')
13.7±0.06ms 13.8±0.07ms 1.01 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude')
13.6±0.05ms 13.7±0.07ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include')
13.7±0.1ms 13.7±0.05ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')
286±1μs 286±1μs 1 mpas_ocean.PointInPolygon.time_face_search_lonlat('120km')
287±1μs 287±0.5μs 1 mpas_ocean.PointInPolygon.time_face_search_lonlat('480km')
267±0.7μs 268±2μs 1 mpas_ocean.PointInPolygon.time_face_search_xyz('120km')
270±2μs 271±1μs 1 mpas_ocean.PointInPolygon.time_face_search_xyz('480km')
213±0.6ms 213±1ms 1 mpas_ocean.RemapDownsample.time_bilinear_remapping
229±0.5ms 230±1ms 1 mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping
4.25±0.02ms 4.22±0.05ms 0.99 mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping
1.19±0.01s 1.19±0.01s 1 mpas_ocean.RemapUpsample.time_bilinear_remapping
35.6±0.4ms 36.4±0.4ms 1.02 mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping
8.82±0.2ms 9.04±0.06ms 1.02 mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping
28.9±0.5ms 29.0±0.3ms 1 mpas_ocean.ZonalAverage.time_zonal_average('120km')
5.96±0.02ms 6.03±0.05ms 1.01 mpas_ocean.ZonalAverage.time_zonal_average('480km')
326M 330M 1.01 quad_hexagon.QuadHexagon.peakmem_open_dataset
325M 326M 1 quad_hexagon.QuadHexagon.peakmem_open_grid
6.55±0.1ms 6.50±0.1ms 0.99 quad_hexagon.QuadHexagon.time_open_dataset
5.49±0.09ms 5.54±0.07ms 1.01 quad_hexagon.QuadHexagon.time_open_grid

@Sevans711 Sevans711 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.

Looking through the plan of action discussed in #1571:

  1. (Looks good to me) Un-deprecate and make compute_face_areas() the public, kwargs-capable entry point. Return face_areas only by default, with return_jacobian=False (numpy-style) for callers that also want the jacobian.
  2. (Looks good to me) Rename _compute_face_areas_compute_face_areas_and_jacobian
  3. (Not completed as planned; see notes below) Add an optional as_dataarray=False flag to return a UxDataArray (result.uxgrid = self) so per-face areas are easy to plot.
  4. (Looks good to me) Keep calculate_total_face_area() but document it as equivalent to compute_face_areas().sum().
  5. (Looks good to me) Mark the low-level area.py routines private
  6. (Mostly good but I have more notes below) Vectorize the core computation.
  7. (Mostly good but I have more notes below) Fix the user guide (area_calc.ipynb).
  8. (Looks good, I did a small commit to help) improve docstrings for face_areas property and compute_face_areas().

More notes / requests / suggestions:

  • (3) Originally we discussed returning a UxDataArray object, not an xarray.DataArray. The current version here returns an xarray.DataArray. Request: make the optional flag return UxDataArray instead. I think UxDataArray is the way to go; it is more convenient, and I think the whole point for this optional flag is to return something more convenient.
    • For example, plotting result from current version looks something like: ux.UxDataArray(uxds.uxgrid.compute_face_areas(as_dataarray=True), uxgrid=uxds.uxgrid).plot(). But, with this requested change, it would become simpler: uxds.uxgrid.compute_face_areas(as_dataarray=True).plot().
    • I wonder if as_dataarray is the wrong name here... Request/suggestion: rename the flag to as_uxarray or as_uxdataarray to clarify it is a UxDataArray instead? Of the three names, my preferences would be as_uxarray (concise & clear) > as_uxdataarray (verbose & clear) > as_dataarray (sounds like xarray not uxarray). (Let me know if you want me to take on the task of setting this up!)
  • (7a) Request: plot face_areas somewhere in area_calc.ipynb, using the new optional flag. Either in a new section, or inside an existing section, e.g. maybe it makes sense to add in the "Calculate Area from Multiple Faces in Spherical Coordinates" section?
  • (7b) In the healpix.ipynb guide, there is a section comparing geometric face area calculations versus theoretical healpix areas. This has been updated to now be equivalent to grid.compute_face_faces() / grid.face_areas. There is also some numerical analysis of how the results vary across the grid. Request: plot this face_areas ratio across the grid. After making changes in response to point (3) above, it should be as simple as: (grid.compute_face_faces(as_uxarray=True) / grid.face_areas).plot().
  • (6) I also measured roughly 6x faster face_areas computations on main versus on this branch. This is a great improvement! Some of the ASV benchmarks show improvements in reducing peak memory usage, but none of them show this significant speedup. Request: add at least one benchmark that gets improved by these changes.

@erogluorhan erogluorhan changed the title Clean up and vectorize the face_areas API (#1571) Clean up and vectorize the face_areas API Jul 17, 2026
@Sevans711

Copy link
Copy Markdown
Collaborator

@rajeeja thank you for adding the total face area error check. I don't want to interrupt you if you're still working on things, but wanted to share my thoughts now while I am thinking about it:

  1. Maybe this check should only be in calculate_total_face_areas(), or in grid.validate()? It seems unusual to impose inefficiency on the low-level _compute_face_areas_and_jacobian method for a validation check.
  2. Flagging that face_areas will be much much larger than 4 * pi for any grids which are not on the unit sphere! From skimming your changes, it is unclear to me whether that would trigger this warning or not. Example non-unit-sphere-grid: 'uxarray/test/meshfiles/mpas/dyamond-30km/gradient_grid_subset.nc'

compute_face_areas() gains as_uxarray (renamed from as_dataarray) returning
a UxDataArray paired with the grid; plot face areas in area_calc.ipynb and
the geometric/theoretical ratio in healpix.ipynb using it; add a FaceAreas
asv benchmark capturing the vectorization speedup.
@rajeeja

rajeeja commented Jul 17, 2026

Copy link
Copy Markdown
Contributor Author

@Sevans711 thanks for the review — fixed the items: compute_face_areas() now returns a UxDataArray via as_uxarray=True (renamed from as_dataarray); face areas plotted in area_calc.ipynb and the geometric/theoretical ratio in healpix.ipynb using it; and added a FaceAreas asv benchmark that captures the speedup.

@rajeeja
rajeeja requested a review from Sevans711 July 17, 2026 18:54
@rajeeja

rajeeja commented Jul 17, 2026

Copy link
Copy Markdown
Contributor Author

@Sevans711 Good catches on both. You're right the check misfires on grids whose face_areas are parsed from the file in physical units — e.g. gradient_grid_subset.nc returns ~1.5e11 m² and would falsely warn, since the property returns file values, not unit-sphere steradians. There's no reliable way to know the unit at that level, so I've removed the check. #425's actual issue (the wrong value) is already fixed independently; I'll leave #425 open for a proper unit-aware validation later rather than ship a check that false-positives.

@Sevans711 Sevans711 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.

Looks good! Thank you for implementing all the requests/suggestions I made. I think it is ready to merge!

I did notice a few small issues with type hinting. We could either fix them as part of this PR, or open a separate issue with a more narrowly-scoped PR to avoid blocking this one. Both of these issues already exist on main; they aren't being newly introduced by this PR.

  1. Type hinting for the result of compute_face_areas() is missing.
  2. Type hinting for input arguments is a bit too permissive. For example, it includes order: int | None = 4, but there is no support for order=None, so I think it should be order: int = 4 instead.

Sidenote, I especially enjoyed seeing the face_areas plot in the Face Area Calculations docs page, I think it looks pretty cool:

Image

@rajeeja

rajeeja commented Jul 20, 2026

Copy link
Copy Markdown
Contributor Author

Thanks for the approval @Sevans711! On the type-hint items — since both issues already exist on main and weren't introduced here, I'd rather not hold up this PR for them. Opened #1590 to track the cleanup separately. Happy to merge whenever Orhan takes a look.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

documentation Improvements or additions to documentation improvement Improvements on existing features or infrastructure redesign Content relating to the ongoing redesign run-benchmark Run ASV benchmark workflow

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Clean up face_areas API

3 participants