Skip to content

warn if calculate_total_face_area() too large; add FESOM areas regression test #1649

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Sevans711 wants to merge 3 commits into from
+24 −4
Open

warn if calculate_total_face_area() too large; add FESOM areas regression test #1649

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
09e13ee
face_areas fesom test, & warning if total too big
Sevans711 Aug 5, 2026
51efa7d
Merge branch 'main' into sevans/area-error-checks
Sevans711 Aug 6, 2026
1a8c198
Merge branch 'main' into sevans/area-error-checks
Sevans711 Aug 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion test/grid/grid/test_areas.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,6 @@ def test_latlon_bounds_populate_bounds_MPAS(gridpath):
uxgrid = ux.open_grid(gridpath("mpas", "QU", "oQU480.231010.nc"))
bounds_xarray = uxgrid.bounds


def _sum_quadrature_jacobians(x, y, z, quadrature_rule, order):
"""Independently sum the Jacobian at each quadrature point of a triangle.

Expand Down Expand Up @@ -163,3 +162,10 @@ def test_calculate_face_area_jacobian_is_quadrature_sum(quadrature_rule, order):

assert jacobian > 0
nt.assert_allclose(jacobian, expected, rtol=1e-12)


def test_face_areas_fesom(gridpath):
"""Ensure correct total area for FESOM grid (~8.3780 sr). Regression test for #425."""
uxgrid = ux.open_grid(gridpath("ugrid", "fesom", "fesom.mesh.diag.nc"))
total_area = uxgrid.calculate_total_face_area()
nt.assert_almost_equal(total_area, 8.3780, decimal=4)
20 changes: 17 additions & 3 deletions uxarray/grid/grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1958,7 +1958,12 @@ def calculate_total_face_area(
"""Calculate the total surface area of all the faces in a mesh.

Equivalent to ``self.compute_face_areas(...).sum()``; provided as a
convenience.
convenience. (Note: for HEALPix grids, when called with default arguments,
this method actually returns ``self.face_areas.sum()`` instead,
which respects HEALPix equal-area property.)

Additionally, raises a warning if the result is larger than
the total area of a sphere (4 * pi * self.sphere_radius**2).

Parameters
----------
Expand Down Expand Up @@ -1986,16 +1991,25 @@ def calculate_total_face_area(
and order == 4
and not latitude_adjusted_area
):
return np.sum(self.face_areas.values)
result = np.sum(self.face_areas.values)

return np.sum(
result = np.sum(
self.compute_face_areas(
quadrature_rule=quadrature_rule,
order=order,
latitude_adjusted_area=latitude_adjusted_area,
)
)

RTOL = 1e-6 # 1e-7 had warnings in existing CI tests (as of 2026-08-05), 1e-6 did not.
if result > 4 * np.pi * self.sphere_radius**2 * (1 + RTOL):
Comment on lines +2004 to +2005

@cmdupuis3 cmdupuis3 Aug 10, 2026

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.

Wouldn't it make sense to also check if result < ... * (1 - RTOL)?

@Sevans711 Sevans711 Aug 10, 2026 •
edited
Loading

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

That check would be like asking "is the total area less than the sphere's surface area?", right?

I wouldn't want a check like that because grids are always free to only cover a smaller region or have holes (e.g., see example in discussion of originally-linked issue).

@cmdupuis3 cmdupuis3 Aug 10, 2026

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.

Yeah. What I mean is that depending on how you approximate the sphere, you may be able to say something about the bounds, so it could be the case that this approximation is strictly larger than the sphere's total area. But if not, we should have tolerance on the lower bound also. But yeah, if we're not worried about preserving area on the lower bound, then we don't need to worry about it.

@Sevans711 Sevans711 Aug 10, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

I think we aren't worried about the lower bound here, I'm mostly understanding the upper bound warning as a quick way to flag the possibility of a malformed grid.

Yeah. What I mean is that depending on how you approximate the sphere, you may be able to say something about the bounds, so it could be the case that this approximation is strictly larger than the sphere's total area. But if not, we should have tolerance on the lower bound also.

I'm not sure what this means, could you clarify a bit further (if still relevant, assuming you agree it is okay to not check a lower bound)?

@cmdupuis3 cmdupuis3 Aug 10, 2026

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.

So like, consider if you were trying to build a grid on the inside surface of the sphere, by placing nodes on the sphere. The area would approximate the sphere's area, but it would inherently be slightly smaller because it's strictly inside. In that case you'd expect the approximate area to be within (1 - RTOL) * sphere area. If we put faces tangent to the sphere, you'd have the opposite situation because the nodes are inherently outside of the sphere. What I'm getting at is depending on the method, you could have some mix of the two, or maybe you have some other way of knowing if the grid area is supposed to be strictly larger or smaller than the sphere.

@Sevans711 Sevans711 Aug 10, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

Ahh, thank you, that helped me understand. I think that could be good to consider in the future, but isn't as relevant here. Here, maybe a clearer way to phrase the intention of the check would be to write something like:

full_sphere_area = 4 * np.pi * self.sphere_radius**2
assert (result < full_sphere_area) or np.isclose(result, full_sphere_area, rtol=RTOL)

I.e., the result is definitely small enough (less than full sphere area) or maybe this is a grid that covers the whole sphere (but account for possible rounding errors in calculations). Right now I am planning to keep the underlying code unchanged, but please let me know if you would want to see it rewritten or maybe with an extra comment added to clarify. (Note, the docstring does already include: Additionally, raises a warning if the result is larger than the total area of a sphere (4 * pi * self.sphere_radius**2).)

warnings.warn(
f"Total face area (={result}) exceeds the surface area of the whole sphere "
f"(={4 * np.pi * self.sphere_radius**2}) (with sphere_radius={self.sphere_radius}).",
)

return result

def compute_face_areas(
self,
quadrature_rule: str = "triangular",
Expand Down