Update units for parsed distances from MPAS (#1179)

* update units when parsing distances from mpas and docstring

* add test, add warning to SpatialHash

* correct phrasing

* add periodic element note

* update note

* update test

Author: philipc2

test/test_mpas.py

@@ -4,6 +4,7 @@
 import uxarray as ux
 import xarray as xr
 from pathlib import Path
+import numpy.testing as nt
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 from uxarray.io._mpas import _replace_padding, _replace_zeros, _to_zero_index, _read_mpas
 
@@ -110,3 +111,14 @@ def test_face_area():
 
     assert "face_areas" in uxgrid_primal._ds
     assert "face_areas" in uxgrid_dual._ds
+
+
+def test_distance_units():
+    xrds = xr.open_dataset(mpas_ocean_mesh)
+    uxgrid = ux.open_grid(mpas_ocean_mesh)
+
+    assert "edge_node_distances" in uxgrid._ds
+    assert "edge_face_distances" in uxgrid._ds
+
+    nt.assert_array_almost_equal(uxgrid['edge_node_distances'].values, (xrds['dvEdge'].values / xrds.attrs['sphere_radius']))
+    nt.assert_array_almost_equal(uxgrid['edge_face_distances'].values, (xrds['dcEdge'].values / xrds.attrs['sphere_radius']))

uxarray/grid/grid.py

@@ -1387,7 +1387,7 @@ def node_face_connectivity(self, value):
 
     @property
     def edge_node_distances(self):
-        """Distances between the two nodes that surround each edge in degrees.
+        """Distances between the two nodes that surround each edge in radians.
 
         Dimensions ``(n_edge, )``
         """
@@ -1404,7 +1404,7 @@ def edge_node_distances(self, value):
     @property
     def edge_face_distances(self):
         """Distances between the centers of the faces that saddle each edge in
-        degrees.
+        radians.
 
         Dimensions ``(n_edge, )``
         """
@@ -1782,6 +1782,10 @@ def get_spatial_hash(
         self._spatialhash : grid.Neighbors.SpatialHash
             SpatialHash instance
 
+        Note
+        ----
+        Does not currently support queries on periodic elements.
+
         Examples
         --------
         Open a grid from a file path:

uxarray/grid/neighbors.py

@@ -792,6 +792,10 @@ class SpatialHash:
         Source grid used to construct the hash grid and hash table
     reconstruct : bool, default=False
         If true, reconstructs the spatial hash
+
+    Note
+    ----
+    Does not currently support queries on periodic elements.
     """
 
     def __init__(
@@ -806,11 +810,21 @@ def __init__(
 
         # Hash grid size
         self._dh = self._hash_cell_size()
+
         # Lower left corner of the hash grid
-        self._xmin = np.deg2rad(self._source_grid.node_lon.min().to_numpy()) - self._dh
-        self._ymin = np.deg2rad(self._source_grid.node_lat.min().to_numpy()) - self._dh
-        self._xmax = np.deg2rad(self._source_grid.node_lon.max().to_numpy()) + self._dh
-        self._ymax = np.deg2rad(self._source_grid.node_lat.max().to_numpy()) + self._dh
+        lon_min = np.deg2rad(self._source_grid.node_lon.min().to_numpy())
+        lat_min = np.deg2rad(self._source_grid.node_lat.min().to_numpy())
+        lon_max = np.deg2rad(self._source_grid.node_lon.max().to_numpy())
+        lat_max = np.deg2rad(self._source_grid.node_lat.max().to_numpy())
+
+        self._xmin = lon_min - self._dh
+        self._ymin = lat_min - self._dh
+        self._xmax = lon_max + self._dh
+        self._ymax = lat_max + self._dh
+
+        print(self._xmin, self._xmax)
+        print(self._ymin, self._ymax)
+
         # Number of x points in the hash grid; used for
         # array flattening
         Lx = self._xmax - self._xmin
@@ -866,10 +880,15 @@ def _initialize_face_hash_table(self):
             )
             i2, j2 = self._hash_index2d(coords)
 
-            for eid in range(self._source_grid.n_face):
-                for j in range(j1[eid], j2[eid] + 1):
-                    for i in range(i1[eid], i2[eid] + 1):
-                        index_to_face[i + self._nx * j].append(eid)
+            try:
+                for eid in range(self._source_grid.n_face):
+                    for j in range(j1[eid], j2[eid] + 1):
+                        for i in range(i1[eid], i2[eid] + 1):
+                            index_to_face[i + self._nx * j].append(eid)
+            except IndexError:
+                raise IndexError(
+                    "list index out of range. This may indicate incorrect `edge_node_distances` values."
+                )
 
             return index_to_face
 

uxarray/io/_mpas.py

@@ -403,7 +403,7 @@ def _parse_edge_faces(in_ds, out_ds, mesh_type):
 
 def _parse_edge_node_distances(in_ds, out_ds):
     """Parses ``edge_node_distances``"""
-    edge_node_distances = in_ds["dvEdge"]
+    edge_node_distances = in_ds["dvEdge"] / in_ds.attrs["sphere_radius"]
 
     out_ds["edge_node_distances"] = edge_node_distances.assign_attrs(
         descriptors.EDGE_NODE_DISTANCES_ATTRS
@@ -412,7 +412,7 @@ def _parse_edge_node_distances(in_ds, out_ds):
 
 def _parse_edge_face_distances(in_ds, out_ds):
     """Parses ``edge_face_distances``"""
-    edge_face_distances = in_ds["dcEdge"]
+    edge_face_distances = in_ds["dcEdge"] / in_ds.attrs["sphere_radius"]
 
     out_ds["edge_face_distances"] = edge_face_distances.assign_attrs(
         descriptors.EDGE_FACE_DISTANCES_ATTRS