Skip to content

Fix gradient returning all NaN on SCRIP-derived grids #1533

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
rajeeja wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Fix gradient returning all NaN on SCRIP-derived grids #1533

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
353f5c6
Skip edges with missing neighbors in gradient instead of NaN-ing enti…
rajeeja Jun 26, 2026
43c3967
Merge branch 'main' into rajeeja/gradient_fix
rajeeja Jun 30, 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
10 changes: 7 additions & 3 deletions test/core/test_vector_calculus.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,13 +23,17 @@ def test_gradient_output_format(self, gridpath, datasetpath):
assert uxds['t2m'].sizes == grad_ds.sizes

def test_gradient_all_boundary_faces(self, gridpath, datasetpath):
"""Quad hexagon grid has 4 faces, each of which are on the boundary, so the expected gradients are zero for both components"""
"""Quad hexagon grid has 4 faces, all touching the boundary.

Each face still has some interior edges, so the gradient should
produce finite (small) values rather than NaN.
"""
uxds = ux.open_dataset(gridpath("ugrid", "quad-hexagon", "grid.nc"), datasetpath("ugrid", "quad-hexagon", "data.nc"))

grad = uxds['t2m'].gradient()

assert np.isnan(grad['meridional_gradient']).all()
assert np.isnan(grad['zonal_gradient']).all()
assert not np.isnan(grad['meridional_gradient']).any()
assert not np.isnan(grad['zonal_gradient']).any()


class TestGradientMPASOcean:
Expand Down
105 changes: 52 additions & 53 deletions uxarray/core/gradient.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -295,59 +295,58 @@ def _compute_gradients_on_faces(
# Parallel across faces
for face_idx in prange(n_face):
gradient = np.zeros(3)

if not _check_face_on_boundary(
face_idx,
face_node_connectivity,
node_edge_connectivity,
edge_face_connectivity,
): # check face is not on boundary
for node_idx in face_node_connectivity[
face_idx
]: # take each node on that face
if node_idx != INT_FILL_VALUE:
for edge_idx in node_edge_connectivity[
node_idx
]: # grab each edge connected to that node
if edge_idx != INT_FILL_VALUE:
if (
face_idx not in edge_face_connectivity[edge_idx]
): # check if edge connected to original face
face1_idx = edge_face_connectivity[edge_idx][0]
face2_idx = edge_face_connectivity[edge_idx][1]

face1_coords = face_coords[face1_idx]
face2_coords = face_coords[face2_idx]

# compute normal pointing outwards from face
cross = np.cross(face1_coords, face2_coords)
norm = np.linalg.norm(cross)
if (
np.dot(cross, face1_coords - face_coords[face_idx])
> 0
):
normal = cross / norm
else:
normal = -cross / norm

# compute arc length between the two faces
arc_length = _compute_arc_length(
face_lat[face1_idx],
face_lat[face2_idx],
face_lon[face1_idx],
face_lon[face2_idx],
)

# compute trapezoidal rule
trapz = (data[face1_idx] + data[face2_idx]) / 2

# add to the gradient (subtract correction term)
gradient = (
gradient
+ (trapz - data[face_idx]) * arc_length * normal
)

else:
has_contribution = False

for node_idx in face_node_connectivity[face_idx]: # take each node on that face
if node_idx != INT_FILL_VALUE:
for edge_idx in node_edge_connectivity[
node_idx
]: # grab each edge connected to that node
if edge_idx != INT_FILL_VALUE:
# Skip edges that lack a second face neighbor
# instead of NaN-ing the entire face. Fixes
# grids where edge_face_connectivity has
# spurious INT_FILL_VALUE entries (e.g. SCRIP-
# derived SE grids like ne120np4). See #1452.
if INT_FILL_VALUE in edge_face_connectivity[edge_idx]:
continue

if (
face_idx not in edge_face_connectivity[edge_idx]
): # check if edge connected to original face
face1_idx = edge_face_connectivity[edge_idx][0]
face2_idx = edge_face_connectivity[edge_idx][1]

face1_coords = face_coords[face1_idx]
face2_coords = face_coords[face2_idx]

# compute normal pointing outwards from face
cross = np.cross(face1_coords, face2_coords)
norm = np.linalg.norm(cross)
if np.dot(cross, face1_coords - face_coords[face_idx]) > 0:
normal = cross / norm
else:
normal = -cross / norm

# compute arc length between the two faces
arc_length = _compute_arc_length(
face_lat[face1_idx],
face_lat[face2_idx],
face_lon[face1_idx],
face_lon[face2_idx],
)

# compute trapezoidal rule
trapz = (data[face1_idx] + data[face2_idx]) / 2

# add to the gradient (subtract correction term)
gradient = (
gradient
+ (trapz - data[face_idx]) * arc_length * normal
)
has_contribution = True

if not has_contribution:
gradient = np.full(3, np.nan)

node_neighbors = face_node_connectivity[face_idx]
Expand Down
Loading