From 19be0abd4f9d7ea0db1f0f3312c33dea3041f7a9 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Mon, 5 Oct 2026 16:31:25 -0500 Subject: [PATCH 1/2] Compute exact face/latitude-band overlap areas for conservative zonal 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 --- test/grid/integrate/test_zonal.py | 124 +++++++++++- uxarray/core/zonal.py | 300 ++++++++++++++++++------------ 2 files changed, 304 insertions(+), 120 deletions(-) diff --git a/test/grid/integrate/test_zonal.py b/test/grid/integrate/test_zonal.py index 950194584..bbd69e378 100644 --- a/test/grid/integrate/test_zonal.py +++ b/test/grid/integrate/test_zonal.py @@ -400,4 +400,126 @@ def test_invalid_lat_input_raises(self): uxda.zonal_anomaly(lat=[10.0, -10.0, 30.0]) -from uxarray.core.zonal import _compute_face_band_weights # noqa: E402 +from uxarray.core.zonal import ( # noqa: E402 + _compute_band_overlap_area, + _compute_face_band_weights, +) +from uxarray.grid.area import calculate_face_area # noqa: E402 +from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array_subset # noqa: E402 + + +def _face_edges(lat, lon): + """Edge array (n_edges, 2, 3) for a face given its vertices in degrees.""" + lat, lon = np.deg2rad(lat), np.deg2rad(lon) + v = np.stack( + [np.cos(lat) * np.cos(lon), np.cos(lat) * np.sin(lon), np.sin(lat)], axis=1 + ) + return np.stack([v, np.roll(v, -1, axis=0)], axis=1) + + +def _band_area(face_edges, lat0, lat1): + return _compute_band_overlap_area( + face_edges, np.sin(np.deg2rad(lat0)), np.sin(np.deg2rad(lat1)) + ) + + +def _zone_area(lat0, lat1): + return 2 * np.pi * (np.sin(np.deg2rad(lat1)) - np.sin(np.deg2rad(lat0))) + + +class TestBandOverlapArea: + """Face/latitude-band overlap areas used as conservative zonal weights.""" + + @pytest.mark.parametrize( + "lat0, lat1", + [(20.0, 60.0), (-10.0, 10.0), (0.0, 30.0), (60.0, 90.0), (85.0, 90.0), + (-90.0, -85.0), (-37.3, 81.2), (12.0, 12.05)], + ) + def test_zone_area_identity(self, gridpath, lat0, lat1): + """The weights of one band on a global grid sum to the zone's area.""" + grid = ux.open_grid(gridpath("ugrid", "outCSne30", "outCSne30.ug")) + ((_, w),) = _compute_face_band_weights(grid, np.array([lat0, lat1])) + nt.assert_allclose(w.sum(), _zone_area(lat0, lat1), rtol=1e-9) + + def test_overlaps_partition_face_area(self, gridpath): + """A face's overlaps with a set of bands tiling the sphere sum to its area.""" + grid = ux.open_grid(gridpath("ugrid", "outCSne30", "outCSne30.ug")) + totals = np.zeros(grid.n_face) + for idx, w in _compute_face_band_weights(grid, np.arange(-90.0, 91.0, 7.5)): + totals[idx] += w + nt.assert_allclose(totals, grid.face_areas.values, rtol=1e-9) + + def test_contained_face_matches_face_area(self, gridpath): + """For a face inside the band, the overlap area is the face area.""" + grid = ux.open_grid(gridpath("ugrid", "outCSne30", "outCSne30.ug")) + faces = grid.get_faces_between_latitudes((20.0, 60.0)) + assert faces.size > 0 + n_nodes = grid.n_nodes_per_face.values + edges = _get_cartesian_face_edge_nodes_array_subset( + faces, + grid.face_node_connectivity.values, + n_nodes, + grid.n_max_face_nodes, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + areas = [ + _band_area(edges[k, : n_nodes[f]], 20.0, 60.0) for k, f in enumerate(faces) + ] + nt.assert_allclose(areas, grid.face_areas.values[faces], rtol=1e-9) + + @pytest.mark.parametrize("lat0, lat1", [(-90.0, 90.0), (12.0, 33.0), (25.0, 28.0)]) + def test_invariant_under_rotation_and_orientation(self, lat0, lat1): + """Rotating a face across the antimeridian or reversing it changes nothing.""" + lat = np.array([10.0, 15.0, 35.0, 30.0]) + lon = np.array([-8.0, 9.0, 7.0, -6.0]) + expected = _band_area(_face_edges(lat, lon), lat0, lat1) + assert expected > 0.0 + for shift in (90.0, 180.0, 183.0, 260.0): + fe = _face_edges(lat, lon + shift) + assert _band_area(fe, lat0, lat1) == pytest.approx(expected, rel=1e-12) + fe_cw = _face_edges(lat[::-1], lon[::-1] + shift) + assert _band_area(fe_cw, lat0, lat1) == pytest.approx(expected, rel=1e-12) + + @pytest.mark.parametrize("hemisphere", [1.0, -1.0]) + @pytest.mark.parametrize("reverse", [False, True]) + def test_face_enclosing_pole(self, hemisphere, reverse): + """A face around a pole (no vertex on it) contains the polar cap.""" + lat = hemisphere * np.array([80.0, 80.0, 80.0, 80.0]) + lon = np.array([45.0, 135.0, 225.0, 315.0]) + if reverse: + lat, lon = lat[::-1], lon[::-1] + fe = _face_edges(lat, lon) + x, y, z = (fe[:, 0, i] for i in range(3)) + face_area, _ = calculate_face_area(x, y, z, order=6) + + def band(a, b): + return _band_area(fe, *sorted((hemisphere * a, hemisphere * b))) + + # Edges peak at ~82.9 deg, so the 85-90 cap lies entirely inside the face. + assert band(85.0, 90.0) == pytest.approx(_zone_area(85.0, 90.0), rel=1e-12) + assert band(60.0, 75.0) == 0.0 + assert band(-90.0, 90.0) == pytest.approx(face_area, rel=1e-9) + parts = band(70.0, 81.0) + band(81.0, 84.0) + band(84.0, 90.0) + assert parts == pytest.approx(face_area, rel=1e-9) + + def test_face_with_vertex_at_pole(self): + """A face with a vertex on the pole turns there between meridians.""" + fe = _face_edges(np.array([90.0, 80.0, 80.0]), np.array([0.0, 0.0, 90.0])) + x, y, z = (fe[:, 0, i] for i in range(3)) + face_area, _ = calculate_face_area(x, y, z, order=6) + + assert _band_area(fe, -90.0, 90.0) == pytest.approx(face_area, rel=1e-9) + # The edge opposite the pole peaks at ~82.9 deg, so the face covers a + # quarter of the 85-90 cap. + quarter_cap = _zone_area(85.0, 90.0) / 4 + assert _band_area(fe, 85.0, 90.0) == pytest.approx(quarter_cap, rel=1e-12) + parts = _band_area(fe, 70.0, 82.0) + _band_area(fe, 82.0, 90.0) + assert parts == pytest.approx(face_area, rel=1e-9) + + def test_face_touching_band_has_zero_overlap(self): + """A face sharing only an edge with the band boundary has zero overlap.""" + fe = _face_edges(np.array([0.0, 0.0, 3.0, 3.0]), np.array([0.0, 3.0, 3.0, 0.0])) + assert _band_area(fe, -10.0, 0.0) == 0.0 + assert _band_area(fe, 0.0, 10.0) > 0.0 diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py index 77ca4bc90..62bb5b068 100644 --- a/uxarray/core/zonal.py +++ b/uxarray/core/zonal.py @@ -1,17 +1,17 @@ import numpy as np from numba import njit +from uxarray.constants import ERROR_TOLERANCE from uxarray.errors import DimensionError -from uxarray.grid.area import calculate_face_area -from uxarray.grid.geometry import _unique_points from uxarray.grid.integrate import _zonal_face_weights, _zonal_face_weights_robust -from uxarray.grid.intersections import ( - gca_const_lat_intersection, - get_number_of_intersections, -) -from uxarray.grid.utils import ( - _get_cartesian_face_edge_nodes_array_subset, - _small_angle_of_2_vectors, +from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array_subset +from uxarray.utils.numba_math import ( + _numba_add3, + _numba_cross3, + _numba_div3_scalar, + _numba_dot3, + _numba_mul3_scalar, + _numba_norm3, ) @@ -101,139 +101,201 @@ def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=Fal @njit(cache=True) -def _sort_points_by_angle(points): - """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 - Parameters - ---------- - points : ndarray - Array of 3D points on unit sphere, shape (n_points, 3) - Returns - ------- - ndarray - Points sorted by longitude angle +@njit(cache=True) +def _lon_delta(a, b): + """Signed longitude change from ``a`` to ``b`` in (-pi, pi]. + + Computed from the horizontal projections of the two points, so it is + unaffected by the antimeridian. """ - n_points = points.shape[0] - if n_points <= 1: - return points.copy() - - # Calculate angles (longitude) - angles = np.empty(n_points, dtype=np.float64) - x_axis = np.array([1.0, 0.0, 0.0]) - for i in range(n_points): - # Project point to xy plane for longitude calculation - point_xy = np.array([points[i, 0], points[i, 1], 0.0]) - point_xy_norm = np.linalg.norm(point_xy) - - if point_xy_norm < 1e-15: - angles[i] = 0.0 # Point at pole - else: - point_xy_unit = point_xy / point_xy_norm - angle = _small_angle_of_2_vectors(x_axis, point_xy_unit) - # Determine sign based on y coordinate - if points[i, 1] < 0: - angle = -angle - angles[i] = angle + return np.arctan2(a[0] * b[1] - a[1] * b[0], a[0] * b[0] + a[1] * b[1]) - # Simple insertion sort (numba-friendly for small arrays) - sorted_points = points.copy() - for i in range(1, n_points): - key_angle = angles[i] - key_point = sorted_points[i].copy() - j = i - 1 - 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. + + With ``n = a x b`` and ``t = n x p`` the (unnormalized) tangent at ``p``, an + antiderivative along the arc is ``-arctan(t_z / n_z)``. The difference of the + two arctangents is evaluated with a single ``arctan2`` so meridian arcs + (``n_z == 0``) give 0 instead of dividing by zero. + """ + n = _numba_cross3(a, b) + ta_z = n[0] * a[1] - n[1] * a[0] + tb_z = n[0] * b[1] - n[1] * b[0] + return -np.arctan2(n[2] * (tb_z - ta_z), n[2] * n[2] + ta_z * tb_z) + - angles[j + 1] = key_angle - sorted_points[j + 1] = key_point +@njit(cache=True) +def _edge_band_integral(a, b, z_min, z_max): + """Line integral of ``(clip(z, z_min, z_max) - z_min) dlon`` along the arc a -> b. - return sorted_points + The arc is split where it crosses ``z = z_min`` and ``z = z_max`` (at most + twice each), so on every piece the integrand is either a constant or ``z`` + itself, both of which integrate exactly. + + Returns the integral and whether any piece of the arc lies strictly inside + the band. + """ + a = (a[0], a[1], a[2]) + b = (b[0], b[1], b[2]) + n = _numba_cross3(a, b) + sin_arc = _numba_norm3(n) + if sin_arc == 0.0: + return 0.0, False + arc = np.arctan2(sin_arc, _numba_dot3(a, b)) + + # Parametrize the arc as p(t) = a cos(t) + u sin(t), t in [0, arc], where u is + # the unit tangent at a pointing toward b. Then z(t) = r cos(t - t0). + u = _numba_div3_scalar(_numba_cross3(n, a), sin_arc) + r = np.hypot(a[2], u[2]) + t0 = np.arctan2(u[2], a[2]) + + ts = np.empty(6, dtype=np.float64) + ts[0] = 0.0 + nt = 1 + for c in (z_min, z_max): + if np.abs(c) < r: + d = np.arccos(c / r) + for t in (t0 - d, t0 + d): + if t > np.pi: + t -= 2.0 * np.pi + elif t <= -np.pi: + t += 2.0 * np.pi + if 0.0 < t < arc: + ts[nt] = t + nt += 1 + ts[nt] = arc + nt += 1 + ts[:nt].sort() + + total = 0.0 + inside = False + p0 = a + for i in range(1, nt): + if i == nt - 1: + p1 = b + else: + p1 = _numba_add3( + _numba_mul3_scalar(a, np.cos(ts[i])), + _numba_mul3_scalar(u, np.sin(ts[i])), + ) + + # Each piece lies entirely below, inside, or above the band; classify it + # by its midpoint. The integrand is continuous in z, so a piece grazing a + # band edge contributes (almost) the same under either classification. + t_mid = 0.5 * (ts[i - 1] + ts[i]) + z_mid = a[2] * np.cos(t_mid) + u[2] * np.sin(t_mid) + if z_mid > z_min: + dlon = _lon_delta(p0, p1) + if z_mid >= z_max: + total += (z_max - z_min) * dlon + else: + total += _gca_z_dlon_integral(p0, p1) - z_min * dlon + inside = True + p0 = p1 + + return total, inside @njit(cache=True) -def _compute_band_overlap_area( - face_edges_xyz, z_min, z_max, quadrature_rule="gaussian", order=4 -): - """Compute overlap area between a face and latitude band using area.py functions. +def _compute_band_overlap_area(face_edges_xyz, z_min, z_max): + """Compute the exact overlap area between a face and a latitude band. - This function finds the intersection polygon between a face and latitude band, - then uses calculate_face_area from area.py with latitude_adjusted_area=True - for accurate area computation when edges lie on constant latitude lines. + The Lambert cylindrical equal-area map ``(lon, z)``, with ``z = sin(lat)``, + preserves area and sends the band to the strip ``z_min <= z <= z_max``. By + Green's theorem the area of the part of a face inside that strip is + + -(closed line integral of (clip(z, z_min, z_max) - z_min) dlon) + + taken around the face boundary, which is evaluated exactly on each + great-circle edge (see :func:`_edge_band_integral`). No intersection + polygon is constructed, so there is no vertex ordering to get wrong, and + longitude only enters through local differences, so faces spanning the + antimeridian need no special handling. + + Pole handling: at a vertex lying on a pole the boundary jumps between + meridians, which contributes the integrand at the pole times that jump. A + face enclosing a pole winds once around it in longitude; the strip region + above (north) or below (south) its boundary then contributes + ``2*pi*(z_max - z_min)`` or ``0`` respectively. Parameters ---------- face_edges_xyz : ndarray - Cartesian coordinates of face edge nodes, shape (n_edges, 2, 3) + Cartesian coordinates of the face's edge nodes on the unit sphere, shape + (n_edges, 2, 3), with the edges in boundary order. z_min, z_max : float - Z-coordinate bounds of the latitude band (z = sin(latitude)) - quadrature_rule : str, optional - Quadrature rule to use ("gaussian" or "triangular", default: "gaussian") - order : int, optional - Quadrature order (default: 4, same as area.py) + Z-coordinate bounds of the latitude band (z = sin(latitude)), + ``z_min <= z_max``. Returns ------- float - Overlap area between face and latitude band + Overlap area between the face and the latitude band. """ - # Pre-allocate for maximum possible intersection points - # Worst case: 2 intersections per edge * 2 boundaries + all vertices - max_points = face_edges_xyz.shape[0] * 4 + face_edges_xyz.shape[0] - polygon_points = np.empty((max_points, 3), dtype=np.float64) - point_count = 0 - - # Find intersections with z_min and z_max boundaries - z_boundaries = np.array([z_min, z_max]) - for z_boundary in z_boundaries: - for e in range(face_edges_xyz.shape[0]): - edge = face_edges_xyz[e] - inter = gca_const_lat_intersection(edge, z_boundary) - nint = get_number_of_intersections(inter) - - for i in range(nint): - if point_count < max_points: - polygon_points[point_count] = inter[i] - point_count += 1 - - # Add face vertices that lie within the band - for e in range(face_edges_xyz.shape[0]): - vertex = face_edges_xyz[e, 0] # First point of edge - if z_min <= vertex[2] <= z_max: - if point_count < max_points: - polygon_points[point_count] = vertex - point_count += 1 - - if point_count < 3: - return 0.0 - - # Remove duplicate points - unique_points = _unique_points(polygon_points[:point_count]) + n_edges = face_edges_xyz.shape[0] + line_integral = 0.0 + total_dlon = 0.0 + z_sum = 0.0 + boundary_in_band = False + + for e in range(n_edges): + a = face_edges_xyz[e, 0] + b = face_edges_xyz[e, 1] + z_sum += a[2] + + if _is_pole_point(a) or _is_pole_point(b): + # Meridian arc to or from a pole: longitude is constant along it, so + # it adds nothing to the integral, and z varies monotonically. + z_lo = min(a[2], b[2]) + z_hi = max(a[2], b[2]) + boundary_in_band |= max(z_lo, z_min) < min(z_hi, z_max) + if _is_pole_point(a): + # Turning at the pole from the incoming to the outgoing meridian + # sweeps longitude at constant z = a_z. + jump = _lon_delta(face_edges_xyz[e - 1, 0], b) + h_pole = min(max(a[2], z_min), z_max) - z_min + line_integral += h_pole * jump + total_dlon += jump + continue - if unique_points.shape[0] < 3: + total_dlon += _lon_delta(a, b) + integral, inside = _edge_band_integral(a, b, z_min, z_max) + line_integral += integral + boundary_in_band |= inside + + area = -line_integral + winding = int(np.round(total_dlon / (2.0 * np.pi))) + if winding == 0: + # Taking the magnitude makes the result independent of whether the + # face is ordered counterclockwise or clockwise. + area = np.abs(area) + else: + # The face encloses a pole. The winding direction combined with which + # pole is enclosed gives the orientation (counterclockwise winds + # eastward around the north pole and westward around the south pole). + north = z_sum > 0.0 + orientation = np.sign(winding) if north else -np.sign(winding) + area *= orientation + if north: + area += 2.0 * np.pi * (z_max - z_min) + + if not boundary_in_band: + # The boundary never enters the band, so the band either misses the + # face or lies entirely inside it (only possible for a face enclosing a + # pole). Return that exactly instead of the near-cancelling sum, so a + # face that merely touches the band gets a weight of exactly zero. + full_band = 2.0 * np.pi * (z_max - z_min) + if winding != 0 and area > 0.5 * full_band: + return full_band return 0.0 - - # Sort points to form a proper polygon - sorted_points = _sort_points_by_angle(unique_points) - - # Use area.py calculate_face_area with latitude adjustment - x = sorted_points[:, 0] - y = sorted_points[:, 1] - z = sorted_points[:, 2] - - area, _ = calculate_face_area( - x, - y, - z, - quadrature_rule=quadrature_rule, - order=order, - latitude_adjusted_area=True, # Key improvement: use latitude adjustment - ) - - return area + return max(area, 0.0) def _compute_face_band_weights(uxgrid, bands): From a4d7e3c6a538bbf8ff667eb6a8ef02b84c119274 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Mon, 5 Oct 2026 16:36:33 -0500 Subject: [PATCH 2/2] Fix inverted sign of the constant-latitude area correction 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 --- test/conftest.py | 4 +++- test/grid/integrate/test_zonal.py | 39 +++++++++++++++++++++++++++++++ uxarray/grid/area.py | 9 +++++-- 3 files changed, 49 insertions(+), 3 deletions(-) diff --git a/test/conftest.py b/test/conftest.py index 50ee2831b..22c680cb9 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -99,7 +99,9 @@ def mesh_constants(): 'NNODES_outRLL1deg': 64442, 'DATAVARS_outCSne30': 4, 'TRI_AREA': 0.02216612469199045, - 'CORRECTED_TRI_AREA': 0.02244844510268421, + # TRI_AREA with the 41.8 deg N top edge taken as a parallel instead of a + # great-circle arc; checked against brute-force integration. + 'CORRECTED_TRI_AREA': 0.021883797973858587, 'MESH30_AREA': 12.566, 'PSI_INTG': 12.566, 'VAR2_INTG': 12.566, diff --git a/test/grid/integrate/test_zonal.py b/test/grid/integrate/test_zonal.py index bbd69e378..80f6a0d22 100644 --- a/test/grid/integrate/test_zonal.py +++ b/test/grid/integrate/test_zonal.py @@ -405,6 +405,10 @@ def test_invalid_lat_input_raises(self): _compute_face_band_weights, ) from uxarray.grid.area import calculate_face_area # noqa: E402 +from uxarray.grid.intersections import ( # noqa: E402 + gca_const_lat_intersection, + get_number_of_intersections, +) from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array_subset # noqa: E402 @@ -518,6 +522,41 @@ def test_face_with_vertex_at_pole(self): parts = _band_area(fe, 70.0, 82.0) + _band_area(fe, 82.0, 90.0) assert parts == pytest.approx(face_area, rel=1e-9) + @pytest.mark.parametrize("hemisphere", [1.0, -1.0]) + def test_matches_latitude_adjusted_clipped_polygon(self, hemisphere): + """Agrees with clipping the face by hand and correcting its parallel edges. + + The clipped ring is (in counterclockwise order) a crossing of the 12 deg + parallel, a face vertex, two crossings of the 28 deg parallel and one more + crossing of 12 deg. Its great-circle area plus the constant-latitude + correction on the two parallel edges is the overlap area. + """ + lat = np.array([10.0, 15.0, 35.0, 30.0]) + lon = np.array([-8.0, 9.0, 7.0, -6.0]) + if hemisphere < 0: + # Mirror across the equator and reverse to stay counterclockwise. + lat, lon = -lat[::-1], lon[::-1] + fe = _face_edges(lat, lon) + + def crossing(edge, lat_deg): + pts = gca_const_lat_intersection(edge, np.sin(np.deg2rad(lat_deg))) + assert get_number_of_intersections(pts) == 1 + return pts[0] + + if hemisphere > 0: + ring = [crossing(fe[0], 12.0), fe[1, 0], crossing(fe[1], 28.0), + crossing(fe[3], 28.0), crossing(fe[3], 12.0)] + else: + ring = [crossing(fe[3], -28.0), crossing(fe[1], -28.0), fe[2, 0], + crossing(fe[2], -12.0), crossing(fe[3], -12.0)] + ring = np.array(ring) + + expected, _ = calculate_face_area( + ring[:, 0], ring[:, 1], ring[:, 2], order=6, latitude_adjusted_area=True + ) + lat0, lat1 = sorted((hemisphere * 12.0, hemisphere * 28.0)) + assert _band_area(fe, lat0, lat1) == pytest.approx(expected, rel=1e-12) + def test_face_touching_band_has_zero_overlap(self): """A face sharing only an edge with the band boundary has zero overlap.""" fe = _face_edges(np.array([0.0, 0.0, 3.0, 3.0]), np.array([0.0, 3.0, 3.0, 0.0])) diff --git a/uxarray/grid/area.py b/uxarray/grid/area.py index f50d9c820..a60f13eba 100644 --- a/uxarray/grid/area.py +++ b/uxarray/grid/area.py @@ -163,8 +163,13 @@ def _face_area_from_quadrature(x, y, z, dG, dW, is_gaussian, latitude_adjusted_a # Calculate the correction term correction = _area_correction(node1, node2) - # Check if the longitude is increasing in the northern hemisphere or decreasing in the southern hemisphere - if (z_sign > 0 and lon_diff > 0) or (z_sign < 0 and lon_diff < 0): + # The great-circle arc between two points on a parallel bulges + # poleward of it. For a counterclockwise face, a westward edge in + # the northern hemisphere (eastward in the southern) has the face + # on its equatorward side, so the bulge lies outside the face and + # the great-circle area overcounts it; the opposite direction + # undercounts it. + if (z_sign > 0 and lon_diff < 0) or (z_sign < 0 and lon_diff > 0): correction = -correction total_correction += correction