[Claude] Conservative zonal mean: partial-face overlap areas are wrong by ~1-5%

[cmdupuis3]

Version

v2026.07.0

How did you install UXarray?

Source

What happened?

Incidental bug found by Claude. Basically, partial face areas are slightly off due to a normalization that was accidentally left unwired like dead code, then removed, and the remaining code subsequently evolved in a way that made the issue opaque.

---

Summary

_compute_band_overlap_area reconstructs the face/latitude-band intersection polygon by sorting its points by longitude rather than by traversal order around the polygon boundary. The resulting ring is generally not the intersection polygon, so the overlap areas used as weights by _compute_face_band_weights are wrong — by a median of 1.3% and up to 3.4% on outCSne30, and the band weights fail the analytic spherical-zone-area identity by up to 6.8%.

This affects every partially-overlapping face in UxDataArray.zonal_mean(..., conservative=True) and UxDataArray.zonal_anomaly(..., conservative=True). Fully-contained faces are unaffected — they take face_areas directly and are exact.

Measured on main @ 73101af7.

Location

• uxarray/core/zonal.py:209 — sorted_points = _sort_points_by_angle(unique_points), the call that discards boundary order.
• uxarray/core/zonal.py:92-143 — _sort_points_by_angle, which sorts by longitude about the z-axis (_small_angle_of_2_vectors(x_axis, point_xy_unit) with a y-sign flip, and angle = 0.0 for points at a pole).
• uxarray/core/zonal.py:222 — latitude_adjusted_area=True, which then applies a constant-latitude correction based on np.isclose of adjacent z values in that scrambled order.

Root cause

The intended algorithm is in the introducing commit as dead code. 2f33780 added:

# Compute centroid and sort by angle
center = np.mean(unique_points, axis=0)
center = center / np.linalg.norm(center)

angles = []
for pt in unique_points:
    # Use longitude for sorting (works well for latitude bands)
    angle = np.arctan2(pt[1], pt[0])

center is computed, normalized, and never used: the comment describes an azimuthal sort about the polygon centroid, but the code sorts by raw longitude about the z-axis. The follow-up refactor 46e0e91 removed the unused center and moved the longitude sort into _sort_points_by_angle, where it stands today.

Sorting by longitude cannot reconstruct the ring:

• For a band-clipped cell it orders points monotonically in longitude, interleaving the z = z_min and z = z_max boundary points instead of traversing lower boundary → up → upper boundary → down.
• It is degenerate across the antimeridian (the y-sign flip inherits the arctan2 branch cut).
• All points at a pole collapse to angle = 0.0.
• The "works well for latitude bands" premise fails even without clipping: on the outCSne30 band below, the sort reorders 1372 of 1380 fully-contained faces.

Evidence

A. The error is in the sort, not in the latitude correction

Ratio to the true face area, for the 1380 outCSne30 faces lying entirely inside the 20°–60° band (so the intersection polygon is the face, and the correct answer is known):

computation	median	min	max
calculate_face_area, original vertex order	1.000000	1.000000	1.000000
calculate_face_area, longitude-sorted order	1.013088	0.975159	1.024841
+ latitude_adjusted_area=True	1.013214	0.975159	1.034282
_compute_band_overlap_area (what the weights use)	1.013214	0.975159	1.034282

The sort contributes ~1.3%; the latitude correction adds ~0.01% on these faces. The latitude correction itself is sound — it is the polygon feeding it that is broken.

Note the spread straddles 1.0 (0.975 → 1.034). This is a two-sided scatter, not the one-sided bias that a deliberate approximation would produce.

B. Band weights violate the spherical-zone-area identity

A global grid tiles the sphere, so the weights of a single band must sum to 2π(sin lat₁ − sin lat₀), independent of any face-area convention:

