Clean up face_areas API

[Sevans711]

Proposed new feature or change:

The public-facing API for Grid face areas is now inconsistent in its support for quadrature kwargs like quadrature_rule and order. The calculate_total_face_area() method supports these kwargs, but the only public way to get face areas (Grid.face_areas property) does not support them. Related: #1567, #1568.

Suggested fix: un-deprecate and improve the public-facing compute_face_areas() method. This would provide public support for quadrature kwargs. The public method was originally deprecated in #1379; pinging @rajeeja and @erogluorhan who worked on that PR. More thoughts here:

• Currently, that method returns (face_areas, jacobian); I would suggest to instead just return face_areas by default, since the method name does not mention jacobian. We could adopt numpy style (e.g. see docstring for np.unique) by adding an optional kwarg like def compute_face_areas(..., return_jacobian=False).
• (new feature) Maybe promote the return type to be UxDataArray (with self as the result.uxgrid) by default, or at least provide a flag to enable that, instead of returning numpy array? This would make it easy to visualize the result!
• Remove _compute_face_areas() (use public method instead) or rename it to _compute_face_areas_and_jacobian() since it also returns the jacobian.

Additional suggestions / questions:

1. Keep the note in compute_face_areas() docstring suggesting to use the face_areas property instead for caching; clarify this note is only relevant if using default kwargs.
2. Clarify in face_areas property docstring that compute_face_areas() can be used to get result with adjusted quadrature behavior.
3. Improve documentation for kwargs: quadrature_rule, order, and latitude_adjusted_area. Minimum: quadrature_rule and order should specify what options are supported. Ideal: explain what each of these things actually means, either directly or by linking to a paper or online documentation somewhere. E.g., what is a quadrature rule? What does "corrects the area of the faces accounting for lines of constant latitude" actually mean?
4. Why do we provide calculate_total_face_area() when it is just np.sum(grid.compute_face_areas()) (or np.sum(grid.face_areas) if using default kwargs)? For me personally, the easiest thought process goes like: "I want total face area --> I should get face areas, then sum them." The existence of calculate_total_face_area() adds work, because when I see it I also think "wait, should I be using that method instead? Is there some quirk to getting the total area such that np.sum() isn't good enough?" until I inspect and see it really is just sum(). I would suggest to remove this method, or at least add near the top of its docstring something like "equivalent to grid.compute_face_areas().sum()"
5. Should the method get_all_face_area_from_coords be marked as private (i.e., leading underscore)? It doesn't seem to appear on the public uxarray docs, and it seems to be very low-level (inputs are x, y, z, face_nodes, face_geometry). Similarly, should the other methods in grid/area.py by marked as private?

[rajeeja]

Yes — it has been a while. I wrote most of this a few years ago and it is time this module got a rejig. Thanks for the detailed issue.

From what I can see, the plan of action should include:

1. 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. Rename _compute_face_areas → _compute_face_areas_and_jacobian so the name matches what it returns, and keep it as the internal worker.
3. Add an optional as_dataarray=False flag to return a UxDataArray (result.uxgrid = self) so per-face areas are easy to plot.
4. Keep calculate_total_face_area() but document it as equivalent to compute_face_areas().sum() now that Honor quadrature kwargs in calculate_total_face_area #1568 has made its quadrature kwargs actually take effect (previously they were silently ignored).
5. Mark the low-level area.py routines private (get_all_face_area_from_coords, the jacobian/quadrature helpers) with leading underscores, since they are not part of the public surface.
6. Vectorize the core computation. get_all_face_area_from_coords currently loops over every face in Python-level serial code inside @njit; it should be parallelized (prange) and made vectorization-friendly. This is the hot path on large meshes and the main performance win in this cleanup, so it should be a first-class goal alongside the API tidy-up, not an afterthought.
7. Fix the user guide (area_calc.ipynb). It currently instructs users to call the private _compute_face_areas() for custom quadrature, and its order=1 vs order=4 narrative predates Honor quadrature kwargs in calculate_total_face_area #1568 — both need updating against the new public method and the fixed kwarg behavior.

I'll take this one and assign it to myself.

[Sevans711]

No problem, and thank you for the detailed response and taking on this issue! That all sounds good to me. I think your plan includes basically everything I was thinking of, plus some extra (especially the possible speedup!). Please keep me posted on relevant changes / PR; I'd be happy to take a close look whenever things are ready for review.

The only point I can think to add to your action plan would be:
8. improve docstrings for face_areas property and compute_face_areas() (as in my original points (2) & (3)).