Repository: UXARRAY/uxarray

Issue:
Move connectivity construction functions to separate location
Move connectivity construction functions such as `Grid._build_edge_node_connectivity()` into a separate location, like `connectivity.py` 

Discussion:

Patch:
diff --git a/docs/examples/003-area-calc.ipynb b/docs/examples/003-area-calc.ipynb
index c60cebc8..62dede7c 100644
--- a/docs/examples/003-area-calc.ipynb
+++ b/docs/examples/003-area-calc.ipynb
@@ -586,7 +586,7 @@
    },
    "outputs": [],
    "source": [
-    "from uxarray.utils.helpers import calculate_face_area"
+    "from uxarray.grid.area import calculate_face_area"
    ]
   },
   {
diff --git a/docs/examples/004-working-with-mpas-grids.ipynb b/docs/examples/004-working-with-mpas-grids.ipynb
index 2fa0d8d7..a94949c1 100644
--- a/docs/examples/004-working-with-mpas-grids.ipynb
+++ b/docs/examples/004-working-with-mpas-grids.ipynb
@@ -447,7 +447,7 @@
     "import matplotlib.pyplot as plt\n",
     "from matplotlib.collections import PolyCollection\n",
     "\n",
-    "from uxarray.utils.helpers import close_face_nodes"
+    "from uxarray.grid.connectivity import close_face_nodes"
    ]
   },
   {
diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
index 8b569296..2daffb7e 100644
--- a/docs/internal_api/index.rst
+++ b/docs/internal_api/index.rst
@@ -1,8 +1,8 @@
 .. currentmodule:: uxarray
 
-########
+############
 Internal API
-########
+############
 
 This page shows already-implemented Uxarray internal API functions. You can also
 check the draft `UXarray API
@@ -23,28 +23,27 @@ Class
 
    UxDataset
 
-
 Attributes
 ----------
 .. autosummary::
    :toctree: _autosummary
+
    UxDataset._source_datasets
    UxDataset._uxgrid
 
 
-
 Methods
 -------
 .. autosummary::
    :toctree: _autosummary
-   __getitem__
-   __setitem__
-   _calculate_binary_op
-   _construct_dataarray
-   _construct_direct
-   _copy
-   _replace
 
+   UxDataset.__getitem__
+   UxDataset.__setitem__
+   UxDataset._calculate_binary_op
+   UxDataset._construct_dataarray
+   UxDataset._construct_direct
+   UxDataset._copy
+   UxDataset._replace
 
 
 UxDataArray
@@ -73,11 +72,9 @@ Methods
 .. autosummary::
    :toctree: _autosummary
 
-   _construct_direct
-   _copy
-   _replace
-
-
+   UxDataArray._construct_direct
+   UxDataArray._copy
+   UxDataArray._replace
 
 Grid
 ===========
@@ -89,45 +86,15 @@ Class
 
    Grid
 
-
-IO
-----------
-.. autosummary::
-   :toctree: _autosummary
-
-   io._exodus._read_exodus
-   io._exodus._encode_exodus
-   io._exodus._get_element_type
-   io._mpas._dual_to_ugrid
-   io._mpas._primal_to_ugrid
-   io._mpas._replace_padding
-   io._mpas._replace_zeros
-   io._mpas._to_zero_index
-   io._mpas._set_global_attrs
-   io._mpas._read_mpas
-   io._ugrid._encode_ugrid
-   io._ugrid._read_ugrid
-   io._scrip._read_scrip
-   io._scrip._encode_scrip
-   io._scrip._to_ugrid
-
-
 Methods
 -------
 .. autosummary::
    :toctree: _autosummary
 
-   __init_grid_var_names__
-   __from_ds__
-   __from_vert__
-   __init_grid_var_attrs__
-   _build_edge_node_connectivity
-   _build_face_edges_connectivity
-   _build_nNodes_per_face
-   _populate_cartesian_xyz_coord
-   _populate_lonlat_coord
-
-
+   Grid.__init_grid_var_names__
+   Grid.__from_ds__
+   Grid.__from_vert__
+   Grid.__init_grid_var_attrs__
 
 Attributes
 ----------
@@ -140,15 +107,34 @@ Operators
 ---------
 .. autosummary::
    :toctree: _autosummary
+
    Grid.__eq__
    Grid.__ne__
 
 Helpers
 ===========
 
-.. currentmodule:: uxarray
+Connectivity
+------------
+.. autosummary::
+   :toctree: _autosummary
+
+   grid.connectivity._replace_fill_values
+   grid.connectivity._build_nNodes_per_face
+   grid.connectivity._build_edge_node_connectivity
+   grid.connectivity._build_face_edges_connectivity
+
+Coordinates
+-----------
+.. autosummary::
+   :toctree: _autosummary
+
+   grid.coordinates._populate_cartesian_xyz_coord
+   grid.coordinates._populate_lonlat_coord
+
+IO
+---
 .. autosummary::
    :toctree: _autosummary
 
-   utils.helpers._is_ugrid
-   utils.helpers._replace_fill_values
+   io.utils._parse_grid_type
diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
index 553a19ef..0c3f4295 100644
--- a/docs/user_api/index.rst
+++ b/docs/user_api/index.rst
@@ -9,8 +9,6 @@ check the draft `UXarray API
 <https://github.com/UXARRAY/uxarray/blob/main/docs/user_api/uxarray_api.md>`_
 documentation to see the tentative whole API and let us know if you have any feedback!
 
-
-
 UxDataset
 =========
 A ``xarray.Dataset``-like, multi-dimensional, in memory, array database.
@@ -20,7 +18,6 @@ dataset operators and attributes through the ``uxgrid`` accessor.
 Below is a list of features explicitly added to `UxDataset` to work on
 Unstructured Grids:
 
-
 Class
 -----
 .. autosummary::
@@ -33,28 +30,24 @@ IO
 .. autosummary::
    :toctree: _autosummary
 
-    open_dataset
-    open_mfdataset
-    enable_jit_cache
-    disable_jit_cache
-    enable_jit
-    disable_jit
+   open_dataset
+   open_mfdataset
 
 Attributes
 ----------
 .. autosummary::
    :toctree: _autosummary
 
-    UxDataset.uxgrid
-    UxDataset.source_datasets
+   UxDataset.uxgrid
+   UxDataset.source_datasets
 
 Methods
 -------
 .. autosummary::
    :toctree: _autosummary
 
-    UxDataset.info
-    UxDataset.integrate
+   UxDataset.info
+   UxDataset.integrate
 
 
 UxDataArray
@@ -139,23 +132,44 @@ Attributes
    Grid.Mesh2_edge_nodes
    Grid.Mesh2_face_edges
 
-
 Helpers
 ===========
 
-.. currentmodule:: uxarray
+Face Area
+----------
+.. autosummary::
+   :toctree: _autosummary
+
+   grid.area.calculate_face_area
+   grid.area.get_all_face_area_from_coords
+   grid.area.calculate_spherical_triangle_jacobian
+   grid.area.calculate_spherical_triangle_jacobian_barycentric
+   grid.area.get_gauss_quadratureDG
+   grid.area.get_tri_quadratureDG
+
+Connectivity
+------------
+.. autosummary::
+   :toctree: _autosummary
+
+   grid.connectivity.close_face_nodes
+
+Coordinates
+-----------
+.. autosummary::
+   :toctree: _autosummary
+
+   grid.coordinates.node_lonlat_rad_to_xyz
+   grid.coordinates.node_xyz_to_lonlat_rad
+   grid.coordinates.normalize_in_place
+   grid.coordinates.grid_center_lat_lon
+
+Numba
+-----
 .. autosummary::
    :toctree: _autosummary
 
-   calculate_face_area
-   calculate_spherical_triangle_jacobian
-   calculate_spherical_triangle_jacobian_barycentric
-   close_face_nodes
-   get_all_face_area_from_coords
-   get_gauss_quadratureDG
-   get_tri_quadratureDG
-   grid_center_lat_lon
-   node_xyz_to_lonlat_rad
-   node_lonlat_rad_to_xyz
-   normalize_in_place
-   parse_grid_type
+   utils.enable_jit_cache
+   utils.disable_jit_cache
+   utils.enable_jit
+   utils.disable_jit
diff --git a/uxarray/__init__.py b/uxarray/__init__.py
index b9aedd24..209ef4df 100644
--- a/uxarray/__init__.py
+++ b/uxarray/__init__.py
@@ -1,11 +1,11 @@
-from uxarray.utils.helpers import *
-from uxarray.utils.constants import (INT_DTYPE, INT_FILL_VALUE)
-from uxarray.core.grid import Grid
+# from uxarray.utils.helpers import *
+from uxarray.constants import (INT_DTYPE, INT_FILL_VALUE)
+from uxarray.grid import Grid
 from uxarray.core.api import (open_grid, open_dataset, open_mfdataset)
 from uxarray.core.dataarray import UxDataArray
 from uxarray.core.dataset import UxDataset
 
-# Sets the version of uxarray currently installed
+# Sets the version of uxarray currently installeds
 # Attempt to import the needed modules
 try:
     from importlib.metadata import version as _version
diff --git a/uxarray/utils/constants.py b/uxarray/constants.py
similarity index 100%
rename from uxarray/utils/constants.py
rename to uxarray/constants.py
diff --git a/uxarray/core/api.py b/uxarray/core/api.py
index e4f09656..1e1d904f 100644
--- a/uxarray/core/api.py
+++ b/uxarray/core/api.py
@@ -1,16 +1,14 @@
 """UXarray dataset module."""
 
 import os
-
-from pathlib import Path
-from typing import Any, Dict, Optional, Union
-
 import numpy as np
 import xarray as xr
 
-import uxarray.utils.constants
-from uxarray.core.grid import Grid
+from pathlib import Path
+from typing import Any, Dict, Optional, Union
 
+import uxarray.constants
+from uxarray.grid import Grid
 from uxarray.core.dataset import UxDataset
 
 
@@ -277,51 +275,3 @@ def open_mfdataset(grid_filename_or_obj: str,
     uxds = UxDataset(ds, uxgrid=uxgrid, source_datasets=str(paths))
 
     return uxds
-
-
-def enable_jit_cache():
-    """Allows Numba's JIT cache to be turned on.
-
-    This cache variable lets @njit cache the machine code generated
-    between runs, allowing for faster run times due to the fact that the
-    code doesn't need to regenerate the machine code every run time. Our
-    use case here was to study performance, in regular usage one might
-    never turn off caching as it will only help if frequently modifying
-    the code or because users have very limited disk space. The default
-    is on (True)
-    """
-    uxarray.utils.constants.ENABLE_JIT_CACHE = True
-
-
-def disable_jit_cache():
-    """Allows Numba's JIT cache to be turned on off.
-
-    This cache variable lets @njit cache the machine code generated
-    between runs, allowing for faster run times due to the fact that the
-    code doesn't need to regenerate the machine code every run time. Our
-    use case here was to study performance, in regular usage one might
-    never turn off caching as it will only help if frequently modifying
-    the code or because users have very limited disk space. The default
-    is on (True)
-    """
-    uxarray.utils.constants.ENABLE_JIT_CACHE = False
-
-
-def enable_jit():
-    """Allows Numba's JIT application to be turned on.
-
-    This lets users choose whether they want machine code to be
-    generated to speed up the performance of the code on large files.
-    The default is on (True)
-    """
-    uxarray.utils.constants.ENABLE_JIT = True
-
-
-def disable_jit():
-    """Allows Numba's JIT application to be turned off.
-
-    This lets users choose whether they want machine code to be
-    generated to speed up the performance of the code on large files.
-    The default is on (True)
-    """
-    uxarray.utils.constants.ENABLE_JIT = False
diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 72f46102..ee98ed8f 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -3,7 +3,7 @@ import xarray as xr
 
 from typing import Optional
 
-from uxarray.core.grid import Grid
+from uxarray.grid.grid import Grid
 
 
 class UxDataArray(xr.DataArray):
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index fdb1097e..4de3369f 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -8,7 +8,7 @@ from collections.abc import Hashable
 from typing import Optional, IO
 
 from uxarray.core.dataarray import UxDataArray
-from uxarray.core.grid import Grid
+from uxarray.grid.grid import Grid
 
 
 class UxDataset(xr.Dataset):
diff --git a/uxarray/grid/__init__.py b/uxarray/grid/__init__.py
new file mode 100644
index 00000000..8e434771
--- /dev/null
+++ b/uxarray/grid/__init__.py
@@ -0,0 +1,1 @@
+from .grid import Grid
diff --git a/uxarray/utils/get_quadratureDG.py b/uxarray/grid/area.py
similarity index 57%
rename from uxarray/utils/get_quadratureDG.py
rename to uxarray/grid/area.py
index 87b4ba7f..30dbb453 100644
--- a/uxarray/utils/get_quadratureDG.py
+++ b/uxarray/grid/area.py
@@ -1,6 +1,318 @@
 import numpy as np
 from numba import njit, config
-from uxarray.utils.constants import ENABLE_JIT_CACHE, ENABLE_JIT
+from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT
+
+from uxarray.grid.coordinates import node_lonlat_rad_to_xyz
+
+config.DISABLE_JIT = not ENABLE_JIT
+
+
+@njit(cache=ENABLE_JIT_CACHE)
+def calculate_face_area(x,
+                        y,
+                        z,
+                        quadrature_rule="gaussian",
+                        order=4,
+                        coords_type="spherical"):
+    """Calculate area of a face on sphere.
+
+    Parameters
+    ----------
+    x : list, required
+        x-coordinate of all the nodes forming the face
+
+    y : list, required
+        y-coordinate of all the nodes forming the face
+
+    z : list, required
+        z-coordinate of all the nodes forming the face
+
+    quadrature_rule : str, optional
+        triangular and Gaussian quadrature supported, expected values: "triangular" or "gaussian"
+
+    order: int, optional
+        Order of the quadrature rule. Default is 4.
+
+        Supported values:
+            - Gaussian Quadrature: 1 to 10
+            - Triangular: 1, 4, 8, 10 and 12
+
+    coords_type : str, optional
+        coordinate type, default is spherical, can be cartesian also.
+    """
+    area = 0.0  # set area to 0
+    order = order
+
+    if quadrature_rule == "gaussian":
+        dG, dW = get_gauss_quadratureDG(order)
+    elif quadrature_rule == "triangular":
+        dG, dW = get_tri_quadratureDG(order)
+    else:
+        raise ValueError(
+            "Invalid quadrature rule, specify gaussian or triangular")
+
+    num_nodes = len(x)
+
+    # num triangles is two less than the total number of nodes
+    num_triangles = num_nodes - 2
+
+    # Using tempestremap GridElements: https://github.com/ClimateGlobalChange/tempestremap/blob/master/src/GridElements.cpp
+    # loop through all sub-triangles of face
+    for j in range(0, num_triangles):
+        node1 = np.array([x[0], y[0], z[0]], dtype=np.float64)
+        node2 = np.array([x[j + 1], y[j + 1], z[j + 1]], dtype=np.float64)
+        node3 = np.array([x[j + 2], y[j + 2], z[j + 2]], dtype=np.float64)
+
+        if (coords_type == "spherical"):
+            node1 = np.array(
+                node_lonlat_rad_to_xyz([np.deg2rad(x[0]),
+                                        np.deg2rad(y[0])]))
+            node2 = np.array(
+                node_lonlat_rad_to_xyz(
+                    [np.deg2rad(x[j + 1]),
+                     np.deg2rad(y[j + 1])]))
+            node3 = np.array(
+                node_lonlat_rad_to_xyz(
+                    [np.deg2rad(x[j + 2]),
+                     np.deg2rad(y[j + 2])]))
+
+        for p in range(len(dW)):
+            if quadrature_rule == "gaussian":
+                for q in range(len(dW)):
+                    dA = dG[0][p]
+                    dB = dG[0][q]
+                    jacobian = calculate_spherical_triangle_jacobian(
+                        node1, node2, node3, dA, dB)
+                    area += dW[p] * dW[q] * jacobian
+            elif quadrature_rule == "triangular":
+                dA = dG[p][0]
+                dB = dG[p][1]
+                jacobian = calculate_spherical_triangle_jacobian_barycentric(
+                    node1, node2, node3, dA, dB)
+                area += dW[p] * jacobian
+
+    return area
+
+
+@njit(cache=ENABLE_JIT_CACHE)
+def get_all_face_area_from_coords(x,
+                                  y,
+                                  z,
+                                  face_nodes,
+                                  face_geometry,
+                                  dim,
+                                  quadrature_rule="triangular",
+                                  order=4,
+                                  coords_type="spherical"):
+    """Given coords, connectivity and other area calculation params, this
+    routine loop over all faces and return an numpy array with areas of each
+    face.
+
+    Parameters
+    ----------
+    x : ndarray, required
+        x-coordinate of all the nodes
+
+    y : ndarray, required
+        y-coordinate of all the nodes
+
+    z : ndarray, required
+        z-coordinate of all the nodes
+
+    face_nodes : 2D ndarray, required
+         node ids of each face
+
+    dim : int, required
+         dimension
+
+    quadrature_rule : str, optional
+        "triangular" or "gaussian". Defaults to triangular
+
+    order : int, optional
+        count or order for Gaussian or spherical resp. Defaults to 4 for spherical.
+
+    coords_type : str, optional
+        coordinate type, default is spherical, can be cartesian also.
+
+    Returns
+    -------
+    area of all faces : ndarray
+    """
+
+    n_face, n_max_face_nodes = face_nodes.shape
+
+    # set initial area of each face to 0
+    area = np.zeros(n_face)
+
+    for face_idx, max_nodes in enumerate(face_geometry):
+        face_x = x[face_nodes[face_idx, 0:max_nodes]]
+        face_y = y[face_nodes[face_idx, 0:max_nodes]]
+
+        # check if z dimension
+
+        if dim > 2:
+            face_z = z[face_nodes[face_idx, 0:max_nodes]]
+        else:
+            face_z = face_x * 0.0
+
+        # After getting all the nodes of a face assembled call the  cal. face area routine
+        face_area = calculate_face_area(face_x, face_y, face_z, quadrature_rule,
+                                        order, coords_type)
+        # store current face area
+        area[face_idx] = face_area
+
+    return area
+
+
+@njit(cache=ENABLE_JIT_CACHE)
+def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB):
+    """Calculate Jacobian of a spherical triangle. This is a helper function
+    for calculating face area.
+
+    Parameters
+    ----------
+    node1 : list, required
+        First node of the triangle
+
+    node2 : list, required
+        Second node of the triangle
+
+    node3 : list, required
+        Third node of the triangle
+
+    dA : float, required
+        quadrature point
+
+    dB : float, required
+        quadrature point
+
+    Returns
+    -------
+    jacobian : float
+    """
+    dF = np.array([
+        (1.0 - dB) * ((1.0 - dA) * node1[0] + dA * node2[0]) + dB * node3[0],
+        (1.0 - dB) * ((1.0 - dA) * node1[1] + dA * node2[1]) + dB * node3[1],
+        (1.0 - dB) * ((1.0 - dA) * node1[2] + dA * node2[2]) + dB * node3[2]
+    ])
+
+    dDaF = np.array([(1.0 - dB) * (node2[0] - node1[0]),
+                     (1.0 - dB) * (node2[1] - node1[1]),
+                     (1.0 - dB) * (node2[2] - node1[2])])
+
+    dDbF = np.array([
+        -(1.0 - dA) * node1[0] - dA * node2[0] + node3[0],
+        -(1.0 - dA) * node1[1] - dA * node2[1] + node3[1],
+        -(1.0 - dA) * node1[2] - dA * node2[2] + node3[2]
+    ])
+
+    dInvR = 1.0 / np.sqrt(dF[0] * dF[0] + dF[1] * dF[1] + dF[2] * dF[2])
+
+    dDaG = np.array([
+        dDaF[0] * (dF[1] * dF[1] + dF[2] * dF[2]) - dF[0] *
+        (dDaF[1] * dF[1] + dDaF[2] * dF[2]),
+        dDaF[1] * (dF[0] * dF[0] + dF[2] * dF[2]) - dF[1] *
+        (dDaF[0] * dF[0] + dDaF[2] * dF[2]),
+        dDaF[2] * (dF[0] * dF[0] + dF[1] * dF[1]) - dF[2] *
+        (dDaF[0] * dF[0] + dDaF[1] * dF[1])
+    ])
+
+    dDbG = np.array([
+        dDbF[0] * (dF[1] * dF[1] + dF[2] * dF[2]) - dF[0] *
+        (dDbF[1] * dF[1] + dDbF[2] * dF[2]),
+        dDbF[1] * (dF[0] * dF[0] + dF[2] * dF[2]) - dF[1] *
+        (dDbF[0] * dF[0] + dDbF[2] * dF[2]),
+        dDbF[2] * (dF[0] * dF[0] + dF[1] * dF[1]) - dF[2] *
+        (dDbF[0] * dF[0] + dDbF[1] * dF[1])
+    ])
+
+    dDenomTerm = dInvR * dInvR * dInvR
+
+    dDaG *= dDenomTerm
+    dDbG *= dDenomTerm
+
+    #  Cross product gives local Jacobian
+    nodeCross = np.cross(dDaG, dDbG)
+    dJacobian = np.sqrt(nodeCross[0] * nodeCross[0] +
+                        nodeCross[1] * nodeCross[1] +
+                        nodeCross[2] * nodeCross[2])
+
+    return dJacobian
+
+
+@njit(cache=ENABLE_JIT_CACHE)
+def calculate_spherical_triangle_jacobian_barycentric(node1, node2, node3, dA,
+                                                      dB):
+    """Calculate Jacobian of a spherical triangle. This is a helper function
+    for calculating face area.
+
+    Parameters
+    ----------
+    node1 : list, required
+        First node of the triangle
+
+    node2 : list, required
+        Second node of the triangle
+
+    node3 : list, required
+        Third node of the triangle
+
+    dA : float, required
+        first component of barycentric coordinates of quadrature point
+
+    dB : float, required
+        second component of barycentric coordinates of quadrature point
+
+    Returns
+    -------
+    jacobian : float
+    """
+
+    dF = np.array([
+        dA * node1[0] + dB * node2[0] + (1.0 - dA - dB) * node3[0],
+        dA * node1[1] + dB * node2[1] + (1.0 - dA - dB) * node3[1],
+        dA * node1[2] + dB * node2[2] + (1.0 - dA - dB) * node3[2]
+    ])
+
+    dDaF = np.array(
+        [node1[0] - node3[0], node1[1] - node3[1], node1[2] - node3[2]])
+
+    dDbF = np.array(
+        [node2[0] - node3[0], node2[1] - node3[1], node2[2] - node3[2]])
+
+    dInvR = 1.0 / np.sqrt(dF[0] * dF[0] + dF[1] * dF[1] + dF[2] * dF[2])
+
+    dDaG = np.array([
+        dDaF[0] * (dF[1] * dF[1] + dF[2] * dF[2]) - dF[0] *
+        (dDaF[1] * dF[1] + dDaF[2] * dF[2]),
+        dDaF[1] * (dF[0] * dF[0] + dF[2] * dF[2]) - dF[1] *
+        (dDaF[0] * dF[0] + dDaF[2] * dF[2]),
+        dDaF[2] * (dF[0] * dF[0] + dF[1] * dF[1]) - dF[2] *
+        (dDaF[0] * dF[0] + dDaF[1] * dF[1])
+    ])
+
+    dDbG = np.array([
+        dDbF[0] * (dF[1] * dF[1] + dF[2] * dF[2]) - dF[0] *
+        (dDbF[1] * dF[1] + dDbF[2] * dF[2]),
+        dDbF[1] * (dF[0] * dF[0] + dF[2] * dF[2]) - dF[1] *
+        (dDbF[0] * dF[0] + dDbF[2] * dF[2]),
+        dDbF[2] * (dF[0] * dF[0] + dF[1] * dF[1]) - dF[2] *
+        (dDbF[0] * dF[0] + dDbF[1] * dF[1])
+    ])
+
+    dDenomTerm = dInvR * dInvR * dInvR
+
+    dDaG *= dDenomTerm
+    dDbG *= dDenomTerm
+
+    #  Cross product gives local Jacobian
+    nodeCross = np.cross(dDaG, dDbG)
+    dJacobian = np.sqrt(nodeCross[0] * nodeCross[0] +
+                        nodeCross[1] * nodeCross[1] +
+                        nodeCross[2] * nodeCross[2])
+
+    return 0.5 * dJacobian
+
 
 config.DISABLE_JIT = not ENABLE_JIT
 
diff --git a/uxarray/grid/connectivity.py b/uxarray/grid/connectivity.py
new file mode 100644
index 00000000..0cea8ada
--- /dev/null
+++ b/uxarray/grid/connectivity.py
@@ -0,0 +1,232 @@
+import numpy as np
+import xarray as xr
+
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
+
+
+def close_face_nodes(Mesh2_face_nodes, nMesh2_face, nMaxMesh2_face_nodes):
+    """Closes (``Mesh2_face_nodes``) by inserting the first node index after
+    the last non-fill-value node.
+
+    Parameters
+    ----------
+    Mesh2_face_nodes : np.ndarray
+        Connectivity array for constructing a face from its nodes
+    nMesh2_face : constant
+        Number of faces
+    nMaxMesh2_face_nodes : constant
+        Max number of nodes that compose a face
+
+    Returns
+    ----------
+    closed : ndarray
+        Closed (padded) Mesh2_face_nodes
+
+    Example
+    ----------
+    Given face nodes with shape [2 x 5]
+        [0, 1, 2, 3, FILL_VALUE]
+        [4, 5, 6, 7, 8]
+    Pads them to the following with shape [2 x 6]
+        [0, 1, 2, 3, 0, FILL_VALUE]
+        [4, 5, 6, 7, 8, 4]
+    """
+
+    # padding to shape [nMesh2_face, nMaxMesh2_face_nodes + 1]
+    closed = np.ones((nMesh2_face, nMaxMesh2_face_nodes + 1),
+                     dtype=INT_DTYPE) * INT_FILL_VALUE
+
+    # set all non-paded values to original face nodee values
+    closed[:, :-1] = Mesh2_face_nodes.copy()
+
+    # instance of first fill value
+    first_fv_idx_2d = np.argmax(closed == INT_FILL_VALUE, axis=1)
+
+    # 2d to 1d index for np.put()
+    first_fv_idx_1d = first_fv_idx_2d + (
+        (nMaxMesh2_face_nodes + 1) * np.arange(0, nMesh2_face))
+
+    # column of first node values
+    first_node_value = Mesh2_face_nodes[:, 0].copy()
+
+    # insert first node column at occurrence of first fill value
+    np.put(closed.ravel(), first_fv_idx_1d, first_node_value)
+
+    return closed
+
+
+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, optional
+        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
+    if original_fill is not None and np.isnan(original_fill):
+        fill_val_idx = np.isnan(grid_var)
+    else:
+        fill_val_idx = grid_var == 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 fill value can be represented with current integer data type
+    if np.issubdtype(new_dtype, np.integer):
+        int_min = np.iinfo(grid_var.dtype).min
+        int_max = np.iinfo(grid_var.dtype).max
+        # ensure new_fill is in range [int_min, int_max]
+        if new_fill < int_min or new_fill > int_max:
+            raise ValueError(f'New fill value: {new_fill} not representable by'
+                             f' integer dtype: {grid_var.dtype}')
+
+    # ensure non-nan fill value can be represented with current float data type
+    elif np.issubdtype(new_dtype, np.floating) and not np.isnan(new_fill):
+        float_min = np.finfo(grid_var.dtype).min
+        float_max = np.finfo(grid_var.dtype).max
+        # ensure new_fill is in range [float_min, float_max]
+        if new_fill < float_min or new_fill > float_max:
+            raise ValueError(f'New fill value: {new_fill} not representable by'
+                             f' float dtype: {grid_var.dtype}')
+    else:
+        raise ValueError(f'Data type {grid_var.dtype} not supported'
+                         f'for grid variables')
+
+    # replace all zeros with a fill value
+    grid_var[fill_val_idx] = new_fill
+
+    return grid_var
+
+
+def _build_nNodes_per_face(grid):
+    """Constructs ``nNodes_per_face``, which contains the number of non- fill-
+    value nodes for each face in ``Mesh2_face_nodes``"""
+
+    # padding to shape [nMesh2_face, nMaxMesh2_face_nodes + 1]
+    closed = np.ones((grid.nMesh2_face, grid.nMaxMesh2_face_nodes + 1),
+                     dtype=INT_DTYPE) * INT_FILL_VALUE
+
+    closed[:, :-1] = grid.Mesh2_face_nodes.copy()
+
+    nNodes_per_face = np.argmax(closed == INT_FILL_VALUE, axis=1)
+
+    # add to internal dataset
+    grid._ds["nNodes_per_face"] = xr.DataArray(
+        data=nNodes_per_face,
+        dims=["nMesh2_face"],
+        attrs={"long_name": "number of non-fill value nodes for each face"})
+
+
+def _build_edge_node_connectivity(grid, repopulate=False):
+    """Constructs the UGRID connectivity variable (``Mesh2_edge_nodes``) and
+    stores it within the internal (``Grid._ds``) and through the attribute
+    (``Grid.Mesh2_edge_nodes``).
+
+    Additionally, the attributes (``inverse_indices``) and
+    (``fill_value_mask``) are stored for constructing other
+    connectivity variables.
+
+    Parameters
+    ----------
+    repopulate : bool, optional
+        Flag used to indicate if we want to overwrite the existed `Mesh2_edge_nodes` and generate a new
+        inverse_indices, default is False
+    """
+
+    # need to derive edge nodes
+    if "Mesh2_edge_nodes" not in grid._ds or repopulate:
+        padded_face_nodes = close_face_nodes(grid.Mesh2_face_nodes.values,
+                                             grid.nMesh2_face,
+                                             grid.nMaxMesh2_face_nodes)
+
+        # array of empty edge nodes where each entry is a pair of indices
+        edge_nodes = np.empty((grid.nMesh2_face * grid.nMaxMesh2_face_nodes, 2),
+                              dtype=INT_DTYPE)
+
+        # first index includes starting node up to non-padded value
+        edge_nodes[:, 0] = padded_face_nodes[:, :-1].ravel()
+
+        # second index includes second node up to padded value
+        edge_nodes[:, 1] = padded_face_nodes[:, 1:].ravel()
+    else:
+        # If "Mesh2_edge_nodes" already exists, directly return the function call
+        return
+
+    # sorted edge nodes
+    edge_nodes.sort(axis=1)
+
+    # unique edge nodes
+    edge_nodes_unique, inverse_indices = np.unique(edge_nodes,
+                                                   return_inverse=True,
+                                                   axis=0)
+    # find all edge nodes that contain a fill value
+    fill_value_mask = np.logical_or(edge_nodes_unique[:, 0] == INT_FILL_VALUE,
+                                    edge_nodes_unique[:, 1] == INT_FILL_VALUE)
+
+    # all edge nodes that do not contain a fill value
+    non_fill_value_mask = np.logical_not(fill_value_mask)
+    edge_nodes_unique = edge_nodes_unique[non_fill_value_mask]
+
+    # Update inverse_indices accordingly
+    indices_to_update = np.where(fill_value_mask)[0]
+
+    remove_mask = np.isin(inverse_indices, indices_to_update)
+    inverse_indices[remove_mask] = INT_FILL_VALUE
+
+    # Compute the indices where inverse_indices exceeds the values in indices_to_update
+    indexes = np.searchsorted(indices_to_update, inverse_indices, side='right')
+    # subtract the corresponding indexes from `inverse_indices`
+    for i in range(len(inverse_indices)):
+        if inverse_indices[i] != INT_FILL_VALUE:
+            inverse_indices[i] -= indexes[i]
+
+    # add Mesh2_edge_nodes to internal dataset
+    grid._ds['Mesh2_edge_nodes'] = xr.DataArray(
+        edge_nodes_unique,
+        dims=["nMesh2_edge", "Two"],
+        attrs={
+            "cf_role": "edge_node_connectivity",
+            "_FillValue": INT_FILL_VALUE,
+            "long_name": "Maps every edge to the two nodes that it connects",
+            "start_index": INT_DTYPE(0),
+            "inverse_indices": inverse_indices,
+            "fill_value_mask": fill_value_mask
+        })
+
+
+def _build_face_edges_connectivity(grid):
+    """Constructs the UGRID connectivity variable (``Mesh2_face_edges``) and
+    stores it within the internal (``Grid._ds``) and through the attribute
+    (``Grid.Mesh2_face_edges``)."""
+    if ("Mesh2_edge_nodes" not in grid._ds or
+            "inverse_indices" not in grid._ds['Mesh2_edge_nodes'].attrs):
+        _build_edge_node_connectivity(grid, repopulate=True)
+
+    inverse_indices = grid._ds['Mesh2_edge_nodes'].inverse_indices
+    inverse_indices = inverse_indices.reshape(grid.nMesh2_face,
+                                              grid.nMaxMesh2_face_nodes)
+
+    grid._ds["Mesh2_face_edges"] = xr.DataArray(
+        data=inverse_indices,
+        dims=["nMesh2_face", "nMaxMesh2_face_edges"],
+        attrs={
+            "cf_role": "face_edges_connectivity",
+            "start_index": INT_DTYPE(0),
+            "long_name": "Maps every edge to the two nodes that it connects",
+        })
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
new file mode 100644
index 00000000..545105f0
--- /dev/null
+++ b/uxarray/grid/coordinates.py
@@ -0,0 +1,288 @@
+import xarray as xr
+import numpy as np
+import math
+
+from numba import njit, config
+
+from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT
+
+config.DISABLE_JIT = not ENABLE_JIT
+
+
+@njit(cache=ENABLE_JIT_CACHE)
+def node_lonlat_rad_to_xyz(node_coord):
+    """Helper function to Convert the node coordinate from 2D
+    longitude/latitude to normalized 3D xyz.
+
+    Parameters
+    ----------
+    node: float list
+        2D coordinates[longitude, latitude] in radiance
+
+    Returns
+    ----------
+    float list
+        the result array of the unit 3D coordinates [x, y, z] vector where :math:`x^2 + y^2 + z^2 = 1`
+
+    Raises
+    ----------
+    RuntimeError
+        The input array doesn't have the size of 3.
+    """
+    if len(node_coord) != 2:
+        raise RuntimeError(
+            "Input array should have a length of 2: [longitude, latitude]")
+    lon = node_coord[0]
+    lat = node_coord[1]
+    return [np.cos(lon) * np.cos(lat), np.sin(lon) * np.cos(lat), np.sin(lat)]
+
+
+@njit(cache=ENABLE_JIT_CACHE)
+def node_xyz_to_lonlat_rad(node_coord):
+    """Calculate the latitude and longitude in radiance for a node represented
+    in the [x, y, z] 3D Cartesian coordinates.
+
+    Parameters
+    ----------
+    node_coord: float list
+        3D Cartesian Coordinates [x, y, z] of the node
+
+    Returns
+    ----------
+    float list
+        the result array of longitude and latitude in radian [longitude_rad, latitude_rad]
+
+    Raises
+    ----------
+    RuntimeError
+        The input array doesn't have the size of 3.
+    """
+    if len(node_coord) != 3:
+        raise RuntimeError("Input array should have a length of 3: [x, y, z]")
+    reference_tolerance = 1.0e-12
+    [dx, dy, dz] = normalize_in_place(node_coord)
+    dx /= np.absolute(dx * dx + dy * dy + dz * dz)
+    dy /= np.absolute(dx * dx + dy * dy + dz * dz)
+    dz /= np.absolute(dx * dx + dy * dy + dz * dz)
+
+    if np.absolute(dz) < (1.0 - reference_tolerance):
+        d_lon_rad = math.atan2(dy, dx)
+        d_lat_rad = np.arcsin(dz)
+
+        if d_lon_rad < 0.0:
+            d_lon_rad += 2.0 * np.pi
+    elif dz > 0.0:
+        d_lon_rad = 0.0
+        d_lat_rad = 0.5 * np.pi
+    else:
+        d_lon_rad = 0.0
+        d_lat_rad = -0.5 * np.pi
+
+    return [d_lon_rad, d_lat_rad]
+
+
+@njit(cache=ENABLE_JIT_CACHE)
+def normalize_in_place(node):
+    """Helper function to project an arbitrary node in 3D coordinates [x, y, z]
+    on the unit sphere. It uses the `np.linalg.norm` internally to calculate
+    the magnitude.
+
+    Parameters
+    ----------
+    node: float list
+        3D Cartesian Coordinates [x, y, z]
+
+    Returns
+    ----------
+    float list
+        the result unit vector [x, y, z] where :math:`x^2 + y^2 + z^2 = 1`
+
+    Raises
+    ----------
+    RuntimeError
+        The input array doesn't have the size of 3.
+    """
+    if len(node) != 3:
+        raise RuntimeError("Input array should have a length of 3: [x, y, z]")
+
+    return np.array(node) / np.linalg.norm(np.array(node), ord=2)
+
+
+def grid_center_lat_lon(ds):
+    """Using scrip file variables ``grid_corner_lat`` and ``grid_corner_lon``,
+    calculates the ``grid_center_lat`` and ``grid_center_lon``.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        Dataset that contains ``grid_corner_lat`` and ``grid_corner_lon``
+        data variables
+
+    Returns
+    -------
+    center_lon : :class:`numpy.ndarray`
+        The calculated center longitudes of the grid box based on the corner
+        points
+    center_lat : :class:`numpy.ndarray`
+        The calculated center latitudes of the grid box based on the corner
+        points
+    """
+
+    # Calculate and create grid center lat/lon
+    scrip_corner_lon = ds['grid_corner_lon']
+    scrip_corner_lat = ds['grid_corner_lat']
+
+    # convert to radians
+    rad_corner_lon = np.deg2rad(scrip_corner_lon)
+    rad_corner_lat = np.deg2rad(scrip_corner_lat)
+
+    # get nodes per face
+    nodes_per_face = rad_corner_lat.shape[1]
+
+    # geographic center of each cell
+    x = np.sum(np.cos(rad_corner_lat) * np.cos(rad_corner_lon),
+               axis=1) / nodes_per_face
+    y = np.sum(np.cos(rad_corner_lat) * np.sin(rad_corner_lon),
+               axis=1) / nodes_per_face
+    z = np.sum(np.sin(rad_corner_lat), axis=1) / nodes_per_face
+
+    center_lon = np.rad2deg(np.arctan2(y, x))
+    center_lat = np.rad2deg(np.arctan2(z, np.sqrt(x**2 + y**2)))
+
+    # Make negative lons positive
+    center_lon[center_lon < 0] += 360
+
+    return center_lat, center_lon
+
+
+def _populate_cartesian_xyz_coord(grid):
+    """A helper function that populates the xyz attribute in UXarray.Grid._ds.
+    This function is called when we need to use the cartesian coordinates for
+    each node to do the calculation but the input data only has the
+    "Mesh2_node_x" and "Mesh2_node_y" in degree.
+
+    Note
+    ----
+    In the UXarray, we abide the UGRID convention and make sure the following attributes will always have its
+    corresponding units as stated below:
+
+    Mesh2_node_x
+     unit:  "degree_east" for longitude
+    Mesh2_node_y
+     unit:  "degrees_north" for latitude
+    Mesh2_node_z
+     unit:  "m"
+    Mesh2_node_cart_x
+     unit:  "m"
+    Mesh2_node_cart_y
+     unit:  "m"
+    Mesh2_node_cart_z
+     unit:  "m"
+    """
+
+    # Check if the cartesian coordinates are already populated
+    if "Mesh2_node_cart_x" in grid._ds.keys():
+        return
+
+    # check for units and create Mesh2_node_cart_x/y/z set to grid._ds
+    nodes_lon_rad = np.deg2rad(grid.Mesh2_node_x.values)
+    nodes_lat_rad = np.deg2rad(grid.Mesh2_node_y.values)
+    nodes_rad = np.stack((nodes_lon_rad, nodes_lat_rad), axis=1)
+    nodes_cart = np.asarray(list(map(node_lonlat_rad_to_xyz, list(nodes_rad))))
+
+    grid._ds["Mesh2_node_cart_x"] = xr.DataArray(
+        data=nodes_cart[:, 0],
+        dims=["nMesh2_node"],
+        attrs={
+            "standard_name": "cartesian x",
+            "units": "m",
+        })
+    grid._ds["Mesh2_node_cart_y"] = xr.DataArray(
+        data=nodes_cart[:, 1],
+        dims=["nMesh2_node"],
+        attrs={
+            "standard_name": "cartesian y",
+            "units": "m",
+        })
+    grid._ds["Mesh2_node_cart_z"] = xr.DataArray(
+        data=nodes_cart[:, 2],
+        dims=["nMesh2_node"],
+        attrs={
+            "standard_name": "cartesian z",
+            "units": "m",
+        })
+
+
+def _populate_lonlat_coord(grid):
+    """Helper function that populates the longitude and latitude and store it
+    into the Mesh2_node_x and Mesh2_node_y. This is called when the input data
+    has "Mesh2_node_x", "Mesh2_node_y", "Mesh2_node_z" in meters. Since we want
+    "Mesh2_node_x" and "Mesh2_node_y" always have the "degree" units. For more
+    details, please read the following.
+
+    Raises
+    ------
+        RuntimeError
+            Mesh2_node_x/y/z are not represented in the cartesian format with the unit 'm'/'meters' when calling this function"
+
+    Note
+    ----
+    In the UXarray, we abide the UGRID convention and make sure the following attributes will always have its
+    corresponding units as stated below:
+
+    Mesh2_node_x
+     unit:  "degree_east" for longitude
+    Mesh2_node_y
+     unit:  "degrees_north" for latitude
+    Mesh2_node_z
+     unit:  "m"
+    Mesh2_node_cart_x
+     unit:  "m"
+    Mesh2_node_cart_y
+     unit:  "m"
+    Mesh2_node_cart_z
+     unit:  "m"
+    """
+
+    # Check if the "Mesh2_node_x" is already in longitude
+    if "degree" in grid._ds.Mesh2_node_x.units:
+        return
+
+    # Check if the input Mesh2_node_xyz" are represented in the cartesian format with the unit "m"
+    if ("m" not in grid._ds.Mesh2_node_x.units) or ("m" not in grid._ds.Mesh2_node_y.units) \
+            or ("m" not in grid._ds.Mesh2_node_z.units):
+        raise RuntimeError(
+            "Expected: Mesh2_node_x/y/z should be represented in the cartesian format with the "
+            "unit 'm' when calling this function")
+
+    # Put the cartesian coordinates inside the proper data structure
+    grid._ds["Mesh2_node_cart_x"] = xr.DataArray(
+        data=grid._ds["Mesh2_node_x"].values)
+    grid._ds["Mesh2_node_cart_y"] = xr.DataArray(
+        data=grid._ds["Mesh2_node_y"].values)
+    grid._ds["Mesh2_node_cart_z"] = xr.DataArray(
+        data=grid._ds["Mesh2_node_z"].values)
+
+    # convert the input cartesian values into the longitude latitude degree
+    nodes_cart = np.stack(
+        (grid._ds["Mesh2_node_x"].values, grid._ds["Mesh2_node_y"].values,
+         grid._ds["Mesh2_node_z"].values),
+        axis=1).tolist()
+    nodes_rad = list(map(node_xyz_to_lonlat_rad, nodes_cart))
+    nodes_degree = np.rad2deg(nodes_rad)
+    grid._ds["Mesh2_node_x"] = xr.DataArray(
+        data=nodes_degree[:, 0],
+        dims=["nMesh2_node"],
+        attrs={
+            "standard_name": "longitude",
+            "long_name": "longitude of mesh nodes",
+            "units": "degrees_east",
+        })
+    grid._ds["Mesh2_node_y"] = xr.DataArray(
+        data=nodes_degree[:, 1],
+        dims=["nMesh2_node"],
+        attrs={
+            "standard_name": "latitude",
+            "long_name": "latitude of mesh nodes",
+            "units": "degrees_north",
+        })
diff --git a/uxarray/core/grid.py b/uxarray/grid/grid.py
similarity index 69%
rename from uxarray/core/grid.py
rename to uxarray/grid/grid.py
index 76f76d0d..a09522cc 100644
--- a/uxarray/core/grid.py
+++ b/uxarray/grid/grid.py
@@ -8,10 +8,17 @@ from uxarray.io._mpas import _read_mpas
 from uxarray.io._ugrid import _read_ugrid, _encode_ugrid
 from uxarray.io._shapefile import _read_shpfile
 from uxarray.io._scrip import _read_scrip, _encode_scrip
-from uxarray.utils.helpers import (get_all_face_area_from_coords,
-                                   parse_grid_type, node_xyz_to_lonlat_rad,
-                                   node_lonlat_rad_to_xyz, close_face_nodes)
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.io.utils import _parse_grid_type
+
+from uxarray.grid.area import get_all_face_area_from_coords
+from uxarray.grid.connectivity import (_build_edge_node_connectivity,
+                                       _build_face_edges_connectivity,
+                                       _build_nNodes_per_face)
+
+from uxarray.grid.coordinates import (_populate_lonlat_coord,
+                                      _populate_cartesian_xyz_coord)
+
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 
 class Grid:
@@ -97,7 +104,7 @@ class Grid:
         # check if initializing from string
         # TODO: re-add gridspec initialization when implemented
         elif isinstance(input_obj, xr.Dataset):
-            self.mesh_type = parse_grid_type(input_obj)
+            self.mesh_type = _parse_grid_type(input_obj)
             self.__from_ds__(dataset=input_obj)
         else:
             raise RuntimeError("Dataset is not a valid input type.")
@@ -297,7 +304,47 @@ class Grid:
         return prefix + original_grid_str + dims_heading + dims_str + coord_heading + coords_str + \
             connectivity_heading + connectivity_str
 
-    # attribute properties
+    def __eq__(self, other):
+        """Two grids are equal if they have matching grid topology variables,
+        coordinates, and dims all of which are equal.
+
+        Parameters
+        ----------
+        other : uxarray.Grid
+            The second grid object to be compared with `self`
+
+        Returns
+        -------
+        If two grids are equal : bool
+        """
+        if other is not None:
+            # Iterate over dict to set access attributes
+            for key, value in self.grid_var_names.items():
+                # Check if all grid variables are equal
+                if self._ds.data_vars is not None:
+                    if value in self._ds.data_vars:
+                        if not self._ds[value].equals(
+                                other._ds[other.grid_var_names[key]]):
+                            return False
+        else:
+            return False
+
+        return True
+
+    def __ne__(self, other):
+        """Two grids are not equal if they have differing grid topology
+        variables, coordinates, or dims.
+
+        Parameters
+        ----------
+        other : uxarray.Grid
+            The second grid object to be compared with `self`
+
+        Returns
+        -------
+        If two grids are not equal : bool
+        """
+        return not self.__eq__(other)
 
     @property
     def parsed_attrs(self):
@@ -310,8 +357,6 @@ class Grid:
         unstructured mesh."""
         return self._ds[self.grid_var_names["Mesh2"]]
 
-    # dimension properties
-
     @property
     def nMesh2_node(self):
         """UGRID Dimension ``nMesh2_node``, which represents the total number
@@ -330,7 +375,7 @@ class Grid:
         of edges."""
 
         if "Mesh2_edge_nodes" not in self._ds:
-            self._build_edge_node_connectivity(repopulate=True)
+            _build_edge_node_connectivity(self, repopulate=True)
 
         return self._ds['Mesh2_edge_nodes'].shape[0]
 
@@ -349,7 +394,7 @@ class Grid:
         """
 
         if "Mesh2_face_edges" not in self._ds:
-            self._build_face_edges_connectivity()
+            _build_face_edges_connectivity(self)
 
         return self._ds["Mesh2_face_edges"].shape[1]
 
@@ -361,11 +406,9 @@ class Grid:
         Dimensions (``nMesh2_nodes``) and DataType ``INT_DTYPE``.
         """
         if "nNodes_per_face" not in self._ds:
-            self._build_nNodes_per_face()
+            _build_nNodes_per_face(self)
         return self._ds["nNodes_per_face"]
 
-    # coordinate properties
-
     @property
     def Mesh2_node_x(self):
         """UGRID Coordinate Variable ``Mesh2_node_x``, which contains the
@@ -383,7 +426,7 @@ class Grid:
         Dimensions (``nMesh2_node``)
         """
         if "Mesh2_node_cart_x" not in self._ds:
-            self._populate_cartesian_xyz_coord()
+            _populate_cartesian_xyz_coord(self)
         return self._ds['Mesh2_node_cart_x']
 
     @property
@@ -415,7 +458,7 @@ class Grid:
         Dimensions (``nMesh2_node``)
         """
         if "Mesh2_node_cart_y" not in self._ds:
-            self._populate_cartesian_xyz_coord()
+            _populate_cartesian_xyz_coord(self)
         return self._ds['Mesh2_node_cart_y']
 
     @property
@@ -459,8 +502,6 @@ class Grid:
             self._populate_cartesian_xyz_coord()
         return self._ds['Mesh2_node_cart_z']
 
-    # connectivity properties
-
     @property
     def Mesh2_face_nodes(self):
         """UGRID Connectivity Variable ``Mesh2_face_nodes``, which maps each
@@ -488,7 +529,7 @@ class Grid:
         Nodes are in arbitrary order.
         """
         if "Mesh2_edge_nodes" not in self._ds:
-            self._build_edge_node_connectivity()
+            _build_edge_node_connectivity(self)
 
         return self._ds['Mesh2_edge_nodes']
 
@@ -501,10 +542,18 @@ class Grid:
         DataType ``INT_DTYPE``.
         """
         if "Mesh2_face_edges" not in self._ds:
-            self._build_face_edges_connectivity()
+            _build_face_edges_connectivity(self)
 
         return self._ds["Mesh2_face_edges"]
 
+    @property
+    def face_areas(self):
+        """Declare face_areas as a property."""
+        # if self._face_areas is not None: it allows for using the cached result
+        if self._face_areas is None:
+            self.compute_face_areas()
+        return self._face_areas
+
     def copy(self):
         """Returns a deep copy of this grid."""
         return Grid(xr.Dataset(self._ds),
@@ -634,57 +683,6 @@ class Grid:
 
         return self._face_areas
 
-    def __eq__(self, other):
-        """Two grids are equal if they have matching grid topology variables,
-        coordinates, and dims all of which are equal.
-
-        Parameters
-        ----------
-        other : uxarray.Grid
-            The second grid object to be compared with `self`
-
-        Returns
-        -------
-        If two grids are equal : bool
-        """
-        if other is not None:
-            # Iterate over dict to set access attributes
-            for key, value in self.grid_var_names.items():
-                # Check if all grid variables are equal
-                if self._ds.data_vars is not None:
-                    if value in self._ds.data_vars:
-                        if not self._ds[value].equals(
-                                other._ds[other.grid_var_names[key]]):
-                            return False
-        else:
-            return False
-
-        return True
-
-    def __ne__(self, other):
-        """Two grids are not equal if they have differing grid topology
-        variables, coordinates, or dims.
-
-        Parameters
-        ----------
-        other : uxarray.Grid
-            The second grid object to be compared with `self`
-
-        Returns
-        -------
-        If two grids are not equal : bool
-        """
-        return not self.__eq__(other)
-
-    # use the property keyword for declaration on face_areas property
-    @property
-    def face_areas(self):
-        """Declare face_areas as a property."""
-        # if self._face_areas is not None: it allows for using the cached result
-        if self._face_areas is None:
-            self.compute_face_areas()
-        return self._face_areas
-
     # TODO: Make a decision on whether to provide Dataset- or DataArray-specific
     # functions from within Grid
     # def integrate(self, var_ds, quadrature_rule="triangular", order=4):
@@ -730,263 +728,3 @@ class Grid:
     #     integral = np.dot(face_areas, face_vals)
     #
     #     return integral
-
-    def _build_edge_node_connectivity(self, repopulate=False):
-        """Constructs the UGRID connectivity variable (``Mesh2_edge_nodes``)
-        and stores it within the internal (``Grid._ds``) and through the
-        attribute (``Grid.Mesh2_edge_nodes``).
-
-        Additionally, the attributes (``inverse_indices``) and
-        (``fill_value_mask``) are stored for constructing other
-        connectivity variables.
-
-        Parameters
-        ----------
-        repopulate : bool, optional
-            Flag used to indicate if we want to overwrite the existed `Mesh2_edge_nodes` and generate a new
-            inverse_indices, default is False
-        """
-
-        # need to derive edge nodes
-        if "Mesh2_edge_nodes" not in self._ds or repopulate:
-            padded_face_nodes = close_face_nodes(self.Mesh2_face_nodes.values,
-                                                 self.nMesh2_face,
-                                                 self.nMaxMesh2_face_nodes)
-
-            # array of empty edge nodes where each entry is a pair of indices
-            edge_nodes = np.empty(
-                (self.nMesh2_face * self.nMaxMesh2_face_nodes, 2),
-                dtype=INT_DTYPE)
-
-            # first index includes starting node up to non-padded value
-            edge_nodes[:, 0] = padded_face_nodes[:, :-1].ravel()
-
-            # second index includes second node up to padded value
-            edge_nodes[:, 1] = padded_face_nodes[:, 1:].ravel()
-        else:
-            # If "Mesh2_edge_nodes" already exists, directly return the function call
-            return
-
-        # sorted edge nodes
-        edge_nodes.sort(axis=1)
-
-        # unique edge nodes
-        edge_nodes_unique, inverse_indices = np.unique(edge_nodes,
-                                                       return_inverse=True,
-                                                       axis=0)
-        # find all edge nodes that contain a fill value
-        fill_value_mask = np.logical_or(
-            edge_nodes_unique[:, 0] == INT_FILL_VALUE,
-            edge_nodes_unique[:, 1] == INT_FILL_VALUE)
-
-        # all edge nodes that do not contain a fill value
-        non_fill_value_mask = np.logical_not(fill_value_mask)
-        edge_nodes_unique = edge_nodes_unique[non_fill_value_mask]
-
-        # Update inverse_indices accordingly
-        indices_to_update = np.where(fill_value_mask)[0]
-
-        remove_mask = np.isin(inverse_indices, indices_to_update)
-        inverse_indices[remove_mask] = INT_FILL_VALUE
-
-        # Compute the indices where inverse_indices exceeds the values in indices_to_update
-        indexes = np.searchsorted(indices_to_update,
-                                  inverse_indices,
-                                  side='right')
-        # subtract the corresponding indexes from `inverse_indices`
-        for i in range(len(inverse_indices)):
-            if inverse_indices[i] != INT_FILL_VALUE:
-                inverse_indices[i] -= indexes[i]
-
-        # add Mesh2_edge_nodes to internal dataset
-        self._ds['Mesh2_edge_nodes'] = xr.DataArray(
-            edge_nodes_unique,
-            dims=["nMesh2_edge", "Two"],
-            attrs={
-                "cf_role":
-                    "edge_node_connectivity",
-                "_FillValue":
-                    INT_FILL_VALUE,
-                "long_name":
-                    "Maps every edge to the two nodes that it connects",
-                "start_index":
-                    INT_DTYPE(0),
-                "inverse_indices":
-                    inverse_indices,
-                "fill_value_mask":
-                    fill_value_mask
-            })
-
-    def _build_face_edges_connectivity(self):
-        """Constructs the UGRID connectivity variable (``Mesh2_face_edges``)
-        and stores it within the internal (``Grid._ds``) and through the
-        attribute (``Grid.Mesh2_face_edges``)."""
-        if ("Mesh2_edge_nodes" not in self._ds or
-                "inverse_indices" not in self._ds['Mesh2_edge_nodes'].attrs):
-            self._build_edge_node_connectivity(repopulate=True)
-
-        inverse_indices = self._ds['Mesh2_edge_nodes'].inverse_indices
-        inverse_indices = inverse_indices.reshape(self.nMesh2_face,
-                                                  self.nMaxMesh2_face_nodes)
-
-        self._ds["Mesh2_face_edges"] = xr.DataArray(
-            data=inverse_indices,
-            dims=["nMesh2_face", "nMaxMesh2_face_edges"],
-            attrs={
-                "cf_role":
-                    "face_edges_connectivity",
-                "start_index":
-                    INT_DTYPE(0),
-                "long_name":
-                    "Maps every edge to the two nodes that it connects",
-            })
-
-    def _populate_cartesian_xyz_coord(self):
-        """A helper function that populates the xyz attribute in
-        UXarray.Grid._ds. This function is called when we need to use the
-        cartesian coordinates for each node to do the calculation but the input
-        data only has the "Mesh2_node_x" and "Mesh2_node_y" in degree.
-
-        Note
-        ----
-        In the UXarray, we abide the UGRID convention and make sure the following attributes will always have its
-        corresponding units as stated below:
-
-        Mesh2_node_x
-         unit:  "degree_east" for longitude
-        Mesh2_node_y
-         unit:  "degrees_north" for latitude
-        Mesh2_node_z
-         unit:  "m"
-        Mesh2_node_cart_x
-         unit:  "m"
-        Mesh2_node_cart_y
-         unit:  "m"
-        Mesh2_node_cart_z
-         unit:  "m"
-        """
-
-        # Check if the cartesian coordinates are already populated
-        if "Mesh2_node_cart_x" in self._ds.keys():
-            return
-
-        # check for units and create Mesh2_node_cart_x/y/z set to self._ds
-        nodes_lon_rad = np.deg2rad(self.Mesh2_node_x.values)
-        nodes_lat_rad = np.deg2rad(self.Mesh2_node_y.values)
-        nodes_rad = np.stack((nodes_lon_rad, nodes_lat_rad), axis=1)
-        nodes_cart = np.asarray(
-            list(map(node_lonlat_rad_to_xyz, list(nodes_rad))))
-
-        self._ds["Mesh2_node_cart_x"] = xr.DataArray(
-            data=nodes_cart[:, 0],
-            dims=["nMesh2_node"],
-            attrs={
-                "standard_name": "cartesian x",
-                "units": "m",
-            })
-        self._ds["Mesh2_node_cart_y"] = xr.DataArray(
-            data=nodes_cart[:, 1],
-            dims=["nMesh2_node"],
-            attrs={
-                "standard_name": "cartesian y",
-                "units": "m",
-            })
-        self._ds["Mesh2_node_cart_z"] = xr.DataArray(
-            data=nodes_cart[:, 2],
-            dims=["nMesh2_node"],
-            attrs={
-                "standard_name": "cartesian z",
-                "units": "m",
-            })
-
-    def _populate_lonlat_coord(self):
-        """Helper function that populates the longitude and latitude and store
-        it into the Mesh2_node_x and Mesh2_node_y. This is called when the
-        input data has "Mesh2_node_x", "Mesh2_node_y", "Mesh2_node_z" in
-        meters. Since we want "Mesh2_node_x" and "Mesh2_node_y" always have the
-        "degree" units. For more details, please read the following.
-
-        Raises
-        ------
-            RuntimeError
-                Mesh2_node_x/y/z are not represented in the cartesian format with the unit 'm'/'meters' when calling this function"
-
-        Note
-        ----
-        In the UXarray, we abide the UGRID convention and make sure the following attributes will always have its
-        corresponding units as stated below:
-
-        Mesh2_node_x
-         unit:  "degree_east" for longitude
-        Mesh2_node_y
-         unit:  "degrees_north" for latitude
-        Mesh2_node_z
-         unit:  "m"
-        Mesh2_node_cart_x
-         unit:  "m"
-        Mesh2_node_cart_y
-         unit:  "m"
-        Mesh2_node_cart_z
-         unit:  "m"
-        """
-
-        # Check if the "Mesh2_node_x" is already in longitude
-        if "degree" in self._ds.Mesh2_node_x.units:
-            return
-
-        # Check if the input Mesh2_node_xyz" are represented in the cartesian format with the unit "m"
-        if ("m" not in self._ds.Mesh2_node_x.units) or ("m" not in self._ds.Mesh2_node_y.units) \
-                or ("m" not in self._ds.Mesh2_node_z.units):
-            raise RuntimeError(
-                "Expected: Mesh2_node_x/y/z should be represented in the cartesian format with the "
-                "unit 'm' when calling this function")
-
-        # Put the cartesian coordinates inside the proper data structure
-        self._ds["Mesh2_node_cart_x"] = xr.DataArray(
-            data=self._ds["Mesh2_node_x"].values)
-        self._ds["Mesh2_node_cart_y"] = xr.DataArray(
-            data=self._ds["Mesh2_node_y"].values)
-        self._ds["Mesh2_node_cart_z"] = xr.DataArray(
-            data=self._ds["Mesh2_node_z"].values)
-
-        # convert the input cartesian values into the longitude latitude degree
-        nodes_cart = np.stack(
-            (self._ds["Mesh2_node_x"].values, self._ds["Mesh2_node_y"].values,
-             self._ds["Mesh2_node_z"].values),
-            axis=1).tolist()
-        nodes_rad = list(map(node_xyz_to_lonlat_rad, nodes_cart))
-        nodes_degree = np.rad2deg(nodes_rad)
-        self._ds["Mesh2_node_x"] = xr.DataArray(
-            data=nodes_degree[:, 0],
-            dims=["nMesh2_node"],
-            attrs={
-                "standard_name": "longitude",
-                "long_name": "longitude of mesh nodes",
-                "units": "degrees_east",
-            })
-        self._ds["Mesh2_node_y"] = xr.DataArray(
-            data=nodes_degree[:, 1],
-            dims=["nMesh2_node"],
-            attrs={
-                "standard_name": "latitude",
-                "long_name": "latitude of mesh nodes",
-                "units": "degrees_north",
-            })
-
-    def _build_nNodes_per_face(self):
-        """Constructs ``nNodes_per_face``, which contains the number of non-
-        fill-value nodes for each face in ``Mesh2_face_nodes``"""
-
-        # padding to shape [nMesh2_face, nMaxMesh2_face_nodes + 1]
-        closed = np.ones((self.nMesh2_face, self.nMaxMesh2_face_nodes + 1),
-                         dtype=INT_DTYPE) * INT_FILL_VALUE
-
-        closed[:, :-1] = self.Mesh2_face_nodes.copy()
-
-        nNodes_per_face = np.argmax(closed == INT_FILL_VALUE, axis=1)
-
-        # add to internal dataset
-        self._ds["nNodes_per_face"] = xr.DataArray(
-            data=nNodes_per_face,
-            dims=["nMesh2_face"],
-            attrs={"long_name": "number of non-fill value nodes for each face"})
diff --git a/uxarray/io/_exodus.py b/uxarray/io/_exodus.py
index fee79b74..5b83e4dd 100644
--- a/uxarray/io/_exodus.py
+++ b/uxarray/io/_exodus.py
@@ -3,8 +3,8 @@ import numpy as np
 from pathlib import PurePath
 from datetime import datetime
 
-from uxarray.utils.helpers import _replace_fill_values
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.grid.connectivity import _replace_fill_values
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 
 # Exodus Number is one-based.
diff --git a/uxarray/io/_mpas.py b/uxarray/io/_mpas.py
index a42f3b9c..fb0cb9cd 100644
--- a/uxarray/io/_mpas.py
+++ b/uxarray/io/_mpas.py
@@ -2,7 +2,7 @@ import xarray as xr
 import numpy as np
 import warnings
 
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 
 def _primal_to_ugrid(in_ds, out_ds):
diff --git a/uxarray/io/_scrip.py b/uxarray/io/_scrip.py
index d19f779f..a0503b2e 100644
--- a/uxarray/io/_scrip.py
+++ b/uxarray/io/_scrip.py
@@ -1,9 +1,9 @@
 import xarray as xr
 import numpy as np
 
-from uxarray.utils.helpers import grid_center_lat_lon, _replace_fill_values
-
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.grid.connectivity import _replace_fill_values
+from uxarray.grid.coordinates import grid_center_lat_lon
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 
 def _to_ugrid(in_ds, out_ds):
diff --git a/uxarray/io/_ugrid.py b/uxarray/io/_ugrid.py
index 42e298aa..8fda938d 100644
--- a/uxarray/io/_ugrid.py
+++ b/uxarray/io/_ugrid.py
@@ -1,7 +1,6 @@
-import xarray as xr
 import numpy as np
-from uxarray.utils.helpers import _replace_fill_values
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.grid.connectivity import _replace_fill_values
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 
 def _read_ugrid(xr_ds, var_names_dict):
@@ -113,3 +112,19 @@ def _standardize_fill_values(ds, var_names_dict):
             '_FillValue'] = INT_FILL_VALUE
 
     return ds
+
+
+def _is_ugrid(ds):
+    """Check mesh topology and dimension."""
+    standard_name = lambda v: v is not None
+    # getkeys_filter_by_attribute(filepath, attr_name, attr_val)
+    # return type KeysView
+    node_coords_dv = ds.filter_by_attrs(node_coordinates=standard_name)
+    face_conn_dv = ds.filter_by_attrs(face_node_connectivity=standard_name)
+    topo_dim_dv = ds.filter_by_attrs(topology_dimension=standard_name)
+    mesh_topo_dv = ds.filter_by_attrs(cf_role="mesh_topology")
+    if len(mesh_topo_dv) != 0 and len(topo_dim_dv) != 0 and len(
+            face_conn_dv) != 0 and len(node_coords_dv) != 0:
+        return True
+    else:
+        return False
diff --git a/uxarray/io/utils.py b/uxarray/io/utils.py
new file mode 100644
index 00000000..ce531748
--- /dev/null
+++ b/uxarray/io/utils.py
@@ -0,0 +1,86 @@
+from uxarray.io._ugrid import _is_ugrid
+
+
+def _parse_grid_type(dataset):
+    """Checks input and contents to determine grid type. Supports detection of
+    UGrid, SCRIP, Exodus and shape file.
+
+    Parameters
+    ----------
+    dataset : Xarray dataset
+       Xarray dataset of the grid
+
+    Returns
+    -------
+    mesh_type : str
+        File type of the file, ug, exo, scrip or shp
+
+    Raises
+    ------
+    RuntimeError
+            If invalid file type
+    ValueError
+        If file is not in UGRID format
+    """
+    # exodus with coord or coordx
+    if "coord" in dataset:
+        mesh_type = "exo"
+    elif "coordx" in dataset:
+        mesh_type = "exo"
+    # scrip with grid_center_lon
+    elif "grid_center_lon" in dataset:
+        mesh_type = "scrip"
+    # ugrid topology
+    elif _is_ugrid(dataset):
+        mesh_type = "ugrid"
+    elif "verticesOnCell" in dataset:
+        mesh_type = "mpas"
+    else:
+        raise RuntimeError(f"Could not recognize dataset format.")
+    return mesh_type
+
+    # check mesh topology and dimension
+    # try:
+    #     standard_name = lambda v: v is not None
+    #     # getkeys_filter_by_attribute(filepath, attr_name, attr_val)
+    #     # return type KeysView
+    #     ext_ds = xr.open_dataset(filepath, mask_and_scale=False)
+    #     node_coords_dv = ext_ds.filter_by_attrs(
+    #         node_coordinates=standard_name).keys()
+    #     face_conn_dv = ext_ds.filter_by_attrs(
+    #         face_node_connectivity=standard_name).keys()
+    #     topo_dim_dv = ext_ds.filter_by_attrs(
+    #         topology_dimension=standard_name).keys()
+    #     mesh_topo_dv = ext_ds.filter_by_attrs(cf_role="mesh_topology").keys()
+    #     if list(mesh_topo_dv)[0] != "" and list(topo_dim_dv)[0] != "" and list(
+    #             face_conn_dv)[0] != "" and list(node_coords_dv)[0] != "":
+    #         mesh_type = "ugrid"
+    #     else:
+    #         raise ValueError(
+    #             "cf_role is other than mesh_topology, the input NetCDF file is not UGRID format"
+    #         )
+    # except KeyError as e:
+    #     msg = str(e) + ': {}'.format(filepath)
+    # except (TypeError, AttributeError) as e:
+    #     msg = str(e) + ': {}'.format(filepath)
+    # except (RuntimeError, OSError) as e:
+    #     # check if this is a shp file
+    #     # we won't use xarray to load that file
+    #     if file_extension == ".shp":
+    #         mesh_type = "shp"
+    #     else:
+    #         msg = str(e) + ': {}'.format(filepath)
+    # except ValueError as e:
+    #     # check if this is a shp file
+    #     # we won't use xarray to load that file
+    #     if file_extension == ".shp":
+    #         mesh_type = "shp"
+    #     else:
+    #         msg = str(e) + ': {}'.format(filepath)
+    # finally:
+    #     if msg != "":
+    #         msg = "Unable to determine file type, mesh file not supported" + ': {}'.format(
+    #             filepath)
+    #         raise ValueError(msg)
+    #
+    # return mesh_type
diff --git a/uxarray/utils/__init__.py b/uxarray/utils/__init__.py
index e69de29b..914b45a1 100644
--- a/uxarray/utils/__init__.py
+++ b/uxarray/utils/__init__.py
@@ -0,0 +1,1 @@
+from .numba_settings import *
diff --git a/uxarray/utils/helpers.py b/uxarray/utils/helpers.py
deleted file mode 100644
index b471c1ba..00000000
--- a/uxarray/utils/helpers.py
+++ /dev/null
@@ -1,676 +0,0 @@
-import numpy as np
-import xarray as xr
-from pathlib import PurePath
-from .get_quadratureDG import get_gauss_quadratureDG, get_tri_quadratureDG
-from numba import njit, config
-from uxarray.utils.constants import ENABLE_JIT_CACHE, ENABLE_JIT
-import math
-
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
-
-config.DISABLE_JIT = not ENABLE_JIT
-
-
-def parse_grid_type(dataset):
-    """Checks input and contents to determine grid type. Supports detection of
-    UGrid, SCRIP, Exodus and shape file.
-
-    Parameters
-    ----------
-    dataset : Xarray dataset
-       Xarray dataset of the grid
-
-    Returns
-    -------
-    mesh_type : str
-        File type of the file, ug, exo, scrip or shp
-
-    Raises
-    ------
-    RuntimeError
-            If invalid file type
-    ValueError
-        If file is not in UGRID format
-    """
-    # exodus with coord or coordx
-    if "coord" in dataset:
-        mesh_type = "exo"
-    elif "coordx" in dataset:
-        mesh_type = "exo"
-    # scrip with grid_center_lon
-    elif "grid_center_lon" in dataset:
-        mesh_type = "scrip"
-    # ugrid topology
-    elif _is_ugrid(dataset):
-        mesh_type = "ugrid"
-    elif "verticesOnCell" in dataset:
-        mesh_type = "mpas"
-    else:
-        raise RuntimeError(f"Could not recognize dataset format.")
-    return mesh_type
-
-    # check mesh topology and dimension
-    try:
-        standard_name = lambda v: v is not None
-        # getkeys_filter_by_attribute(filepath, attr_name, attr_val)
-        # return type KeysView
-        ext_ds = xr.open_dataset(filepath, mask_and_scale=False)
-        node_coords_dv = ext_ds.filter_by_attrs(
-            node_coordinates=standard_name).keys()
-        face_conn_dv = ext_ds.filter_by_attrs(
-            face_node_connectivity=standard_name).keys()
-        topo_dim_dv = ext_ds.filter_by_attrs(
-            topology_dimension=standard_name).keys()
-        mesh_topo_dv = ext_ds.filter_by_attrs(cf_role="mesh_topology").keys()
-        if list(mesh_topo_dv)[0] != "" and list(topo_dim_dv)[0] != "" and list(
-                face_conn_dv)[0] != "" and list(node_coords_dv)[0] != "":
-            mesh_type = "ugrid"
-        else:
-            raise ValueError(
-                "cf_role is other than mesh_topology, the input NetCDF file is not UGRID format"
-            )
-    except KeyError as e:
-        msg = str(e) + ': {}'.format(filepath)
-    except (TypeError, AttributeError) as e:
-        msg = str(e) + ': {}'.format(filepath)
-    except (RuntimeError, OSError) as e:
-        # check if this is a shp file
-        # we won't use xarray to load that file
-        if file_extension == ".shp":
-            mesh_type = "shp"
-        else:
-            msg = str(e) + ': {}'.format(filepath)
-    except ValueError as e:
-        # check if this is a shp file
-        # we won't use xarray to load that file
-        if file_extension == ".shp":
-            mesh_type = "shp"
-        else:
-            msg = str(e) + ': {}'.format(filepath)
-    finally:
-        if msg != "":
-            msg = "Unable to determine file type, mesh file not supported" + ': {}'.format(
-                filepath)
-            raise ValueError(msg)
-
-    return mesh_type
-
-
-# Calculate the area of all faces.
-@njit(cache=ENABLE_JIT_CACHE)
-def calculate_face_area(x,
-                        y,
-                        z,
-                        quadrature_rule="gaussian",
-                        order=4,
-                        coords_type="spherical"):
-    """Calculate area of a face on sphere.
-
-    Parameters
-    ----------
-    x : list, required
-        x-coordinate of all the nodes forming the face
-
-    y : list, required
-        y-coordinate of all the nodes forming the face
-
-    z : list, required
-        z-coordinate of all the nodes forming the face
-
-    quadrature_rule : str, optional
-        triangular and Gaussian quadrature supported, expected values: "triangular" or "gaussian"
-
-    order: int, optional
-        Order of the quadrature rule. Default is 4.
-
-        Supported values:
-            - Gaussian Quadrature: 1 to 10
-            - Triangular: 1, 4, 8, 10 and 12
-
-    coords_type : str, optional
-        coordinate type, default is spherical, can be cartesian also.
-    """
-    area = 0.0  # set area to 0
-    order = order
-
-    if quadrature_rule == "gaussian":
-        dG, dW = get_gauss_quadratureDG(order)
-    elif quadrature_rule == "triangular":
-        dG, dW = get_tri_quadratureDG(order)
-    else:
-        raise ValueError(
-            "Invalid quadrature rule, specify gaussian or triangular")
-
-    num_nodes = len(x)
-
-    # num triangles is two less than the total number of nodes
-    num_triangles = num_nodes - 2
-
-    # Using tempestremap GridElements: https://github.com/ClimateGlobalChange/tempestremap/blob/master/src/GridElements.cpp
-    # loop through all sub-triangles of face
-    for j in range(0, num_triangles):
-        node1 = np.array([x[0], y[0], z[0]], dtype=np.float64)
-        node2 = np.array([x[j + 1], y[j + 1], z[j + 1]], dtype=np.float64)
-        node3 = np.array([x[j + 2], y[j + 2], z[j + 2]], dtype=np.float64)
-
-        if (coords_type == "spherical"):
-            node1 = np.array(
-                node_lonlat_rad_to_xyz([np.deg2rad(x[0]),
-                                        np.deg2rad(y[0])]))
-            node2 = np.array(
-                node_lonlat_rad_to_xyz(
-                    [np.deg2rad(x[j + 1]),
-                     np.deg2rad(y[j + 1])]))
-            node3 = np.array(
-                node_lonlat_rad_to_xyz(
-                    [np.deg2rad(x[j + 2]),
-                     np.deg2rad(y[j + 2])]))
-
-        for p in range(len(dW)):
-            if quadrature_rule == "gaussian":
-                for q in range(len(dW)):
-                    dA = dG[0][p]
-                    dB = dG[0][q]
-                    jacobian = calculate_spherical_triangle_jacobian(
-                        node1, node2, node3, dA, dB)
-                    area += dW[p] * dW[q] * jacobian
-            elif quadrature_rule == "triangular":
-                dA = dG[p][0]
-                dB = dG[p][1]
-                jacobian = calculate_spherical_triangle_jacobian_barycentric(
-                    node1, node2, node3, dA, dB)
-                area += dW[p] * jacobian
-
-    return area
-
-
-@njit(cache=ENABLE_JIT_CACHE)
-def get_all_face_area_from_coords(x,
-                                  y,
-                                  z,
-                                  face_nodes,
-                                  face_geometry,
-                                  dim,
-                                  quadrature_rule="triangular",
-                                  order=4,
-                                  coords_type="spherical"):
-    """Given coords, connectivity and other area calculation params, this
-    routine loop over all faces and return an numpy array with areas of each
-    face.
-
-    Parameters
-    ----------
-    x : ndarray, required
-        x-coordinate of all the nodes
-
-    y : ndarray, required
-        y-coordinate of all the nodes
-
-    z : ndarray, required
-        z-coordinate of all the nodes
-
-    face_nodes : 2D ndarray, required
-         node ids of each face
-
-    dim : int, required
-         dimension
-
-    quadrature_rule : str, optional
-        "triangular" or "gaussian". Defaults to triangular
-
-    order : int, optional
-        count or order for Gaussian or spherical resp. Defaults to 4 for spherical.
-
-    coords_type : str, optional
-        coordinate type, default is spherical, can be cartesian also.
-
-    Returns
-    -------
-    area of all faces : ndarray
-    """
-
-    n_face, n_max_face_nodes = face_nodes.shape
-
-    # set initial area of each face to 0
-    area = np.zeros(n_face)
-
-    for face_idx, max_nodes in enumerate(face_geometry):
-        face_x = x[face_nodes[face_idx, 0:max_nodes]]
-        face_y = y[face_nodes[face_idx, 0:max_nodes]]
-
-        # check if z dimension
-
-        if dim > 2:
-            face_z = z[face_nodes[face_idx, 0:max_nodes]]
-        else:
-            face_z = face_x * 0.0
-
-        # After getting all the nodes of a face assembled call the  cal. face area routine
-        face_area = calculate_face_area(face_x, face_y, face_z, quadrature_rule,
-                                        order, coords_type)
-        # store current face area
-        area[face_idx] = face_area
-
-    return area
-
-
-@njit(cache=ENABLE_JIT_CACHE)
-def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB):
-    """Calculate Jacobian of a spherical triangle. This is a helper function
-    for calculating face area.
-
-    Parameters
-    ----------
-    node1 : list, required
-        First node of the triangle
-
-    node2 : list, required
-        Second node of the triangle
-
-    node3 : list, required
-        Third node of the triangle
-
-    dA : float, required
-        quadrature point
-
-    dB : float, required
-        quadrature point
-
-    Returns
-    -------
-    jacobian : float
-    """
-    dF = np.array([
-        (1.0 - dB) * ((1.0 - dA) * node1[0] + dA * node2[0]) + dB * node3[0],
-        (1.0 - dB) * ((1.0 - dA) * node1[1] + dA * node2[1]) + dB * node3[1],
-        (1.0 - dB) * ((1.0 - dA) * node1[2] + dA * node2[2]) + dB * node3[2]
-    ])
-
-    dDaF = np.array([(1.0 - dB) * (node2[0] - node1[0]),
-                     (1.0 - dB) * (node2[1] - node1[1]),
-                     (1.0 - dB) * (node2[2] - node1[2])])
-
-    dDbF = np.array([
-        -(1.0 - dA) * node1[0] - dA * node2[0] + node3[0],
-        -(1.0 - dA) * node1[1] - dA * node2[1] + node3[1],
-        -(1.0 - dA) * node1[2] - dA * node2[2] + node3[2]
-    ])
-
-    dInvR = 1.0 / np.sqrt(dF[0] * dF[0] + dF[1] * dF[1] + dF[2] * dF[2])
-
-    dDaG = np.array([
-        dDaF[0] * (dF[1] * dF[1] + dF[2] * dF[2]) - dF[0] *
-        (dDaF[1] * dF[1] + dDaF[2] * dF[2]),
-        dDaF[1] * (dF[0] * dF[0] + dF[2] * dF[2]) - dF[1] *
-        (dDaF[0] * dF[0] + dDaF[2] * dF[2]),
-        dDaF[2] * (dF[0] * dF[0] + dF[1] * dF[1]) - dF[2] *
-        (dDaF[0] * dF[0] + dDaF[1] * dF[1])
-    ])
-
-    dDbG = np.array([
-        dDbF[0] * (dF[1] * dF[1] + dF[2] * dF[2]) - dF[0] *
-        (dDbF[1] * dF[1] + dDbF[2] * dF[2]),
-        dDbF[1] * (dF[0] * dF[0] + dF[2] * dF[2]) - dF[1] *
-        (dDbF[0] * dF[0] + dDbF[2] * dF[2]),
-        dDbF[2] * (dF[0] * dF[0] + dF[1] * dF[1]) - dF[2] *
-        (dDbF[0] * dF[0] + dDbF[1] * dF[1])
-    ])
-
-    dDenomTerm = dInvR * dInvR * dInvR
-
-    dDaG *= dDenomTerm
-    dDbG *= dDenomTerm
-
-    #  Cross product gives local Jacobian
-    nodeCross = np.cross(dDaG, dDbG)
-    dJacobian = np.sqrt(nodeCross[0] * nodeCross[0] +
-                        nodeCross[1] * nodeCross[1] +
-                        nodeCross[2] * nodeCross[2])
-
-    return dJacobian
-
-
-@njit(cache=ENABLE_JIT_CACHE)
-def calculate_spherical_triangle_jacobian_barycentric(node1, node2, node3, dA,
-                                                      dB):
-    """Calculate Jacobian of a spherical triangle. This is a helper function
-    for calculating face area.
-
-    Parameters
-    ----------
-    node1 : list, required
-        First node of the triangle
-
-    node2 : list, required
-        Second node of the triangle
-
-    node3 : list, required
-        Third node of the triangle
-
-    dA : float, required
-        first component of barycentric coordinates of quadrature point
-
-    dB : float, required
-        second component of barycentric coordinates of quadrature point
-
-    Returns
-    -------
-    jacobian : float
-    """
-
-    dF = np.array([
-        dA * node1[0] + dB * node2[0] + (1.0 - dA - dB) * node3[0],
-        dA * node1[1] + dB * node2[1] + (1.0 - dA - dB) * node3[1],
-        dA * node1[2] + dB * node2[2] + (1.0 - dA - dB) * node3[2]
-    ])
-
-    dDaF = np.array(
-        [node1[0] - node3[0], node1[1] - node3[1], node1[2] - node3[2]])
-
-    dDbF = np.array(
-        [node2[0] - node3[0], node2[1] - node3[1], node2[2] - node3[2]])
-
-    dInvR = 1.0 / np.sqrt(dF[0] * dF[0] + dF[1] * dF[1] + dF[2] * dF[2])
-
-    dDaG = np.array([
-        dDaF[0] * (dF[1] * dF[1] + dF[2] * dF[2]) - dF[0] *
-        (dDaF[1] * dF[1] + dDaF[2] * dF[2]),
-        dDaF[1] * (dF[0] * dF[0] + dF[2] * dF[2]) - dF[1] *
-        (dDaF[0] * dF[0] + dDaF[2] * dF[2]),
-        dDaF[2] * (dF[0] * dF[0] + dF[1] * dF[1]) - dF[2] *
-        (dDaF[0] * dF[0] + dDaF[1] * dF[1])
-    ])
-
-    dDbG = np.array([
-        dDbF[0] * (dF[1] * dF[1] + dF[2] * dF[2]) - dF[0] *
-        (dDbF[1] * dF[1] + dDbF[2] * dF[2]),
-        dDbF[1] * (dF[0] * dF[0] + dF[2] * dF[2]) - dF[1] *
-        (dDbF[0] * dF[0] + dDbF[2] * dF[2]),
-        dDbF[2] * (dF[0] * dF[0] + dF[1] * dF[1]) - dF[2] *
-        (dDbF[0] * dF[0] + dDbF[1] * dF[1])
-    ])
-
-    dDenomTerm = dInvR * dInvR * dInvR
-
-    dDaG *= dDenomTerm
-    dDbG *= dDenomTerm
-
-    #  Cross product gives local Jacobian
-    nodeCross = np.cross(dDaG, dDbG)
-    dJacobian = np.sqrt(nodeCross[0] * nodeCross[0] +
-                        nodeCross[1] * nodeCross[1] +
-                        nodeCross[2] * nodeCross[2])
-
-    return 0.5 * dJacobian
-
-
-def _is_ugrid(ds):
-    """Check mesh topology and dimension."""
-    standard_name = lambda v: v is not None
-    # getkeys_filter_by_attribute(filepath, attr_name, attr_val)
-    # return type KeysView
-    node_coords_dv = ds.filter_by_attrs(node_coordinates=standard_name)
-    face_conn_dv = ds.filter_by_attrs(face_node_connectivity=standard_name)
-    topo_dim_dv = ds.filter_by_attrs(topology_dimension=standard_name)
-    mesh_topo_dv = ds.filter_by_attrs(cf_role="mesh_topology")
-    if len(mesh_topo_dv) != 0 and len(topo_dim_dv) != 0 and len(
-            face_conn_dv) != 0 and len(node_coords_dv) != 0:
-        return True
-    else:
-        return False
-
-
-def grid_center_lat_lon(ds):
-    """Using scrip file variables ``grid_corner_lat`` and ``grid_corner_lon``,
-    calculates the ``grid_center_lat`` and ``grid_center_lon``.
-
-    Parameters
-    ----------
-    ds : xarray.Dataset
-        Dataset that contains ``grid_corner_lat`` and ``grid_corner_lon``
-        data variables
-
-    Returns
-    -------
-    center_lon : :class:`numpy.ndarray`
-        The calculated center longitudes of the grid box based on the corner
-        points
-    center_lat : :class:`numpy.ndarray`
-        The calculated center latitudes of the grid box based on the corner
-        points
-    """
-
-    # Calculate and create grid center lat/lon
-    scrip_corner_lon = ds['grid_corner_lon']
-    scrip_corner_lat = ds['grid_corner_lat']
-
-    # convert to radians
-    rad_corner_lon = np.deg2rad(scrip_corner_lon)
-    rad_corner_lat = np.deg2rad(scrip_corner_lat)
-
-    # get nodes per face
-    nodes_per_face = rad_corner_lat.shape[1]
-
-    # geographic center of each cell
-    x = np.sum(np.cos(rad_corner_lat) * np.cos(rad_corner_lon),
-               axis=1) / nodes_per_face
-    y = np.sum(np.cos(rad_corner_lat) * np.sin(rad_corner_lon),
-               axis=1) / nodes_per_face
-    z = np.sum(np.sin(rad_corner_lat), axis=1) / nodes_per_face
-
-    center_lon = np.rad2deg(np.arctan2(y, x))
-    center_lat = np.rad2deg(np.arctan2(z, np.sqrt(x**2 + y**2)))
-
-    # Make negative lons positive
-    center_lon[center_lon < 0] += 360
-
-    return center_lat, center_lon
-
-
-@njit(cache=ENABLE_JIT_CACHE)
-def node_lonlat_rad_to_xyz(node_coord):
-    """Helper function to Convert the node coordinate from 2D
-    longitude/latitude to normalized 3D xyz.
-
-    Parameters
-    ----------
-    node: float list
-        2D coordinates[longitude, latitude] in radiance
-
-    Returns
-    ----------
-    float list
-        the result array of the unit 3D coordinates [x, y, z] vector where :math:`x^2 + y^2 + z^2 = 1`
-
-    Raises
-    ----------
-    RuntimeError
-        The input array doesn't have the size of 3.
-    """
-    if len(node_coord) != 2:
-        raise RuntimeError(
-            "Input array should have a length of 2: [longitude, latitude]")
-    lon = node_coord[0]
-    lat = node_coord[1]
-    return [np.cos(lon) * np.cos(lat), np.sin(lon) * np.cos(lat), np.sin(lat)]
-
-
-@njit(cache=ENABLE_JIT_CACHE)
-def node_xyz_to_lonlat_rad(node_coord):
-    """Calculate the latitude and longitude in radiance for a node represented
-    in the [x, y, z] 3D Cartesian coordinates.
-
-    Parameters
-    ----------
-    node_coord: float list
-        3D Cartesian Coordinates [x, y, z] of the node
-
-    Returns
-    ----------
-    float list
-        the result array of longitude and latitude in radian [longitude_rad, latitude_rad]
-
-    Raises
-    ----------
-    RuntimeError
-        The input array doesn't have the size of 3.
-    """
-    if len(node_coord) != 3:
-        raise RuntimeError("Input array should have a length of 3: [x, y, z]")
-    reference_tolerance = 1.0e-12
-    [dx, dy, dz] = normalize_in_place(node_coord)
-    dx /= np.absolute(dx * dx + dy * dy + dz * dz)
-    dy /= np.absolute(dx * dx + dy * dy + dz * dz)
-    dz /= np.absolute(dx * dx + dy * dy + dz * dz)
-
-    if np.absolute(dz) < (1.0 - reference_tolerance):
-        d_lon_rad = math.atan2(dy, dx)
-        d_lat_rad = np.arcsin(dz)
-
-        if d_lon_rad < 0.0:
-            d_lon_rad += 2.0 * np.pi
-    elif dz > 0.0:
-        d_lon_rad = 0.0
-        d_lat_rad = 0.5 * np.pi
-    else:
-        d_lon_rad = 0.0
-        d_lat_rad = -0.5 * np.pi
-
-    return [d_lon_rad, d_lat_rad]
-
-
-@njit(cache=ENABLE_JIT_CACHE)
-def normalize_in_place(node):
-    """Helper function to project an arbitrary node in 3D coordinates [x, y, z]
-    on the unit sphere. It uses the `np.linalg.norm` internally to calculate
-    the magnitude.
-
-    Parameters
-    ----------
-    node: float list
-        3D Cartesian Coordinates [x, y, z]
-
-    Returns
-    ----------
-    float list
-        the result unit vector [x, y, z] where :math:`x^2 + y^2 + z^2 = 1`
-
-    Raises
-    ----------
-    RuntimeError
-        The input array doesn't have the size of 3.
-    """
-    if len(node) != 3:
-        raise RuntimeError("Input array should have a length of 3: [x, y, z]")
-
-    return 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, optional
-        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
-    if original_fill is not None and np.isnan(original_fill):
-        fill_val_idx = np.isnan(grid_var)
-    else:
-        fill_val_idx = grid_var == 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 fill value can be represented with current integer data type
-    if np.issubdtype(new_dtype, np.integer):
-        int_min = np.iinfo(grid_var.dtype).min
-        int_max = np.iinfo(grid_var.dtype).max
-        # ensure new_fill is in range [int_min, int_max]
-        if new_fill < int_min or new_fill > int_max:
-            raise ValueError(f'New fill value: {new_fill} not representable by'
-                             f' integer dtype: {grid_var.dtype}')
-
-    # ensure non-nan fill value can be represented with current float data type
-    elif np.issubdtype(new_dtype, np.floating) and not np.isnan(new_fill):
-        float_min = np.finfo(grid_var.dtype).min
-        float_max = np.finfo(grid_var.dtype).max
-        # ensure new_fill is in range [float_min, float_max]
-        if new_fill < float_min or new_fill > float_max:
-            raise ValueError(f'New fill value: {new_fill} not representable by'
-                             f' float dtype: {grid_var.dtype}')
-    else:
-        raise ValueError(f'Data type {grid_var.dtype} not supported'
-                         f'for grid variables')
-
-    # replace all zeros with a fill value
-    grid_var[fill_val_idx] = new_fill
-
-    return grid_var
-
-
-def close_face_nodes(Mesh2_face_nodes, nMesh2_face, nMaxMesh2_face_nodes):
-    """Closes (``Mesh2_face_nodes``) by inserting the first node index after
-    the last non-fill-value node.
-
-    Parameters
-    ----------
-    Mesh2_face_nodes : np.ndarray
-        Connectivity array for constructing a face from its nodes
-    nMesh2_face : constant
-        Number of faces
-    nMaxMesh2_face_nodes : constant
-        Max number of nodes that compose a face
-
-    Returns
-    ----------
-    closed : ndarray
-        Closed (padded) Mesh2_face_nodes
-
-    Example
-    ----------
-    Given face nodes with shape [2 x 5]
-        [0, 1, 2, 3, FILL_VALUE]
-        [4, 5, 6, 7, 8]
-    Pads them to the following with shape [2 x 6]
-        [0, 1, 2, 3, 0, FILL_VALUE]
-        [4, 5, 6, 7, 8, 4]
-    """
-
-    # padding to shape [nMesh2_face, nMaxMesh2_face_nodes + 1]
-    closed = np.ones((nMesh2_face, nMaxMesh2_face_nodes + 1),
-                     dtype=INT_DTYPE) * INT_FILL_VALUE
-
-    # set all non-paded values to original face nodee values
-    closed[:, :-1] = Mesh2_face_nodes.copy()
-
-    # instance of first fill value
-    first_fv_idx_2d = np.argmax(closed == INT_FILL_VALUE, axis=1)
-
-    # 2d to 1d index for np.put()
-    first_fv_idx_1d = first_fv_idx_2d + (
-        (nMaxMesh2_face_nodes + 1) * np.arange(0, nMesh2_face))
-
-    # column of first node values
-    first_node_value = Mesh2_face_nodes[:, 0].copy()
-
-    # insert first node column at occurrence of first fill value
-    np.put(closed.ravel(), first_fv_idx_1d, first_node_value)
-
-    return closed
diff --git a/uxarray/utils/numba_settings.py b/uxarray/utils/numba_settings.py
new file mode 100644
index 00000000..ea7e28ff
--- /dev/null
+++ b/uxarray/utils/numba_settings.py
@@ -0,0 +1,49 @@
+import uxarray.constants
+
+
+def enable_jit_cache():
+    """Allows Numba's JIT cache to be turned on.
+
+    This cache variable lets @njit cache the machine code generated
+    between runs, allowing for faster run times due to the fact that the
+    code doesn't need to regenerate the machine code every run time. Our
+    use case here was to study performance, in regular usage one might
+    never turn off caching as it will only help if frequently modifying
+    the code or because users have very limited disk space. The default
+    is on (True)
+    """
+    uxarray.constants.ENABLE_JIT_CACHE = True
+
+
+def disable_jit_cache():
+    """Allows Numba's JIT cache to be turned on off.
+
+    This cache variable lets @njit cache the machine code generated
+    between runs, allowing for faster run times due to the fact that the
+    code doesn't need to regenerate the machine code every run time. Our
+    use case here was to study performance, in regular usage one might
+    never turn off caching as it will only help if frequently modifying
+    the code or because users have very limited disk space. The default
+    is on (True)
+    """
+    uxarray.constants.ENABLE_JIT_CACHE = False
+
+
+def enable_jit():
+    """Allows Numba's JIT application to be turned on.
+
+    This lets users choose whether they want machine code to be
+    generated to speed up the performance of the code on large files.
+    The default is on (True)
+    """
+    uxarray.constants.ENABLE_JIT = True
+
+
+def disable_jit():
+    """Allows Numba's JIT application to be turned off.
+
+    This lets users choose whether they want machine code to be
+    generated to speed up the performance of the code on large files.
+    The default is on (True)
+    """
+    uxarray.constants.ENABLE_JIT = False


Test patch:
diff --git a/test/test_exodus.py b/test/test_exodus.py
index dda8c660..a07bbb65 100644
--- a/test/test_exodus.py
+++ b/test/test_exodus.py
@@ -5,7 +5,7 @@ from unittest import TestCase
 from pathlib import Path
 
 import uxarray as ux
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
diff --git a/test/test_grid.py b/test/test_grid.py
index 50bd3353..10ef2e70 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -8,6 +8,10 @@ from pathlib import Path
 
 import uxarray as ux
 
+from uxarray.grid.connectivity import _build_edge_node_connectivity, _build_face_edges_connectivity
+
+from uxarray.grid.coordinates import _populate_cartesian_xyz_coord, _populate_lonlat_coord
+
 try:
     import constants
 except ImportError:
@@ -429,7 +433,7 @@ class TestPopulateCoordinates(TestCase):
         verts_degree = np.stack((lon_deg, lat_deg), axis=1)
 
         vgrid = ux.open_grid(verts_degree, islatlon=False)
-        vgrid._populate_cartesian_xyz_coord()
+        _populate_cartesian_xyz_coord(vgrid)
 
         for i in range(0, vgrid.nMesh2_node):
             nt.assert_almost_equal(vgrid._ds["Mesh2_node_cart_x"].values[i],
@@ -470,7 +474,7 @@ class TestPopulateCoordinates(TestCase):
         verts_cart = np.stack((cart_x, cart_y, cart_z), axis=1)
 
         vgrid = ux.open_grid(verts_cart, islatlon=False)
-        vgrid._populate_lonlat_coord()
+        _populate_lonlat_coord(vgrid)
         # The connectivity in `__from_vert__()` will be formed in a reverse order
         lon_deg, lat_deg = zip(*reversed(list(zip(lon_deg, lat_deg))))
         for i in range(0, vgrid.nMesh2_node):
@@ -665,7 +669,7 @@ class TestConnectivity(TestCase):
         edge_nodes_expected = np.unique(edge_nodes_expected, axis=0)
 
         # construct edge nodes
-        mpas_grid_ux._build_edge_node_connectivity(repopulate=True)
+        _build_edge_node_connectivity(mpas_grid_ux, repopulate=True)
         edge_nodes_output = mpas_grid_ux._ds['Mesh2_edge_nodes'].values
 
         self.assertTrue(np.array_equal(edge_nodes_expected, edge_nodes_output))
@@ -688,7 +692,7 @@ class TestConnectivity(TestCase):
 
             mesh2_face_nodes = tgrid._ds["Mesh2_face_nodes"]
 
-            tgrid._build_face_edges_connectivity()
+            _build_face_edges_connectivity(tgrid)
             mesh2_face_edges = tgrid._ds.Mesh2_face_edges
             mesh2_edge_nodes = tgrid._ds.Mesh2_edge_nodes
 
@@ -723,7 +727,7 @@ class TestConnectivity(TestCase):
 
         mesh2_face_nodes = tgrid._ds["Mesh2_face_nodes"]
 
-        tgrid._build_face_edges_connectivity()
+        _build_face_edges_connectivity(tgrid)
         mesh2_face_edges = tgrid._ds.Mesh2_face_edges
         mesh2_edge_nodes = tgrid._ds.Mesh2_edge_nodes
 
@@ -746,7 +750,7 @@ class TestConnectivity(TestCase):
             self.f5_deg, self.f6_deg
         ]
         uds = ux.open_grid(verts)
-        uds._build_face_edges_connectivity()
+        _build_face_edges_connectivity(uds)
         n_face = len(uds._ds["Mesh2_face_edges"].values)
         n_node = uds.nMesh2_node
         n_edge = len(uds._ds["Mesh2_edge_nodes"].values)
diff --git a/test/test_helpers.py b/test/test_helpers.py
index fe2caee3..a14c52e0 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -9,8 +9,8 @@ from pathlib import Path
 
 import uxarray as ux
 
-from uxarray.utils.helpers import _replace_fill_values
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.grid.connectivity import _replace_fill_values
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 try:
     import constants
@@ -38,13 +38,8 @@ class TestIntegrate(TestCase):
         face_nodes = np.array([[0, 1, 2]]).astype(INT_DTYPE)
         face_dimension = np.array([3], dtype=INT_DTYPE)
 
-        area = ux.get_all_face_area_from_coords(x,
-                                                y,
-                                                z,
-                                                face_nodes,
-                                                face_dimension,
-                                                3,
-                                                coords_type="cartesian")
+        area = ux.grid.area.get_all_face_area_from_coords(
+            x, y, z, face_nodes, face_dimension, 3, coords_type="cartesian")
 
         nt.assert_almost_equal(area, constants.TRI_AREA, decimal=1)
 
@@ -56,20 +51,21 @@ class TestIntegrate(TestCase):
         y = np.array([-5.77350269e-01, 5.77350269e-01, 5.77350269e-01])
         z = np.array([-0.57735027, -0.57735027, -0.57735027])
 
-        area = ux.calculate_face_area(x, y, z, "gaussian", 5, "cartesian")
+        area = ux.grid.area.calculate_face_area(x, y, z, "gaussian", 5,
+                                                "cartesian")
 
         nt.assert_almost_equal(area, constants.TRI_AREA, decimal=3)
 
     def test_quadrature(self):
         order = 1
-        dG, dW = ux.get_tri_quadratureDG(order)
+        dG, dW = ux.grid.area.get_tri_quadratureDG(order)
         G = np.array([[0.33333333, 0.33333333, 0.33333333]])
         W = np.array([1.0])
 
         np.testing.assert_array_almost_equal(G, dG)
         np.testing.assert_array_almost_equal(W, dW)
 
-        dG, dW = ux.get_gauss_quadratureDG(order)
+        dG, dW = ux.grid.area.get_gauss_quadratureDG(order)
 
         G = np.array([[0.5]])
         W = np.array([1.0])
@@ -90,7 +86,7 @@ class TestGridCenter(TestCase):
         scrip_center_lat = ds_scrip_CSne8['grid_center_lat']
 
         # Calculate the center_lat/lon using same dataset's corner_lat/lon
-        calc_center = ux.grid_center_lat_lon(ds_scrip_CSne8)
+        calc_center = ux.grid.coordinates.grid_center_lat_lon(ds_scrip_CSne8)
         calc_lat = calc_center[0]
         calc_lon = calc_center[1]
 
@@ -102,7 +98,7 @@ class TestGridCenter(TestCase):
 class TestCoordinatesConversion(TestCase):
 
     def test_normalize_in_place(self):
-        [x, y, z] = ux.utils.helpers.normalize_in_place(
+        [x, y, z] = ux.grid.coordinates.normalize_in_place(
             [random.random(), random.random(),
              random.random()])
 
@@ -110,15 +106,15 @@ class TestCoordinatesConversion(TestCase):
                              err_tolerance)
 
     def test_node_xyz_to_lonlat_rad(self):
-        [x, y, z] = ux.utils.helpers.normalize_in_place([
+        [x, y, z] = ux.grid.coordinates.normalize_in_place([
             random.uniform(-1, 1),
             random.uniform(-1, 1),
             random.uniform(-1, 1)
         ])
 
-        [lon, lat] = ux.utils.helpers.node_xyz_to_lonlat_rad([x, y, z])
+        [lon, lat] = ux.grid.coordinates.node_xyz_to_lonlat_rad([x, y, z])
         [new_x, new_y,
-         new_z] = ux.utils.helpers.node_lonlat_rad_to_xyz([lon, lat])
+         new_z] = ux.grid.coordinates.node_lonlat_rad_to_xyz([lon, lat])
 
         self.assertLessEqual(np.absolute(new_x - x), err_tolerance)
         self.assertLessEqual(np.absolute(new_y - y), err_tolerance)
@@ -130,9 +126,10 @@ class TestCoordinatesConversion(TestCase):
             random.uniform(-0.5 * np.pi, 0.5 * np.pi)
         ]
 
-        [x, y, z] = ux.utils.helpers.node_lonlat_rad_to_xyz([lon, lat])
+        [x, y, z] = ux.grid.coordinates.node_lonlat_rad_to_xyz([lon, lat])
 
-        [new_lon, new_lat] = ux.utils.helpers.node_xyz_to_lonlat_rad([x, y, z])
+        [new_lon,
+         new_lat] = ux.grid.coordinates.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)
diff --git a/test/test_mpas.py b/test/test_mpas.py
index d95dd363..2c401d6b 100644
--- a/test/test_mpas.py
+++ b/test/test_mpas.py
@@ -7,7 +7,7 @@ import numpy as np
 import os
 from pathlib import Path
 
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
diff --git a/test/test_scrip.py b/test/test_scrip.py
index b9fd79ac..ab48296f 100644
--- a/test/test_scrip.py
+++ b/test/test_scrip.py
@@ -1,126 +1,103 @@
-import uxarray as ux
+import os
 import xarray as xr
+
 from unittest import TestCase
-import numpy as np
+from pathlib import Path
+import warnings
 import numpy.testing as nt
 
-import os
-from pathlib import Path
+import uxarray as ux
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
-from uxarray.io._scrip import _read_scrip, _encode_scrip
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+try:
+    import constants
+except ImportError:
+    from . import constants
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
-gridfile_ne30 = current_path / 'meshfiles' / "ugrid" / "outCSne30" / 'outCSne30.ug'
-gridfile_ne8 = current_path / 'meshfiles' / "scrip" / "outCSne8" / 'outCSne8.nc'
+gridfile_ne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
+gridfile_RLL1deg = current_path / "meshfiles" / "ugrid" / "outRLL1deg" / "outRLL1deg.ug"
+gridfile_RLL10deg_ne4 = current_path / "meshfiles" / "ugrid" / "ov_RLL10deg_CSne4" / "ov_RLL10deg_CSne4.ug"
 
-ds_ne30 = xr.open_dataset(gridfile_ne30, decode_times=False,
-                          engine='netcdf4')  # mesh2_node_x/y
-ds_ne8 = xr.open_dataset(gridfile_ne8, decode_times=False,
-                         engine='netcdf4')  # grid_corner_lat/lon
+gridfile_exo_ne8 = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g"
 
 
-class TestScrip(TestCase):
+class TestUgrid(TestCase):
 
-    def test_exception_nonSCRIP(self):
-        """Checks that exception is raised if non-SCRIP formatted file is
-        passed to function."""
+    def test_read_ugrid(self):
+        """Reads a ugrid file."""\
 
-        self.assertRaises(TypeError, _read_scrip(ds_ne30))
+        uxgrid_ne30 = ux.open_grid(str(gridfile_ne30))
+        uxgrid_RLL1deg = ux.open_grid(str(gridfile_RLL1deg))
+        uxgrid_RLL10deg_ne4 = ux.open_grid(str(gridfile_RLL10deg_ne4))
 
-    def test_scrip_is_not_ugrid(self):
-        """tests that function has correctly created a ugrid function and no
-        longer uses SCRIP variable names (grid_corner_lat), the function will
-        raise an exception."""
-        new_ds = _read_scrip(ds_ne8)
+        nt.assert_equal(uxgrid_ne30.Mesh2_node_x.size,
+                        constants.NNODES_outCSne30)
+        nt.assert_equal(uxgrid_RLL1deg.Mesh2_node_x.size,
+                        constants.NNODES_outRLL1deg)
+        nt.assert_equal(uxgrid_RLL10deg_ne4.Mesh2_node_x.size,
+                        constants.NNODES_ov_RLL10deg_CSne4)
 
-        assert ds_ne8['grid_corner_lat'].any()
+    def test_read_ugrid_opendap(self):
+        """Read an ugrid model from an OPeNDAP URL."""
 
-        with self.assertRaises(KeyError):
-            new_ds['grid_corner_lat']
+        try:
+            # make sure we can read the ugrid file from the OPeNDAP URL
+            url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc"
+            uxgrid_url = ux.open_grid(url, drop_variables="siglay")
 
-    def test_encode_scrip(self):
-        """Read a UGRID dataset and encode that as a SCRIP format.
+        except OSError:
+            # print warning and pass if we can't connect to the OPeNDAP server
+            warnings.warn(f'Could not connect to OPeNDAP server: {url}')
+            pass
 
-        Look for specific variable names in the returned dataset to
-        determine encoding was successful.
-        """
+        else:
+
+            assert isinstance(getattr(uxgrid_url, "Mesh2_node_x"), xr.DataArray)
+            assert isinstance(getattr(uxgrid_url, "Mesh2_node_y"), xr.DataArray)
+            assert isinstance(getattr(uxgrid_url, "Mesh2_face_nodes"),
+                              xr.DataArray)
+
+    def test_encode_ugrid(self):
+        """Read an Exodus dataset and encode that as a UGRID format."""
 
-        # Create UGRID from SCRIP file
-        scrip_in_ds = _read_scrip(ds_ne8)
-        new_path = current_path / "meshfiles" / "scrip_to_ugrid.ug"
-        scrip_in_ds.to_netcdf(str(new_path))  # Save as new file
-
-        # Use xarray open_dataset, create a uxarray grid object to then create SCRIP file from new UGRID file
-        ugrid_out = ux.open_grid(new_path)
-
-        scrip_encode_ds = _encode_scrip(ugrid_out.Mesh2_face_nodes,
-                                        ugrid_out.Mesh2_node_x,
-                                        ugrid_out.Mesh2_node_y,
-                                        ugrid_out.face_areas)
-
-        # Test newly created SCRIP is same as original SCRIP
-        nt.assert_array_almost_equal(scrip_encode_ds['grid_corner_lat'],
-                                     ds_ne8['grid_corner_lat'])
-        nt.assert_array_almost_equal(scrip_encode_ds['grid_corner_lon'],
-                                     ds_ne8['grid_corner_lon'])
-
-        # Tests that calculated center lat/lon values are equivalent to original
-        nt.assert_array_almost_equal(scrip_encode_ds['grid_center_lon'],
-                                     ds_ne8['grid_center_lon'])
-        nt.assert_array_almost_equal(scrip_encode_ds['grid_center_lat'],
-                                     ds_ne8['grid_center_lat'])
-
-        # Tests that calculated face area values are equivalent to original
-        nt.assert_array_almost_equal(scrip_encode_ds['grid_area'],
-                                     ds_ne8['grid_area'])
-
-        # Tests that calculated grid imask values are equivalent to original
-        nt.assert_array_almost_equal(scrip_encode_ds['grid_imask'],
-                                     ds_ne8['grid_imask'])
-
-        # Test that "mesh" variables are not in new file
-        with self.assertRaises(KeyError):
-            assert scrip_encode_ds['Mesh2_node_x']
-            assert scrip_encode_ds['Mesh2_node_y']
-
-    def test_scrip_variable_names(self):
-        """Tests that returned dataset from writer function has all required
-        SCRIP variables."""
-        uxds_ne30 = ux.open_grid(gridfile_ne30)
-        scrip30 = _encode_scrip(uxds_ne30.Mesh2_face_nodes,
-                                uxds_ne30.Mesh2_node_x, uxds_ne30.Mesh2_node_y,
-                                uxds_ne30.face_areas)
-
-        # List of relevant variable names for a scrip file
-        var_list = [
-            'grid_corner_lat', 'grid_dims', 'grid_imask', 'grid_area',
-            'grid_center_lon'
-        ]
-
-        for i in range(len(var_list) - 1):
-            assert scrip30[var_list[i]].any()
-
-    def test_ugrid_variable_names(self):
-        """Tests that returned dataset uses UGRID compliant variables."""
-        mesh08 = _read_scrip(ds_ne8)
-
-        # Create a flattened and unique array for comparisons
-        corner_lon = ds_ne8['grid_corner_lon'].values
-        corner_lon = corner_lon.flatten()
-        strip_lon = np.unique(corner_lon)
-
-        assert strip_lon.all() == mesh08['Mesh2_node_x'].all()
+        ux_grid = ux.open_grid(gridfile_exo_ne8)
+        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."""
 
-        ux_grid_01 = ux.open_grid(gridfile_ne30)
-        ux_grid_02 = ux.open_grid(gridfile_ne8)
+        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"
 
-        grids = [ux_grid_01, ux_grid_02]
-        for grid in grids:
+        ux_grid1 = ux.open_grid(ug_filename1)
+        ux_grid2 = ux.open_grid(ug_filename2)
+        ux_grid3 = ux.open_grid(ug_filename3)
+
+        # check for correct dtype and fill value
+        grids_with_fill = [ux_grid2]
+        for grid in grids_with_fill:
             assert grid.Mesh2_face_nodes.dtype == INT_DTYPE
             assert grid.Mesh2_face_nodes._FillValue == INT_FILL_VALUE
+            assert INT_FILL_VALUE in grid.Mesh2_face_nodes.values
+
+        grids_without_fill = [ux_grid1, ux_grid3]
+        for grid in grids_without_fill:
+            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" / "outRLL1deg" / "outRLL1deg.ug"
+        ux_grid = ux.open_grid(ug_filename)
+
+        assert ux_grid.Mesh2_face_nodes.dtype == INT_DTYPE
+        assert ux_grid.Mesh2_face_nodes._FillValue == INT_FILL_VALUE
+        assert INT_FILL_VALUE in ux_grid.Mesh2_face_nodes.values
diff --git a/test/test_ugrid.py b/test/test_ugrid.py
index c1c9b660..ab48296f 100644
--- a/test/test_ugrid.py
+++ b/test/test_ugrid.py
@@ -7,7 +7,7 @@ import warnings
 import numpy.testing as nt
 
 import uxarray as ux
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
 try:
     import constants