Fix edge_face_distance functions in neighbors.py. Add unit test in test_neighbors.py (#1293)

* changed from using nodes to face values

* added a test for construct_edge_face_distances

* adjusted to np arrays in test_neighbors.py

* fixed expected values in test_quad_hex in test_gradient

---------

Co-authored-by: Philip Chmielowiec <[email protected]>

Author: egallmeier

test/test_gradient.py

@@ -53,7 +53,7 @@ def test_quad_hex():
         else:
             assert grad.values[i] != 0
 
-    expected_values = np.array([27.95, 20.79, 28.96, 0, 0, 0, 0, 60.64, 0, 86.45, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+    expected_values = np.array([28.963, 13.100, 14.296, 0, 0, 0, 0, 67.350, 0, 85.9397, 0, 0, 0, 0, 0, 0, 0, 0, 0])
     nt.assert_almost_equal(grad.values, expected_values, 1e-2)
 
 def test_normalization():

test/test_neighbors.py

@@ -4,6 +4,7 @@
 import xarray as xr
 from pathlib import Path
 import uxarray as ux
+from uxarray.grid.neighbors import _construct_edge_face_distances
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
@@ -157,3 +158,30 @@ def test_kdtree_multi_point_query():
         for i, cur_c in enumerate(c):
             single_ind = uxgrid.get_kd_tree(coordinates="nodes").query_radius(cur_c, 45)
             assert np.array_equal(single_ind, multi_ind[i])
+
+def test_construct_edge_face_distances():
+    """
+    Test _construct_edge_face_distances by verifying known great-circle distances
+    between face centers on a unit sphere.
+    """
+    face_lon = np.array([0, 0, 90, 90, -45])
+    face_lat = np.array([0, 90, 0, 90, 0])
+    edge_faces = np.array([
+            [0, 1],  # from (0,0) to (0,90)
+            [0, 2],  # from (0,0) to (90,0)
+            [1, 3],  # from (0,90) to (90,90) — both poles, same point
+            [2, 4],  # from (90,0) to (-45,0)
+        ])
+
+
+    # Expected great-circle distances in radians
+    expected = np.array([
+        np.pi / 2,  # 0 → 1
+        np.pi / 2,  # 0 → 2
+        0.0,  # 1 → 3
+        3 * np.pi / 4  # 2 → 4
+    ])
+
+    # Run the function under test
+    calculated = _construct_edge_face_distances(face_lon, face_lat, edge_faces)
+    np.testing.assert_array_almost_equal(calculated, expected, decimal=5)

uxarray/grid/neighbors.py

@@ -1097,7 +1097,7 @@ def _construct_edge_node_distances(node_lon, node_lat, edge_nodes):
 def _populate_edge_face_distances(grid):
     """Populates ``edge_face_distances``"""
     edge_face_distances = _construct_edge_face_distances(
-        grid.node_lon.values, grid.node_lat.values, grid.edge_face_connectivity.values
+        grid.face_lon.values, grid.face_lat.values, grid.edge_face_connectivity.values
     )
 
     grid._ds["edge_face_distances"] = xr.DataArray(
@@ -1110,19 +1110,19 @@ def _populate_edge_face_distances(grid):
 
 
 @njit(cache=True)
-def _construct_edge_face_distances(node_lon, node_lat, edge_faces):
+def _construct_edge_face_distances(face_lon, face_lat, edge_faces):
     """Helper for computing the arc-distance between faces that saddle a given
     edge."""
 
     saddle_mask = edge_faces[:, 1] != INT_FILL_VALUE
 
     edge_face_distances = np.zeros(edge_faces.shape[0])
 
-    edge_lon_a = np.deg2rad((node_lon[edge_faces[saddle_mask, 0]]))
-    edge_lon_b = np.deg2rad((node_lon[edge_faces[saddle_mask, 1]]))
+    edge_lon_a = np.deg2rad((face_lon[edge_faces[saddle_mask, 0]]))
+    edge_lon_b = np.deg2rad((face_lon[edge_faces[saddle_mask, 1]]))
 
-    edge_lat_a = np.deg2rad((node_lat[edge_faces[saddle_mask, 0]]))
-    edge_lat_b = np.deg2rad((node_lat[edge_faces[saddle_mask, 1]]))
+    edge_lat_a = np.deg2rad((face_lat[edge_faces[saddle_mask, 0]]))
+    edge_lat_b = np.deg2rad((face_lat[edge_faces[saddle_mask, 1]]))
 
     # arc length
     edge_face_distances[saddle_mask] = np.arccos(