grid	band	sum(weights)	exact	error
healpix z3	−90° → 90°	12.56637061	12.56637061	+0.0000%
healpix z3	20° → 60°	3.51665036	3.29242215	+6.8104%
healpix z3	−10° → 10°	2.18697496	2.18212736	+0.2222%
healpix z3	0° → 30°	3.15774644	3.14159265	+0.5142%
healpix z3	60° → 90°	0.84202744	0.84178721	+0.0285%
outCSne30	−90° → 90°	12.56637061	12.56637061	+0.0000%
outCSne30	20° → 60°	3.28569147	3.29242215	−0.2044%
outCSne30	−10° → 10°	2.17686972	2.18212736	−0.2409%
outCSne30	0° → 30°	3.13947029	3.14159265	−0.0676%
outCSne30	60° → 90°	0.84647905	0.84178721	+0.5574%

The sign varies with grid and band. The full-sphere rows are exact to every printed digit precisely because every face is fully contained there, so _compute_band_overlap_area is never called.

C. The error does not converge under grid refinement

Field f = sin(lat), whose exact area-weighted mean over a zone is (z₀ + z₁)/2. With 10° bands:

grid	n_face	max abs. band-mean error	max zone-area error
healpix z3	768	3.394e-03	23.601%
healpix z4	3072	3.705e-03	9.418%
healpix z5	12288	2.672e-03	5.417%

A discretization error would fall roughly 4× per refinement level. Over a 16× increase in face count the band-mean error is flat, which is the signature of a weight defect rather than resolution.

What did you expect to happen?

Suggested fix

Walk the face boundary against z = z_min and z = z_max in traversal order (Sutherland–Hodgman style). Two benefits:

1. Ring order is correct by construction, including antimeridian-spanning and pole-containing faces, which no longitude- or centroid-based sort handles robustly.
2. It records which edges lie on a constant-latitude boundary, so latitude_adjusted_area can be applied to exactly those edges instead of inferred from np.isclose on adjacent z values.

Acceptance criteria, all currently failing:

• Table B's zone-area identity holds to quadrature tolerance for every band on a global grid.
• For a face lying entirely inside a band, _compute_band_overlap_area reproduces face_areas[f]. This one also removes a discontinuity in the weight field: a face that flips between the fully-contained and partial classification currently jumps by that amount.
• Band-mean error converges under refinement.

Can you provide a MCVE to repoduce the bug?

import numpy as np
import uxarray as ux
from uxarray.core.zonal import _compute_face_band_weights

for zoom in (3, 4, 5):
    g = ux.Grid.from_healpix(zoom=zoom)
    for lat0, lat1 in [(-90.0, 90.0), (20.0, 60.0), (0.0, 30.0)]:
        (idx, w), = _compute_face_band_weights(g, np.array([lat0, lat1]))
        exact = 2 * np.pi * (np.sin(np.deg2rad(lat1)) - np.sin(np.deg2rad(lat0)))
        print(f"healpix z{zoom} band({lat0:>6},{lat1:>5}): "
              f"sum(w)={w.sum():.8f} exact={exact:.8f} err={(w.sum()/exact - 1) * 100:+.4f}%")

[cmdupuis3]

Pinging @rajeeja @hongyuchen1030

[rajeeja]

Looks legit! Let me see if the fix is valid or not, I will get back to you soon.

[hongyuchen1030]

Good work on this. The ring construction is the correct approach for constructing the new polygon representing the partially covered face.

The general workflow should be:

1. Compute the intersections between the constant-latitude line and all great-circle-arc edges of the original face. Use the resulting points to construct the new partial face, calculate its area while temporarily treating every edge as a great-circle arc, and record which edge should instead be interpreted as a constant-latitude segment.
2. Apply the area-correction formula to the recorded constant-latitude edge.

Please also add test cases that cover the acceptance criteria. At minimum, include a polygon that crosses the 0° longitude line.

[Sevans711]

Great find! Flagging that any PR fixing this should carefully consider if/how it affects the relevant docs page (https://uxarray.readthedocs.io/en/latest/user-guide/zonal-average.html) and make any updates there too, if needed