Repository navigation
Conversation
… weights _compute_band_overlap_area rebuilt the face/band intersection polygon by sorting its points by longitude, which does not recover the boundary order. Partial-face overlap areas were off by ~1.3% median (up to 3.4%) on outCSne30 and far worse near the poles, so conservative zonal_mean and zonal_anomaly weights failed the spherical-zone-area identity. Replace the polygon reconstruction with an exact line integral. In the equal-area (lon, sin(lat)) coordinates a latitude band is a strip, so by Green's theorem the overlap area is -closed_integral(clip(z, z_min, z_max) - z_min) dlon around the face boundary. On each great-circle edge, split at the band crossings, the integral has a closed form. There is no vertex ordering to get wrong, the antimeridian needs no special case, and pole-enclosing faces and vertices on a pole are handled through the longitude winding. A face whose boundary never enters the band gets an exact 0 (or the full polar cap), so faces that only touch a band keep zero weight. On outCSne30 the band weights now satisfy the zone-area identity to ~2.5e-11 relative, fully-contained faces match face_areas to 1e-10, and the band-mean error of sin(lat) converges under healpix refinement. Remove the now-unused _sort_points_by_angle and add regression tests for the zone-area identity, per-face partition, contained faces, rotation and orientation invariance, pole-enclosing faces, pole vertices and touching faces. Co-Authored-By: Claude Opus 5.5 <[email protected]>
With latitude_adjusted_area=True, _face_area_from_quadrature applied the constant-latitude edge correction with the wrong sign for counterclockwise faces: it added the area between a parallel and its poleward-bulging great-circle arc where it should subtract it, and vice versa. The magnitude from _area_correction was right. For the reference triangle in the tests (top edge on the 41.8 deg N parallel) the corrected area was 0.022448 instead of 0.021884; the CORRECTED_TRI_AREA constant had been pinned to the wrong value. The new value is checked against brute-force integration and agrees with the exact band-overlap line integral. Add a test that clips a face by hand, applies calculate_face_area with latitude_adjusted_area=True, and checks it matches _compute_band_overlap_area in both hemispheres. Co-Authored-By: Claude Opus 5.5 <[email protected]>
cmdupuis3
added
accuracy
ASV BenchmarkingBenchmark Comparison ResultsBenchmarks that have stayed the same:
Benchmarks that have got worse:
|
|
thanks for the work @cmdupuis3 I will look into it in details later this week. |
| """Sort points by longitude angle for proper polygon formation. | ||
| def _is_pole_point(p): | ||
| """Whether ``p`` lies (numerically) on a pole, where longitude is undefined.""" | ||
| return p[0] * p[0] + p[1] * p[1] < ERROR_TOLERANCE * ERROR_TOLERANCE |
There was a problem hiding this comment.
This checks only p[0]**2 + p[1]2 < tol2, i.e. that p is near the z-axis. That counts as “on a pole” only if p lies exactly on the unit sphere. Floating-point points generally don’t, and newly constructed points (intersections, interpolated points) may already carry error. The function also has no guard and no documented precondition, so it is easy to misuse elsewhere.
If the intended meaning really is just “x and y are ~0, so longitude is undefined”, state that in the docstring along with the precondition (unit-norm input, and the tolerance it assumes). Keep it private to this module so it isn’t reused for general pole tests.
Add tests with points slightly off the unit sphere and points near each pole.
|
|
||
| while j >= 0 and angles[j] > key_angle: | ||
| angles[j + 1] = angles[j] | ||
| sorted_points[j + 1] = sorted_points[j] | ||
| j -= 1 | ||
| @njit(cache=True) | ||
| def _gca_z_dlon_integral(a, b): | ||
| """Exact value of the line integral of ``z dlon`` along the minor great-circle arc a -> b. |
There was a problem hiding this comment.
I thought we already have this function implemented
Line 17 in 0287eb0
There was a problem hiding this comment.
The existed solution already tried to avoid the catastrophic cancellation as much as possible, while in this solution the following lines
ta_z = n[0] * a[1] - n[1] * a[0]
tb_z = n[0] * b[1] - n[1] * b[0]
and (tb_z - ta_z) all shows potential risks of the catastrophic cancellation
There was a problem hiding this comment.
consider using sth like this
from uxarray.utils.computing import diff_of_products
from uxarray.grid.utils import _small_angle_of_2_vectors
@njit(cache=True)
def _lon_delta(a, b):
"""Signed longitude change from ``a`` to ``b`` in (-pi, pi]."""
cz_hi, cz_lo = diff_of_products(a[0], b[1], a[1], b[0])
a_xy = (a[0], a[1], 0.0)
b_xy = (b[0], b[1], 0.0)
return np.sign(cz_hi + cz_lo) * _small_angle_of_2_vectors(a_xy, b_xy)| ndarray | ||
| Points sorted by longitude angle | ||
| @njit(cache=True) | ||
| def _lon_delta(a, b): |
There was a problem hiding this comment.
_lon_delta is just the angle between a and b after both are slid along their meridians to the equator. Something like the following should do the job
@njit(cache=True)
def _lon_delta(a, b):
"""Signed longitude change from a to b in (-pi, pi].
Precondition: neither a nor b is at a pole."""
# slide both points to the equator (drop z)
a_eq = (a[0], a[1], 0.0)
b_eq = (b[0], b[1], 0.0)
# magnitude: Kahan angle between the two equator points
angle = _small_angle_of_2_vectors(a_eq, b_eq)
# sign: accurate z-component of a_eq x b_eq (positive = eastward)
s_hi, s_lo = diff_of_products(a[0], b[1], a[1], b[0])
return np.sign(s_hi + s_lo) * angle| angles[j + 1] = key_angle | ||
| sorted_points[j + 1] = key_point | ||
| @njit(cache=True) | ||
| def _edge_band_integral(a, b, z_min, z_max): |
There was a problem hiding this comment.
I agree with the mathematical treatment here: the Green's theorem formulation in (λ, z = sin φ), together with the LCL terms, reduces to the GCA area plus LCL area correction in Chen et al. (2026, GMD), so mathematically it is correct. However, several issues remain. First, the floating-point accuracy handling: the cross products, intersection, and angle computations are evaluated naively instead of with the accurate kernels already in uxarray. Second, the spherical geometry and topology handling。
| b = face_edges_xyz[e, 1] | ||
| z_sum += a[2] | ||
|
|
||
| if _is_pole_point(a) or _is_pole_point(b): |
There was a problem hiding this comment.
I understand the motivation behind the longitude-based winding-number calculation and the additional corrections for north- and south-pole configurations. However, I think we should reconsider this approach and reuse the existing geometric predicates in UXarray.
The incremental-angle algorithm for computing winding numbers is already well established (Weiler, 1994). As discussed in our paper, Chen et al. (2026), Accurate and robust geometric algorithms for regridding on the sphere, for the class of simple, non-self-intersecting faces typically used in geoscientific meshes, winding-number point-in-face testing can be reduced to an equivalent ray-casting even–odd test (Hormann and Agathos, 2001).
In general, we do not recommend implementing winding-number tests using the incremental-angle formulation, since the repeated trigonometric evaluations are computationally expensive and can introduce additional floating-point inaccuracies. Orientation-based ray-casting predicates are generally preferable for these mesh configurations.
I understand that the current implementation attempts to account for the different orientations around the north and south poles. However, the longitude-based winding calculation does not directly distinguish between the two poles, which introduces additional classification and correction logic. In particular, north = z_sum > 0 is not a general geometric criterion for pole containment, since the sum of vertex z-coordinates depends on the number and distribution of vertices and is not invariant under GCA edge subdivision.
More importantly, I don't think we need to introduce another pole-containment calculation here.
In the existing UXarray workflow, face latitude–longitude bounds are already constructed before the overlap-area calculation. Specifically, uxarray/grid/bounds.py already evaluates:
has_north_pole = pole_point_inside_polygon(
1, face_edges_xyz, face_edges_lonlat
)
has_south_pole = pole_point_inside_polygon(
-1, face_edges_xyz, face_edges_lonlat
)Furthermore, _compute_face_band_weights() accesses uxgrid.face_bounds_lat to identify candidate faces before calling _compute_band_overlap_area().
Therefore, the pole-containment predicates have already been evaluated during face-bounds construction. The resulting latitude–longitude bounds also encode relevant polar information, such as latitude bounds reaching ±90° and longitude bounds spanning the full 360°.
I would suggest checking whether we can directly reuse this existing information to determine pole containment. If the current bounds representation does not preserve the distinction between pole containment and merely touching a pole, we could retain the corresponding boolean flags during bounds construction rather than recomputing them.
This would allow us to preserve the proposed Green's-theorem-based area integration while avoiding a redundant, coordinate-dependent winding-number implementation and maintaining consistency with the existing spherical geometry framework.
For a detailed discussion of winding-number algorithms, their relationship to ray casting, and the associated geometric predicates, please refer to 5.6 Face-level: Point-in-face predicate on the spherehttps://egusphere.copernicus.org/preprints/2026/egusphere-2026-636/
| r = np.hypot(a[2], u[2]) | ||
| t0 = np.arctan2(u[2], a[2]) | ||
|
|
||
| ts = np.empty(6, dtype=np.float64) |
There was a problem hiding this comment.
I would strongly recommend not reinventing the wheel here. UXarray already has an established
uxarray/uxarray/grid/intersections.py
Line 507 in 6d6c6da
uxarray/uxarray/grid/intersections.py
Line 594 in 6d6c6da
There are three main concerns with the current approach.
First, mixing angular and Cartesian representations in geometric calculations introduces unnecessary complexity and potential numerical accuracy issues, particularly when repeatedly converting between the two representations.
Second, as I mentioned previously, we generally try to avoid trigonometric functions in geometric predicates and intersection calculations whenever possible. Functions such as atan2, acos, sin, and cos are computationally expensive and introduce additional numerical errors. In particular, computing the intersection parameter using acos(c / r) is especially sensitive near tangential intersections, where c / r approaches ±1. This can affect not only the accuracy of the intersection coordinates but also the correctness of the intersection classification.
Third, the existing GCA–LCL intersection algorithm was specifically designed to support SIMD-friendly vectorization and parallelization. The proposed implementation introduces multiple trigonometric evaluations and conditional branches into the intersection calculation. Transcendental functions are generally more expensive and harder to vectorize efficiently than basic Cartesian arithmetic, while extensive branching can further reduce SIMD efficiency due to divergent execution paths. Consequently, this implementation risks sacrificing both computational performance and numerical robustness compared with the existing algorithm.
We can preserve the proposed Green's-theorem-based integration while reusing the established geometric intersection routine rather than maintaining two independent implementations of the same operation.
2 participants
Closes #1644
Overview
1. Exact face/latitude-band overlap areas (
uxarray/core/zonal.py)_compute_band_overlap_arearebuilt the face/band intersection polygon by sorting its points by longitude, which loses boundary order. This made partial-face weights wrong inzonal_mean(..., conservative=True)andzonal_anomaly(..., conservative=True). OnoutCSne30the error was a median of 1.3% and up to 3.4%. Near the poles it was much larger: −54% for the 85°–90° band on HEALPix z3.The overlap area is now computed exactly, without building a polygon:
Clipping plus latitude correction was considered and rejected. It can't represent faces that enclose a pole, which are common on MPAS, HEALPix and cubed-sphere grids. On ordinary faces it computes the same quantity, and a test asserts the two agree.
Results:
outCSne30)face_areasIt's also about 2.5× faster.
2. Sign of the constant-latitude area correction (
uxarray/grid/area.py)Found while adding the equivalence test. With
latitude_adjusted_area=True, the correction for parallel edges was applied with the wrong sign for counterclockwise faces; its size was right.CORRECTED_TRI_AREAwas pinned to the wrong value (0.022448 → 0.021884, verified by brute-force integration). This changesGrid.calculate_total_face_area(latitude_adjusted_area=True)results.Note: HEALPix weighting (not changed here)
HEALPix
face_areasare the theoretical equal areas 4π/N, but uxarray's great-circle polygons for those pixels differ from 4π/N by up to 10% near the poles, at every resolution. Within a band, fully-contained faces carry 4π/N while partial faces carry polygon geometry. So the zone-area identity holds only approximately on HEALPix: 0.004% for 20°–60° at z5, 1.7% for 85°–90°. One option is to scale each partial overlap byface_areas[f]divided by that face's polygon area. That changes what the weights mean, so it's left for a follow-up.🤖 Generated with Claude Code
PR Checklist
General
Testing & Benchmarking
Documentation and Examples
docs/api.rst; internal (private) function names start with an underscore (_)AI Disclosure
AI Usage: Claude Opus 5.5