Repository: UXARRAY/uxarray

Issue:
Add Standardized Data Types and Fill Values for UGRID connectivity variables
Use a standardized data type and fill values for connectivity variables. Continutation of #189 

Discussion:

Patch:
diff --git a/uxarray/_exodus.py b/uxarray/_exodus.py
index f3070622..d6af036e 100644
--- a/uxarray/_exodus.py
+++ b/uxarray/_exodus.py
@@ -3,7 +3,8 @@ import numpy as np
 from pathlib import PurePath
 from datetime import datetime
 
-int_dtype = np.uint32
+from uxarray.helpers import _replace_fill_values
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 
 # Exodus Number is one-based.
@@ -47,7 +48,7 @@ def _read_exodus(ext_ds, ds_var_names):
             # TODO: Use the data here for Mesh2 construct, if required.
             pass
         elif key == "coord":
-            ds.Mesh2.attrs['topology_dimension'] = int_dtype(
+            ds.Mesh2.attrs['topology_dimension'] = INT_DTYPE(
                 ext_ds.dims['num_dim'])
             ds["Mesh2_node_x"] = xr.DataArray(
                 data=ext_ds.coord[0],
@@ -120,16 +121,23 @@ def _read_exodus(ext_ds, ds_var_names):
 
     # outside the k,v for loop
     # set the face nodes data compiled in "connect" section
+
+    # standardize fill values and data type face nodes
+    face_nodes = _replace_fill_values(grid_var=conn[:] - 1,
+                                      original_fill=-1,
+                                      new_fill=INT_FILL_VALUE,
+                                      new_dtype=INT_DTYPE)
+
     ds["Mesh2_face_nodes"] = xr.DataArray(
-        data=(conn[:] - 1),
+        data=face_nodes,
         dims=["nMesh2_face", "nMaxMesh2_face_nodes"],
         attrs={
             "cf_role":
                 "face_node_connectivity",
             "_FillValue":
-                -1,
+                INT_FILL_VALUE,
             "start_index":
-                int_dtype(
+                INT_DTYPE(
                     0)  # NOTE: This might cause an error if numbering has holes
         })
     print("Finished reading exodus file.")
@@ -168,7 +176,7 @@ def _encode_exodus(ds, ds_var_names, outfile=None):
     now = datetime.now()
     date = now.strftime("%Y:%m:%d")
     time = now.strftime("%H:%M:%S")
-    fp_word = int_dtype(8)
+    fp_word = INT_DTYPE(8)
     exo_version = np.float32(5.0)
     api_version = np.float32(5.0)
 
@@ -228,7 +236,7 @@ def _encode_exodus(ds, ds_var_names, outfile=None):
     conn_nofill = []
 
     # store the number of faces in an array
-    for row in ds[ds_var_names["Mesh2_face_nodes"]].astype(int_dtype).data:
+    for row in ds[ds_var_names["Mesh2_face_nodes"]].astype(INT_DTYPE).data:
 
         # find out -1 in each row, this indicates lower than max face nodes
         arr = np.where(row == -1)
@@ -283,14 +291,15 @@ def _encode_exodus(ds, ds_var_names, outfile=None):
         # assign Data variables
         # convert list to np.array, sorted list guarantees we have the correct info
         conn_blk = conn_nofill[start:start + num_faces]
-        conn_np = np.array([np.array(xi, dtype="i8") for xi in conn_blk])
+        conn_np = np.array([np.array(xi, dtype=INT_DTYPE) for xi in conn_blk])
         exo_ds[str_connect] = xr.DataArray(data=xr.DataArray((conn_np[:] + 1)),
                                            dims=[str_el_in_blk, str_nod_per_el],
                                            attrs={"elem_type": element_type})
 
         # edge type
         exo_ds[str_edge_type] = xr.DataArray(
-            data=xr.DataArray(np.zeros((num_faces, num_nodes), "i8")),
+            data=xr.DataArray(np.zeros((num_faces, num_nodes),
+                                       dtype=INT_DTYPE)),
             dims=[str_el_in_blk, str_nod_per_el])
 
         # global id
@@ -315,7 +324,7 @@ def _encode_exodus(ds, ds_var_names, outfile=None):
                                       attrs={"name": "ID"})
     # eb_status
     exo_ds["eb_status"] = xr.DataArray(data=xr.DataArray(
-        np.ones([num_blks], dtype="i8")),
+        np.ones([num_blks], dtype=INT_DTYPE)),
                                        dims=["num_el_blk"])
 
     # eb_names
diff --git a/uxarray/_mpas.py b/uxarray/_mpas.py
index 0e2af012..5aa8436b 100644
--- a/uxarray/_mpas.py
+++ b/uxarray/_mpas.py
@@ -2,9 +2,7 @@ import xarray as xr
 import numpy as np
 import warnings
 
-# edit once PR #241 is merged
-INT_DTYPE = np.uint32
-INT_FILL_VALUE = np.iinfo(INT_DTYPE).max
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 
 def _primal_to_ugrid(in_ds, out_ds):
diff --git a/uxarray/_scrip.py b/uxarray/_scrip.py
index 3464fd14..7ab4b3bf 100644
--- a/uxarray/_scrip.py
+++ b/uxarray/_scrip.py
@@ -1,9 +1,10 @@
 import xarray as xr
 import numpy as np
 
-from .helpers import grid_center_lat_lon
+from uxarray.helpers import grid_center_lat_lon
 
-int_dtype = np.uint32
+from uxarray.helpers import _replace_fill_values
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 
 def _to_ugrid(in_ds, out_ds):
@@ -70,17 +71,23 @@ def _to_ugrid(in_ds, out_ds):
         out_ds['Mesh2_face_x'] = in_ds['grid_center_lon']
         out_ds['Mesh2_face_y'] = in_ds['grid_center_lat']
 
+        # standardize fill values and data type face nodes
+        face_nodes = _replace_fill_values(unq_inv,
+                                          original_fill=-1,
+                                          new_fill=INT_FILL_VALUE,
+                                          new_dtype=INT_DTYPE)
+
         # set the face nodes data compiled in "connect" section
         out_ds["Mesh2_face_nodes"] = xr.DataArray(
-            data=unq_inv,
+            data=face_nodes,
             dims=["nMesh2_face", "nMaxMesh2_face_nodes"],
             attrs={
                 "cf_role":
                     "face_node_connectivity",
                 "_FillValue":
-                    -1,
+                    INT_FILL_VALUE,
                 "start_index":
-                    int_dtype(
+                    INT_DTYPE(
                         0
                     )  # NOTE: This might cause an error if numbering has holes
             })
@@ -184,7 +191,7 @@ def _encode_scrip(mesh2_face_nodes, mesh2_node_x, mesh2_node_y, face_areas):
     ds = xr.Dataset()
 
     # Make grid corner lat/lon
-    f_nodes = mesh2_face_nodes.values.astype(int_dtype).ravel()
+    f_nodes = mesh2_face_nodes.values.astype(INT_DTYPE).ravel()
 
     # Create arrays to hold lat/lon data
     lat_nodes = mesh2_node_y[f_nodes].values
diff --git a/uxarray/_ugrid.py b/uxarray/_ugrid.py
index 55b41ae3..45f0e807 100644
--- a/uxarray/_ugrid.py
+++ b/uxarray/_ugrid.py
@@ -1,5 +1,8 @@
 import xarray as xr
 
+from uxarray.helpers import _replace_fill_values
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
+
 
 def _read_ugrid(xr_ds, var_names_dict):
     """UGRID file reader.
@@ -49,11 +52,15 @@ def _read_ugrid(xr_ds, var_names_dict):
             var_names_dict["Mesh2_node_x"], var_names_dict["Mesh2_node_y"],
             var_names_dict["Mesh2_node_z"]
         ])
+
+    # standardize fill values and data type for face nodes
+    xr_ds = _standardize_fill_values(xr_ds, var_names_dict)
+
     return xr_ds, var_names_dict
 
 
 def _encode_ugrid(ds):
-    """Encodes UGRID file.
+    """Encodes UGRID file .
     Parameters
     ----------
     ds : xarray.Dataset
@@ -62,3 +69,45 @@ def _encode_ugrid(ds):
     Uses to_netcdf from xarray object.
     """
     return ds
+
+
+def _standardize_fill_values(ds, var_names_dict):
+    """Standardizes the fill values and data type of index variables.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        Input Dataset
+    ds_vars_dict : dict
+        UGRID Variable Name Dictionary
+
+    Returns
+    ----------
+    ds : xarray.Dataset
+        Input Dataset with correct index variables
+    """
+
+    # original face nodes
+    face_nodes = ds[var_names_dict['Mesh2_face_nodes']].values
+
+    # original fill value, if one exists
+    if "_FillValue" in ds[var_names_dict['Mesh2_face_nodes']].attrs:
+        current_fv = ds[var_names_dict['Mesh2_face_nodes']]._FillValue
+    else:
+        current_fv = None
+
+    # if current dtype and fill value are not standardized
+    if face_nodes.dtype != INT_DTYPE or current_fv != INT_FILL_VALUE:
+        # replace fill values and set correct dtype
+        new_face_nodes = _replace_fill_values(grid_var=face_nodes,
+                                              original_fill=current_fv,
+                                              new_fill=INT_FILL_VALUE,
+                                              new_dtype=INT_DTYPE)
+        # reassign data to use updated face nodes
+        ds[var_names_dict['Mesh2_face_nodes']].data = new_face_nodes
+
+        # use new fill value
+        ds[var_names_dict['Mesh2_face_nodes']].attrs[
+            '_FillValue'] = INT_FILL_VALUE
+
+    return ds
diff --git a/uxarray/constants.py b/uxarray/constants.py
new file mode 100644
index 00000000..b3156327
--- /dev/null
+++ b/uxarray/constants.py
@@ -0,0 +1,4 @@
+import numpy as np
+
+INT_DTYPE = np.uintp
+INT_FILL_VALUE = np.iinfo(INT_DTYPE).max
diff --git a/uxarray/grid.py b/uxarray/grid.py
index c1abd6a8..8e31114f 100644
--- a/uxarray/grid.py
+++ b/uxarray/grid.py
@@ -11,8 +11,6 @@ from ._scrip import _read_scrip, _encode_scrip
 from ._mpas import _read_mpas
 from .helpers import get_all_face_area_from_coords, parse_grid_type, _convert_node_xyz_to_lonlat_rad, _convert_node_lonlat_rad_to_xyz
 
-int_dtype = np.uint32
-
 
 class Grid:
     """
@@ -309,7 +307,7 @@ class Grid:
             if not "degree" in self.Mesh2_node_x.units:
                 coords_type = "cartesian"
 
-            face_nodes = self.Mesh2_face_nodes.data.astype(int_dtype)
+            face_nodes = self.Mesh2_face_nodes.data
             dim = self.Mesh2.attrs['topology_dimension']
 
             # initialize z
diff --git a/uxarray/helpers.py b/uxarray/helpers.py
index e8e21476..9d311f79 100644
--- a/uxarray/helpers.py
+++ b/uxarray/helpers.py
@@ -5,8 +5,9 @@ from .get_quadratureDG import get_gauss_quadratureDG, get_tri_quadratureDG
 from numba import njit, config
 import math
 
+from .constants import INT_DTYPE
+
 config.DISABLE_JIT = False
-int_dtype = np.uint32
 
 
 def parse_grid_type(dataset):
@@ -228,7 +229,7 @@ def get_all_face_area_from_coords(x,
     num_faces = face_nodes.shape[0]
     area = np.zeros(num_faces)  # set area of each face to 0
 
-    face_nodes = face_nodes[:].astype(int_dtype)
+    face_nodes = face_nodes[:].astype(INT_DTYPE)
 
     for i in range(num_faces):
         face_z = np.zeros(len(face_nodes[i]))
@@ -558,3 +559,45 @@ def _normalize_in_place(node):
         raise RuntimeError("Input array should have a length of 3: [x, y, z]")
 
     return list(np.array(node) / np.linalg.norm(np.array(node), ord=2))
+
+
+def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
+    """Replaces all instances of the the current fill value (``original_fill``)
+    in (``grid_var``) with (``new_fill``) and converts to the dtype defined by
+    (``new_dtype``)
+
+    Parameters
+    ----------
+    grid_var : np.ndarray
+        grid variable to be modified
+    original_fill : constant
+        original fill value used in (``grid_var``)
+    new_fill : constant
+        new fill value to be used in (``grid_var``)
+    new_dtype : np.dtype
+        new data type to convert (``grid_var``) to
+
+    Returns
+    ----------
+    grid_var : xarray.Dataset
+        Input Dataset with correct fill value and dtype
+    """
+
+    # locations of fill values
+    fill_val_idx = grid_var.ravel() == original_fill
+
+    # convert to new data type
+    if new_dtype != grid_var.dtype and new_dtype is not None:
+        grid_var = grid_var.astype(new_dtype)
+
+    # ensure select fill value can be represented with current dtype
+    dtype_min = np.iinfo(grid_var.dtype).min
+    dtype_max = np.iinfo(grid_var.dtype).max
+    if new_fill not in range(dtype_min, dtype_max + 1):
+        raise ValueError(f'New fill value: {new_fill} not representable by'
+                         f' dtype: {grid_var.dtype}')
+
+    # replace all zeros with a fill value
+    grid_var.ravel()[fill_val_idx] = new_fill
+
+    return grid_var


Test patch:
diff --git a/test/constants.py b/test/constants.py
index f505274f..a2749ad7 100644
--- a/test/constants.py
+++ b/test/constants.py
@@ -11,4 +11,3 @@ TRI_AREA = 1.047
 MESH30_AREA = 12.566
 PSI_INTG = 12.566
 VAR2_INTG = 12.566
-int_dtype = np.uint32
diff --git a/test/test_exodus.py b/test/test_exodus.py
index 6137c7af..68375ab9 100644
--- a/test/test_exodus.py
+++ b/test/test_exodus.py
@@ -6,6 +6,7 @@ from pathlib import Path
 
 import xarray as xr
 import uxarray as ux
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
@@ -39,3 +40,14 @@ class TestExodus(TestCase):
         tgrid.encode_as("ugrid")
         outfile = current_path / "write_test_mixed.exo"
         tgrid.encode_as("exodus")
+
+    def test_standardized_dtype_and_fill(self):
+        """Test to see if Mesh2_Face_Nodes uses the expected integer datatype
+        and expected fill value as set in constants.py."""
+
+        exo2_filename = current_path / "meshfiles" / "exodus" / "mixed" / "mixed.exo"
+        xr_exo_ds = xr.open_dataset(exo2_filename)
+        ux_grid = ux.Grid(xr_exo_ds)
+
+        assert ux_grid.Mesh2_face_nodes.dtype == INT_DTYPE
+        assert ux_grid.Mesh2_face_nodes._FillValue == INT_FILL_VALUE
diff --git a/test/test_helpers.py b/test/test_helpers.py
index ad5d43e9..72c494de 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -9,6 +9,9 @@ from pathlib import Path
 
 import uxarray as ux
 
+from uxarray.helpers import _replace_fill_values
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
+
 try:
     import constants
 except ImportError:
@@ -30,7 +33,7 @@ class TestIntegrate(TestCase):
         x = np.array([0.57735027, 0.57735027, -0.57735027])
         y = np.array([-5.77350269e-01, 5.77350269e-01, 5.77350269e-01])
         z = np.array([-0.57735027, -0.57735027, -0.57735027])
-        face_nodes = np.array([[0, 1, 2]]).astype(constants.int_dtype)
+        face_nodes = np.array([[0, 1, 2]]).astype(INT_DTYPE)
         area = ux.get_all_face_area_from_coords(x,
                                                 y,
                                                 z,
@@ -119,3 +122,46 @@ class TestCoordinatesConversion(TestCase):
          new_lat] = ux.helpers._convert_node_xyz_to_lonlat_rad([x, y, z])
         self.assertLessEqual(np.absolute(new_lon - lon), err_tolerance)
         self.assertLessEqual(np.absolute(new_lat - lat), err_tolerance)
+
+
+class TestConstants(TestCase):
+    # DTYPE as set in constants.py
+    expected_int_dtype = INT_DTYPE
+
+    # INT_FILL_VALUE as set in constants.py
+    fv = INT_FILL_VALUE
+
+    def test_replace_fill_values(self):
+        """Tests _replace_fill_values() helper function across multiple
+        different dtype arrays used as face_nodes."""
+
+        # expected output from _replace_fill_values()
+        face_nodes_gold = np.array(
+            [[1, 2, self.fv], [self.fv, self.fv, self.fv]], dtype=INT_DTYPE)
+
+        # test different datatypes for face_nodes
+        dtypes = [np.int32, np.int64, np.float32, np.float64]
+        for dtype in dtypes:
+            # test face nodes with set dtype
+            face_nodes = np.array([[1, 2, -1], [-1, -1, -1]], dtype=dtype)
+
+            # output of _replace_fill_values()
+            face_nodes_test = _replace_fill_values(grid_var=face_nodes,
+                                                   original_fill=-1,
+                                                   new_fill=INT_FILL_VALUE,
+                                                   new_dtype=INT_DTYPE)
+
+            assert np.array_equal(face_nodes_test, face_nodes_gold)
+
+    def test_replace_fill_values_invalid(self):
+        """Tests _replace_fill_values() helper function attempting to use a
+        fill value that is not representable by the current dtype."""
+
+        face_nodes = np.array([[1, 2, -1], [-1, -1, -1]], dtype=np.uint32)
+        # invalid fill value with dtype should raise a valueError
+        with self.assertRaises(ValueError):
+            # INT_FILL_VALUE (max(uint32) not representable by int16)
+            face_nodes_test = _replace_fill_values(grid_var=face_nodes,
+                                                   original_fill=-1,
+                                                   new_fill=INT_FILL_VALUE,
+                                                   new_dtype=np.int16)
diff --git a/test/test_mpas.py b/test/test_mpas.py
index 7741fa8c..f5b1d46b 100644
--- a/test/test_mpas.py
+++ b/test/test_mpas.py
@@ -7,6 +7,8 @@ import numpy as np
 import os
 from pathlib import Path
 
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
+
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
 
@@ -18,8 +20,7 @@ class TestMPAS(TestCase):
     mpas_xr_ds = xr.open_dataset(mpas_grid_path)
 
     # fill value (remove once there is a unified approach in uxarray)
-    int_dtype = np.uint32
-    fv = np.iinfo(int_dtype).max
+    fv = INT_FILL_VALUE
 
     def test_read_mpas(self):
         """Tests execution of _read_mpas()"""
@@ -79,7 +80,7 @@ class TestMPAS(TestCase):
 
         # two cells with 2, 3 and 2 padded faces respectively
         verticesOnCell = np.array([[1, 2, 1, 1], [3, 4, 5, 3], [6, 7, 0, 0]],
-                                  dtype=self.int_dtype)
+                                  dtype=INT_DTYPE)
 
         # cell has 2, 3 and 2 nodes respectively
         nEdgesOnCell = np.array([2, 3, 2])
@@ -87,7 +88,7 @@ class TestMPAS(TestCase):
         # expected output of _add_fill_values()
         gold_output = np.array([[0, 1, self.fv, self.fv], [2, 3, 4, self.fv],
                                 [5, 6, self.fv, self.fv]],
-                               dtype=self.int_dtype)
+                               dtype=INT_DTYPE)
 
         # test data output
         verticesOnCell = _replace_padding(verticesOnCell, nEdgesOnCell)
diff --git a/test/test_scrip.py b/test/test_scrip.py
index 74c6e404..34da85a3 100644
--- a/test/test_scrip.py
+++ b/test/test_scrip.py
@@ -1,4 +1,5 @@
 from uxarray._scrip import _read_scrip, _encode_scrip
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 import uxarray as ux
 import xarray as xr
 from unittest import TestCase
@@ -110,3 +111,18 @@ class TestScrip(TestCase):
         strip_lon = np.unique(corner_lon)
 
         assert strip_lon.all() == mesh08['Mesh2_node_x'].all()
+
+    def test_standardized_dtype_and_fill(self):
+        """Test to see if Mesh2_Face_Nodes uses the expected integer datatype
+        and expected fill value as set in constants.py."""
+
+        xr_ne30 = xr.open_dataset(ne30)
+        ux_grid_01 = ux.Grid(xr_ne30)
+
+        xr_ne8 = xr.open_dataset(ne8)
+        ux_grid_02 = ux.Grid(xr_ne8)
+
+        grids = [ux_grid_01, ux_grid_02]
+        for grid in grids:
+            assert grid.Mesh2_face_nodes.dtype == INT_DTYPE
+            assert grid.Mesh2_face_nodes._FillValue == INT_FILL_VALUE
diff --git a/test/test_ugrid.py b/test/test_ugrid.py
index cb94cc42..ebdedb65 100644
--- a/test/test_ugrid.py
+++ b/test/test_ugrid.py
@@ -6,6 +6,7 @@ from pathlib import Path
 import warnings
 
 import uxarray as ux
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 try:
     import constants
@@ -63,3 +64,38 @@ class TestUgrid(TestCase):
         xr_grid = xr.open_dataset(str(exo2_filename))
         ux_grid = ux.Grid(xr_grid)
         ux_grid.encode_as("ugrid")
+
+    def test_standardized_dtype_and_fill(self):
+        """Test to see if Mesh2_Face_Nodes uses the expected integer datatype
+        and expected fill value as set in constants.py."""
+
+        ug_filename1 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
+        ug_filename2 = current_path / "meshfiles" / "ugrid" / "outRLL1deg" / "outRLL1deg.ug"
+        ug_filename3 = current_path / "meshfiles" / "ugrid" / "ov_RLL10deg_CSne4" / "ov_RLL10deg_CSne4.ug"
+
+        xr_grid1 = xr.open_dataset(str(ug_filename1))
+        xr_grid2 = xr.open_dataset(str(ug_filename2))
+        xr_grid3 = xr.open_dataset(str(ug_filename3))
+
+        ux_grid1 = ux.Grid(xr_grid1)
+        ux_grid2 = ux.Grid(xr_grid2)
+        ux_grid3 = ux.Grid(xr_grid3)
+
+        # check for correct dtype and fill value
+        grids = [ux_grid1, ux_grid2, ux_grid3]
+        for grid in grids:
+            assert grid.Mesh2_face_nodes.dtype == INT_DTYPE
+            assert grid.Mesh2_face_nodes._FillValue == INT_FILL_VALUE
+
+    def test_standardized_dtype_and_fill_dask(self):
+        """Test to see if Mesh2_Face_Nodes uses the expected integer datatype
+        and expected fill value as set in constants.py.
+
+        with dask chunking
+        """
+        ug_filename = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
+        xr_grid = xr.open_dataset(str(ug_filename), chunks={'nMesh2_node': 100})
+        ux_grid = ux.Grid(xr_grid)
+        assert ux_grid.Mesh2_face_nodes.dtype == INT_DTYPE
+        assert ux_grid.Mesh2_face_nodes._FillValue == INT_FILL_VALUE
+        pass