Repository: UXARRAY/uxarray

Issue:
Use `grid_var_names` everywhere
**Use `grid_var_names` everywhere**
UXARRY uses `grid_var_names` to map mesh topology variable names, but it is not used for all variables currently. For example, `Mesh2_edge_nodes` is not referred to through `grid_var_names`. Not using mapped names becomes an issue when a dataset uses a different topology information variable name from default `Mesh2`.

**Describe the solution you'd like**
Use `grid_var_names` for all topology variables.

Discussion:
erogluorhan: > good work!
> 
> let's merge this for other PR's to go through in this iteration.

This iteration is almost over, and tomorrow will be a release day for UXarray. I'd like to give this PR a review. We can always get a second release if we want in the same calendar month though.
philipc2: > > good work!
> > 
> > let's merge this for other PR's to go through in this iteration.
> 
> This iteration is almost over, and tomorrow will be a release day for UXarray. I'd like to give this PR a review. We can always get a second release if we want in the same calendar month though.

That's good with me! There's no rush and I agree that a second, smaller, release is definitely doable. 
philipc2: @erogluorhan 

Should be ready for a second pass through. I believe I've addressed most of the concerns, with the major change being the following:

Each grid now has a `_source_dims_dict`, which provides the following mapping {"Original Dimension Name" : "UGRID Dimension Name"}, instead of storing all the connectivity and coordinates.

This is now used in our API functions when opening a dataset to remap the original dimension names to their UGRID equivalent, if they were parsed. It can be specified using the "remap_dims" param.

Let me know what you think of this design, and this actually addresses #363 directly now 
erogluorhan: > Each grid now has a `_source_dims_dict`, which provides the following mapping {"Original Dimension Name" : "UGRID Dimension Name"}, instead of storing all the connectivity and coordinates.

Shouldn't it be {"Our standardized UGRID dim name" : "Original Dimension Name from the source"}. Maybe you already meant that?


philipc2: > > Each grid now has a `_source_dims_dict`, which provides the following mapping {"Original Dimension Name" : "UGRID Dimension Name"}, instead of storing all the connectivity and coordinates.
> 
> Shouldn't it be {"Our standardized UGRID dim name" : "Original Dimension Name from the source"}. Maybe you already meant that?

This was originally the case, however since we want to rename the dimensions within our `UxDataset`, we need the mapping described in the comment.

https://docs.xarray.dev/en/stable/generated/xarray.Dataset.rename_dims.html
review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/UXARRAY/uxarray/pull/395"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>
philipc2: @aaronzedwick @rajeeja 

Have either of you encountered this error before with Numba? After these changes, the notebook is starting to fail

![image](https://github.com/UXARRAY/uxarray/assets/67855069/65079fb0-236a-4c08-9ec8-5dd6ddd46c72)

Interestingly, running the same code snippet in testing doesn't cause the issue. Any thoughts? 
aaronzedwick: > @aaronzedwick @rajeeja
> 
> Have either of you encountered this error before with Numba? After these changes, the notebook is starting to fail
> 
> ![image](https://user-images.githubusercontent.com/67855069/268337550-65079fb0-236a-4c08-9ec8-5dd6ddd46c72.png)
> 
> Interestingly, running the same code snippet in testing doesn't cause the issue. Any thoughts?

Interesting. I have never encountered this specific error before. It looks like your last commits removed stuff from the notebook, so maybe that caused it? If it works inside your IDE it seems likely it could be something is wrong with the notebook, although looking at it I can’t see what. Maybe it got corrupted like Hongyu’s?
philipc2: > > @aaronzedwick @rajeeja
> > Have either of you encountered this error before with Numba? After these changes, the notebook is starting to fail
> > ![image](https://user-images.githubusercontent.com/67855069/268337550-65079fb0-236a-4c08-9ec8-5dd6ddd46c72.png)
> > Interestingly, running the same code snippet in testing doesn't cause the issue. Any thoughts?
> 
> Interesting. I have never encountered this specific error before. It looks like your last commits removed stuff from the notebook, so maybe that caused it? If it works inside your IDE it seems likely it could be something is wrong with the notebook, although looking at it I can’t see what. Maybe it got corrupted like Hongyu’s?

It was failing prior to updating the notebook, so I don't believe that's the issue. 


philipc2: > > @aaronzedwick @rajeeja
> > Have either of you encountered this error before with Numba? After these changes, the notebook is starting to fail
> > ![image](https://user-images.githubusercontent.com/67855069/268337550-65079fb0-236a-4c08-9ec8-5dd6ddd46c72.png)
> > Interestingly, running the same code snippet in testing doesn't cause the issue. Any thoughts?
> 
> Interesting. I have never encountered this specific error before. It looks like your last commits removed stuff from the notebook, so maybe that caused it? If it works inside your IDE it seems likely it could be something is wrong with the notebook, although looking at it I can’t see what. Maybe it got corrupted like Hongyu’s?

It's resolved now. I had to update the notebook to use use `ux.open_grid` instead of using the old `Grid` constructor approach to opening the grid. 
philipc2: > and are you planning to "always remap"? Should I expect that or approve without it?

I've made the change to always remap the dims. 
erogluorhan: > > and are you planning to "always remap"? Should I expect that or approve without it?
> 
> I've made the change to always remap the dims.

Can't see these changes
philipc2: > > > and are you planning to "always remap"? Should I expect that or approve without it?
> > 
> > 
> > I've made the change to always remap the dims.
> 
> Can't see these changes

Whoops! Forgot to push them.

Patch:
diff --git a/docs/examples/003-area-calc.ipynb b/docs/examples/003-area-calc.ipynb
index 7cf43601..a487467f 100644
--- a/docs/examples/003-area-calc.ipynb
+++ b/docs/examples/003-area-calc.ipynb
@@ -36,8 +36,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.775400Z",
-     "start_time": "2023-06-09T10:04:28.705400Z"
+     "start_time": "2023-09-16T09:11:10.389453Z",
+     "end_time": "2023-09-16T09:11:11.551466Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -56,8 +56,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.817401Z",
-     "start_time": "2023-06-09T10:04:28.718401Z"
+     "start_time": "2023-09-16T09:11:11.554467Z",
+     "end_time": "2023-09-16T09:11:11.616484Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -75,7 +75,7 @@
     "# Load grid files and create UXarray Grid objects\n",
     "ugrid_ds = xr.open_dataset(ugrid_path)\n",
     "\n",
-    "ugrid = ux.Grid(ugrid_ds)\n",
+    "ugrid = ux.open_grid(ugrid_ds)\n",
     "ugrid"
    ]
   },
@@ -92,8 +92,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.820400Z",
-     "start_time": "2023-06-09T10:04:28.732401Z"
+     "start_time": "2023-09-16T09:11:11.618485Z",
+     "end_time": "2023-09-16T09:11:17.203518Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -127,8 +127,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.821401Z",
-     "start_time": "2023-06-09T10:04:28.777401Z"
+     "start_time": "2023-09-16T09:11:17.219517Z",
+     "end_time": "2023-09-16T09:11:17.222938Z"
     }
    },
    "outputs": [],
@@ -165,8 +165,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.827402Z",
-     "start_time": "2023-06-09T10:04:28.810401Z"
+     "start_time": "2023-09-16T09:11:17.223939Z",
+     "end_time": "2023-09-16T09:11:17.228708Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -190,8 +190,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.904401Z",
-     "start_time": "2023-06-09T10:04:28.896401Z"
+     "start_time": "2023-09-16T09:11:17.227709Z",
+     "end_time": "2023-09-16T09:11:17.293300Z"
     }
    },
    "outputs": [],
@@ -213,8 +213,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.950400Z",
-     "start_time": "2023-06-09T10:04:28.904401Z"
+     "start_time": "2023-09-16T09:11:17.293300Z",
+     "end_time": "2023-09-16T09:11:17.297206Z"
     }
    },
    "outputs": [],
@@ -250,8 +250,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.958400Z",
-     "start_time": "2023-06-09T10:04:28.921401Z"
+     "start_time": "2023-09-16T09:11:17.298206Z",
+     "end_time": "2023-09-16T09:11:17.307131Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -265,10 +265,7 @@
     "         [-0.57735027, 5.77350269e-01, -0.57735027]]]\n",
     "\n",
     "# load our vertices into a UXarray Grid object\n",
-    "vgrid = ux.Grid(verts,\n",
-    "                vertices=True,\n",
-    "                islatlon=False,\n",
-    "                concave=False)\n",
+    "vgrid = ux.open_grid(verts, islatlon=False)\n",
     "\n",
     "vgrid"
    ]
@@ -278,8 +275,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:28.967401Z",
-     "start_time": "2023-06-09T10:04:28.951402Z"
+     "start_time": "2023-09-16T09:11:17.308131Z",
+     "end_time": "2023-09-16T09:11:18.468328Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -303,8 +300,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:29.014400Z",
-     "start_time": "2023-06-09T10:04:28.967401Z"
+     "start_time": "2023-09-16T09:12:05.404404Z",
+     "end_time": "2023-09-16T09:12:05.412402Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -313,10 +310,10 @@
    },
    "outputs": [],
    "source": [
-    "# Set correct units for the x and y coordinates\n",
-    "vgrid.Mesh2_node_x.attrs[\"units\"] = \"km\"\n",
-    "vgrid.Mesh2_node_y.attrs[\"units\"] = \"km\"\n",
-    "vgrid._Mesh2_node_z.attrs[\"units\"] = \"km\"  # This is just a placeholder, UXarray does not support 3D meshes \n",
+    "# # Set correct units for the x and y coordinates\n",
+    "# vgrid.Mesh2_node_x.attrs[\"units\"] = \"km\"\n",
+    "# vgrid.Mesh2_node_y.attrs[\"units\"] = \"km\"\n",
+    "# vgrid.Mesh2_node_z.attrs[\"units\"] = \"km\"  # This is just a placeholder, UXarray does not support 3D meshes\n",
     "\n",
     "# Calculate the area of the triangle\n",
     "area_gaussian = vgrid.calculate_total_face_area(\n",
@@ -343,8 +340,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:29.022400Z",
-     "start_time": "2023-06-09T10:04:28.984401Z"
+     "start_time": "2023-09-16T09:12:07.151001Z",
+     "end_time": "2023-09-16T09:12:07.153518Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -375,8 +372,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:29.046402Z",
-     "start_time": "2023-06-09T10:04:29.001401Z"
+     "start_time": "2023-09-16T09:12:29.640457Z",
+     "end_time": "2023-09-16T09:12:29.645793Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -385,10 +382,8 @@
    },
    "outputs": [],
    "source": [
-    "verts_grid = ux.Grid(faces_verts_ndarray,\n",
-    "                     vertices=True,\n",
-    "                     islatlon=True,\n",
-    "                     concave=False)\n",
+    "verts_grid = ux.open_grid(faces_verts_ndarray, islatlon=False)\n",
+    "\n",
     "verts_grid"
    ]
   },
@@ -397,8 +392,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:29.078400Z",
-     "start_time": "2023-06-09T10:04:29.034401Z"
+     "start_time": "2023-09-16T09:13:10.699552Z",
+     "end_time": "2023-09-16T09:13:11.201163Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -407,7 +402,7 @@
    },
    "outputs": [],
    "source": [
-    "verts_grid.face_areas"
+    "verts_grid.compute_face_areas(latlon=False)"
    ]
   },
   {
@@ -424,8 +419,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:29.086401Z",
-     "start_time": "2023-06-09T10:04:29.048402Z"
+     "start_time": "2023-09-16T09:13:16.252149Z",
+     "end_time": "2023-09-16T09:13:16.258158Z"
     },
     "collapsed": false,
     "jupyter": {
@@ -449,8 +444,8 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:29.086401Z",
-     "start_time": "2023-06-09T10:04:29.064401Z"
+     "start_time": "2023-09-16T09:13:16.915445Z",
+     "end_time": "2023-09-16T09:13:16.921445Z"
     }
    },
    "outputs": [],
@@ -474,21 +469,14 @@
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-06-09T10:04:29.095401Z",
-     "start_time": "2023-06-09T10:04:29.081402Z"
+     "start_time": "2023-09-16T09:13:17.236352Z",
+     "end_time": "2023-09-16T09:13:18.687158Z"
     }
    },
    "outputs": [],
    "source": [
     "calculate_face_area(cart_x, cart_y, cart_z, coords_type=\"cartesian\")"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
index 3c76c423..c634036e 100644
--- a/docs/internal_api/index.rst
+++ b/docs/internal_api/index.rst
@@ -19,6 +19,7 @@ Class
 -----
 .. autosummary::
    :toctree: _autosummary
+
    UxDataset
 
 
@@ -83,25 +84,6 @@ Class
 
    Grid
 
-Methods
--------
-.. autosummary::
-   :toctree: _autosummary
-
-   Grid.__init_grid_var_names__
-   Grid.__from_ds__
-   Grid.__from_vert__
-   Grid.__init_grid_var_attrs__
-
-
-
-Attributes
-----------
-.. autosummary::
-   :toctree: _autosummary
-
-   Grid._Mesh2_node_z
-
 Operators
 ---------
 .. autosummary::
@@ -141,9 +123,12 @@ Coordinates
 .. autosummary::
    :toctree: _autosummary
 
+   grid.coordinates._get_lonlat_from_xyz
+   grid.coordinates._get_xyz_from_lonlat
    grid.coordinates._populate_cartesian_xyz_coord
    grid.coordinates._populate_lonlat_coord
 
+
 Lines
 -----
 .. autosummary::
@@ -164,6 +149,7 @@ UGRID
    io._ugrid._encode_ugrid
    io._ugrid._standardize_fill_value
    io._ugrid._is_ugrid
+   io._ugrid._validate_minimum_ugrid
 
 MPAS
 ----
@@ -205,9 +191,23 @@ Shapefile
 
    io._shapefile._read_shpfile
 
+Vertices
+--------
+.. autosummary::
+   :toctree: _autosummary
+
+   io._vertices._read_face_vertices
+
 Utils
 -----
 .. autosummary::
    :toctree: _autosummary
 
    io.utils._parse_grid_type
+
+Core Utils
+----------
+.. autosummary::
+   :toctree: _autosummary
+
+   core.utils._map_dims_to_ugrid
diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
index fa15f6da..00f9c2db 100644
--- a/docs/user_api/index.rst
+++ b/docs/user_api/index.rst
@@ -105,6 +105,8 @@ IO
    :toctree: _autosummary
 
    open_grid
+   Grid.from_dataset
+   Grid.from_face_vertices
    Grid.to_geodataframe
    Grid.to_polycollection
    Grid.to_shapely_polygons
@@ -128,6 +130,7 @@ Attributes
 .. autosummary::
    :toctree: _autosummary
 
+   Grid.grid_spec
    Grid.Mesh2
    Grid.parsed_attrs
    Grid.nMesh2_node
@@ -138,6 +141,9 @@ Attributes
    Grid.nNodes_per_face
    Grid.Mesh2_node_x
    Grid.Mesh2_node_y
+   Mesh2_node_cart_x
+   Mesh2_node_cart_y
+   Mesh2_node_cart_z
    Grid.Mesh2_face_x
    Grid.Mesh2_face_y
    Grid.Mesh2_face_nodes
@@ -167,6 +173,7 @@ Query Methods
    grid.BallTree.query_radius
 
 
+
 Helpers
 ===========
 
diff --git a/uxarray/core/api.py b/uxarray/core/api.py
index 1e1d904f..56aca8d9 100644
--- a/uxarray/core/api.py
+++ b/uxarray/core/api.py
@@ -4,20 +4,20 @@ import os
 import numpy as np
 import xarray as xr
 
-from pathlib import Path
+from pathlib import Path, PurePath
 from typing import Any, Dict, Optional, Union
 
 import uxarray.constants
 from uxarray.grid import Grid
 from uxarray.core.dataset import UxDataset
+from uxarray.core.utils import _map_dims_to_ugrid
+
+from warnings import warn
 
 
 def open_grid(grid_filename_or_obj: Union[str, Path, xr.DataArray, np.ndarray,
                                           list, tuple],
-              gridspec: Optional[str] = None,
-              vertices: Optional[list] = None,
-              islatlon: Optional[bool] = False,
-              isconcave: Optional[bool] = False,
+              latlon: Optional[bool] = False,
               use_dual: Optional[bool] = False,
               **kwargs: Dict[str, Any]) -> Grid:
     """Creates a ``uxarray.Grid`` object from a grid topology definition.
@@ -32,23 +32,11 @@ def open_grid(grid_filename_or_obj: Union[str, Path, xr.DataArray, np.ndarray,
         ``xr.DataArray``, ``np.ndarray``, ``list``, or ``tuple`` as a vertices
         object to define the grid.
 
-    islatlon : bool, optional
+    latlon : bool, optional
             Specify if the grid is lat/lon based
 
-    isconcave: bool, optional
-        Specify if this grid has concave elements (internal checks for this are possible)
-
-    gridspec: str, optional
-        Specifies gridspec
-
-    vertices: bool, optional
-        Whether to create grid from vertices
-
-    source_grid: str, optional
-        Path or URL to the source grid file. For diagnostic/reporting purposes only.
-
     use_dual: bool, optional
-            Specify whether to use the primal (use_dual=False) or dual (use_dual=True) mesh if the file type is mpas
+        Specify whether to use the primal (use_dual=False) or dual (use_dual=True) mesh if the file type is mpas
 
     **kwargs : Dict[str, Any]
         Additional arguments passed on to ``xarray.open_dataset``. Refer to the
@@ -71,38 +59,36 @@ def open_grid(grid_filename_or_obj: Union[str, Path, xr.DataArray, np.ndarray,
     >>> uxgrid = ux.open_grid("grid_filename.g")
     """
 
-    # Grid definition
-    if isinstance(grid_filename_or_obj,
-                  (list, tuple, np.ndarray, xr.DataArray)):
-        uxgrid = Grid(grid_filename_or_obj,
-                      gridspec=gridspec,
-                      vertices=vertices,
-                      islatlon=islatlon,
-                      isconcave=isconcave,
-                      source_grid=str(grid_filename_or_obj),
-                      use_dual=use_dual)
-    else:
+    if 'source_grid' in kwargs.keys():
+        warn('source_grid is no longer a supported kwarg',
+             DeprecationWarning,
+             stacklevel=2)
+
+    # construct Grid from dataset
+    if isinstance(grid_filename_or_obj, xr.Dataset):
+        uxgrid = Grid.from_dataset(grid_filename_or_obj, use_dual=use_dual)
+
+    # construct Grid from path
+    elif isinstance(grid_filename_or_obj, (str, Path, PurePath)):
         grid_ds = xr.open_dataset(grid_filename_or_obj,
                                   decode_times=False,
-                                  **kwargs)  # type: ignore
+                                  **kwargs)
+
+        uxgrid = Grid.from_dataset(grid_ds, use_dual=use_dual)
 
-        uxgrid = Grid(grid_ds,
-                      gridspec=gridspec,
-                      vertices=vertices,
-                      islatlon=islatlon,
-                      isconcave=isconcave,
-                      source_grid=str(grid_filename_or_obj),
-                      use_dual=use_dual)
+    elif isinstance(grid_filename_or_obj,
+                    (list, tuple, np.ndarray, xr.DataArray)):
+        uxgrid = Grid.from_face_vertices(grid_filename_or_obj, latlon=latlon)
+
+    else:
+        raise ValueError  # TODO: invalid input
 
     return uxgrid
 
 
 def open_dataset(grid_filename_or_obj: str,
                  filename_or_obj: str,
-                 gridspec: Optional[str] = None,
-                 vertices: Optional[list] = None,
-                 islatlon: Optional[bool] = False,
-                 isconcave: Optional[bool] = False,
+                 latlon: Optional[bool] = False,
                  use_dual: Optional[bool] = False,
                  grid_kwargs: Optional[Dict[str, Any]] = {},
                  **kwargs: Dict[str, Any]) -> UxDataset:
@@ -122,23 +108,11 @@ def open_dataset(grid_filename_or_obj: str,
         stores the actual data set. It is the same ``filename_or_obj`` in
         ``xarray.open_dataset``.
 
-    islatlon : bool, optional
+    latlon : bool, optional
             Specify if the grid is lat/lon based
 
-    isconcave: bool, optional
-        Specify if this grid has concave elements (internal checks for this are possible)
-
-    gridspec: str, optional
-        Specifies gridspec
-
-    vertices: bool, optional
-        Whether to create grid from vertices
-
-    source_grid: str, optional
-        Path or URL to the source grid file. For diagnostic/reporting purposes only.
-
     use_dual: bool, optional
-            Specify whether to use the primal (use_dual=False) or dual (use_dual=True) mesh if the file type is mpas
+        Specify whether to use the primal (use_dual=False) or dual (use_dual=True) mesh if the file type is mpas
 
     grid_kwargs : Dict[str, Any], optional
         Additional arguments passed on to ``xarray.open_dataset`` when opening up a Grid File. Refer to the
@@ -168,18 +142,24 @@ def open_dataset(grid_filename_or_obj: str,
     :param grid_kwargs:
     """
 
+    if 'source_grid' in kwargs.keys():
+        warn('source_grid is no longer a supported kwarg',
+             DeprecationWarning,
+             stacklevel=2)
+
     # Grid definition
     uxgrid = open_grid(grid_filename_or_obj,
-                       gridspec=gridspec,
-                       vertices=vertices,
-                       islatlon=islatlon,
-                       isconcave=isconcave,
+                       latlon=latlon,
                        use_dual=use_dual,
                        **grid_kwargs)
 
     # UxDataset
     ds = xr.open_dataset(filename_or_obj, decode_times=False,
                          **kwargs)  # type: ignore
+
+    # map each dimension to its UGRID equivalent
+    ds = _map_dims_to_ugrid(ds, uxgrid._source_dims_dict)
+
     uxds = UxDataset(ds, uxgrid=uxgrid, source_datasets=str(filename_or_obj))
 
     return uxds
@@ -187,10 +167,7 @@ def open_dataset(grid_filename_or_obj: str,
 
 def open_mfdataset(grid_filename_or_obj: str,
                    paths: Union[str, os.PathLike],
-                   gridspec: Optional[str] = None,
-                   vertices: Optional[list] = None,
-                   islatlon: Optional[bool] = False,
-                   isconcave: Optional[bool] = False,
+                   latlon: Optional[bool] = False,
                    use_dual: Optional[bool] = False,
                    grid_kwargs: Optional[Dict[str, Any]] = {},
                    **kwargs: Dict[str, Any]) -> UxDataset:
@@ -209,21 +186,9 @@ def open_mfdataset(grid_filename_or_obj: str,
         Either a string glob in the form "path/to/my/files/*.nc" or an explicit
         list of files to open. It is the same ``paths`` in ``xarray.open_mfdataset``.
 
-    islatlon : bool, optional
+    latlon : bool, optional
             Specify if the grid is lat/lon based
 
-    isconcave: bool, optional
-        Specify if this grid has concave elements (internal checks for this are possible)
-
-    gridspec: str, optional
-        Specifies gridspec
-
-    vertices: bool, optional
-        Whether to create grid from vertices
-
-    source_grid: str, optional
-        Path or URL to the source grid file. For diagnostic/reporting purposes only.
-
     use_dual: bool, optional
         Specify whether to use the primal (use_dual=False) or dual (use_dual=True) mesh if the file type is mpas
 
@@ -261,17 +226,23 @@ def open_mfdataset(grid_filename_or_obj: str,
     >>> ux_ds = ux.open_mfdataset("grid_filename.g", "grid_filename_vortex_*.nc")
     """
 
+    if 'source_grid' in kwargs.keys():
+        warn('source_grid is no longer a supported kwarg',
+             DeprecationWarning,
+             stacklevel=2)
+
     # Grid definition
     uxgrid = open_grid(grid_filename_or_obj,
-                       gridspec=gridspec,
-                       vertices=vertices,
-                       islatlon=islatlon,
-                       isconcave=isconcave,
+                       latlon=latlon,
                        use_dual=use_dual,
                        **grid_kwargs)
 
     # UxDataset
     ds = xr.open_mfdataset(paths, decode_times=False, **kwargs)  # type: ignore
+
+    # map each dimension to its UGRID equivalent
+    ds = _map_dims_to_ugrid(ds, uxgrid._source_dims_dict)
+
     uxds = UxDataset(ds, uxgrid=uxgrid, source_datasets=str(paths))
 
     return uxds
diff --git a/uxarray/core/utils.py b/uxarray/core/utils.py
new file mode 100644
index 00000000..ee5b8e20
--- /dev/null
+++ b/uxarray/core/utils.py
@@ -0,0 +1,18 @@
+def _map_dims_to_ugrid(ds, _source_dims_dict):
+    """Given a dataset containing variables residing on an unstructured grid, remaps the original dimension name
+    to match the UGRID conventions (i.e. "nCell" : "nMesh2_Face")"""
+
+    keys_to_drop = []
+    # obtain all dimensions not present in the original dataset
+    for key in _source_dims_dict.keys():
+        if key not in ds.dims:
+            keys_to_drop.append(key)
+
+    # drop dimensions not present in the original dataset
+    for key in keys_to_drop:
+        _source_dims_dict.pop(key)
+
+    # rename dimensions to follow the UGRID conventions
+    ds = ds.rename_dims(_source_dims_dict)
+
+    return ds
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
index 99fd5852..f0e823d9 100644
--- a/uxarray/grid/coordinates.py
+++ b/uxarray/grid/coordinates.py
@@ -107,6 +107,17 @@ def normalize_in_place(node):
     return list(np.array(node) / np.linalg.norm(np.array(node), ord=2))
 
 
+def _get_xyz_from_lonlat(node_lon, node_lat):
+
+    # check for units and create Mesh2_node_cart_x/y/z set to grid._ds
+    nodes_lon_rad = np.deg2rad(node_lon)
+    nodes_lat_rad = np.deg2rad(node_lat)
+    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))))
+
+    return nodes_cart[:, 0], nodes_cart[:, 1], nodes_cart[:, 2]
+
+
 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
@@ -136,28 +147,26 @@ def _populate_cartesian_xyz_coord(grid):
     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))))
+    # get Cartesian (x, y, z) coordinates from lon/lat
+    x, y, z = _get_xyz_from_lonlat(grid.Mesh2_node_x.values,
+                                   grid.Mesh2_node_y.values)
 
     grid._ds["Mesh2_node_cart_x"] = xr.DataArray(
-        data=nodes_cart[:, 0],
+        data=x,
         dims=["nMesh2_node"],
         attrs={
             "standard_name": "cartesian x",
             "units": "m",
         })
     grid._ds["Mesh2_node_cart_y"] = xr.DataArray(
-        data=nodes_cart[:, 1],
+        data=y,
         dims=["nMesh2_node"],
         attrs={
             "standard_name": "cartesian y",
             "units": "m",
         })
     grid._ds["Mesh2_node_cart_z"] = xr.DataArray(
-        data=nodes_cart[:, 2],
+        data=z,
         dims=["nMesh2_node"],
         attrs={
             "standard_name": "cartesian z",
@@ -165,6 +174,14 @@ def _populate_cartesian_xyz_coord(grid):
         })
 
 
+def _get_lonlat_from_xyz(x, y, z):
+    nodes_cart = np.stack((x, y, z), axis=1).tolist()
+    nodes_rad = list(map(node_xyz_to_lonlat_rad, nodes_cart))
+    nodes_degree = np.rad2deg(nodes_rad)
+
+    return nodes_degree[:, 0], nodes_degree[:, 1]
+
+
 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
@@ -196,34 +213,14 @@ def _populate_lonlat_coord(grid):
      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)
+    # get lon/lat coordinates from Cartesian (x, y, z)
+    lon, lat = _get_lonlat_from_xyz(grid.Mesh2_node_cart_x.values,
+                                    grid.Mesh2_node_cart_y.values,
+                                    grid.Mesh2_node_cart_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)
+    # populate dataset
     grid._ds["Mesh2_node_x"] = xr.DataArray(
-        data=nodes_degree[:, 0],
+        data=lon,
         dims=["nMesh2_node"],
         attrs={
             "standard_name": "longitude",
@@ -231,7 +228,7 @@ def _populate_lonlat_coord(grid):
             "units": "degrees_east",
         })
     grid._ds["Mesh2_node_y"] = xr.DataArray(
-        data=nodes_degree[:, 1],
+        data=lat,
         dims=["nMesh2_node"],
         attrs={
             "standard_name": "latitude",
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 27f45387..82d5c889 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -7,87 +7,108 @@ from typing import Any, Dict, Optional, Union
 # reader and writer imports
 from uxarray.io._exodus import _read_exodus, _encode_exodus
 from uxarray.io._mpas import _read_mpas
-from uxarray.io._ugrid import _read_ugrid, _encode_ugrid
+from uxarray.io._ugrid import _read_ugrid, _encode_ugrid, _validate_minimum_ugrid
 from uxarray.io._shapefile import _read_shpfile
 from uxarray.io._scrip import _read_scrip, _encode_scrip
+from uxarray.io._vertices import _read_face_vertices
 
 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,
-                                       _build_node_faces_connectivity,
-                                       _face_nodes_to_sparse_matrix)
+from uxarray.grid.connectivity import (
+    _build_edge_node_connectivity,
+    _build_face_edges_connectivity,
+    _build_nNodes_per_face,
+    _build_node_faces_connectivity,
+)
 
-from uxarray.grid.coordinates import (normalize_in_place,
-                                      _populate_lonlat_coord,
+from uxarray.grid.coordinates import (_populate_lonlat_coord,
                                       _populate_cartesian_xyz_coord)
 
-from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
-
-from uxarray.grid.geometry import (_build_polygon_shells,
-                                   _build_corrected_polygon_shells,
-                                   _build_antimeridian_face_indices,
+from uxarray.grid.geometry import (_build_antimeridian_face_indices,
                                    _grid_to_polygon_geodataframe,
                                    _grid_to_matplotlib_polycollection,
                                    _grid_to_polygons)
 
 from uxarray.grid.neighbors import BallTree
 
+from warnings import warn
+
 
 class Grid:
-    """Unstructured grid topology definition.
+    """Represents a two-dimensional unstructured grid encoded following the
+    UGRID conventions and provides grid-specific functionality.
+
+    Can be used standalone to work with unstructured grids, or can be paired with either a ``ux.UxDataArray`` or
+    ``ux.UxDataset`` and accessed through the ``.uxgrid`` attribute.
+
+    For constructing a grid from non-UGRID datasets or other types of supported data, see our ``ux.open_grid`` method or
+    specific class methods (``Grid.from_dataset``, ``Grid.from_face_verticies``, etc.)
 
-    Can be used standalone to explore an unstructured grid topology, or
-    can be seen as the property of ``uxarray.UxDataset`` and ``uxarray.DataArray``
-    to make them unstructured grid-aware data sets and arrays.
 
     Parameters
     ----------
-    input_obj : xarray.Dataset, ndarray, list, tuple, required
-        Input ``xarray.Dataset`` or vertex coordinates that form faces.
-
-    Other Parameters
-    ----------------
-    gridspec: bool, optional
-        Specifies gridspec
-    islatlon : bool, optional
-        Specify if the grid is lat/lon based
-    isconcave: bool, optional
-        Specify if this grid has concave elements (internal checks for this are possible)
-    use_dual: bool, optional
-        Specify whether to use the primal (use_dual=False) or dual (use_dual=True) mesh if the file type is mpas
-
-    Raises
-    ------
-        RuntimeError
-            If specified file not found or recognized
+    grid_ds : xr.Dataset
+        ``xarray.Dataset`` encoded in the UGRID conventions
+
+    source_grid_spec : str, default="UGRID"
+        Original unstructured grid format (i.e. UGRID, MPAS, etc.)
+
+    source_dims_dict : dict, default={}
+        Mapping of dimensions from the source dataset to their UGRID equivalent (i.e. {nCell : nMesh2_face})
 
     Examples
     ----------
 
     >>> import uxarray as ux
+    >>> grid_path = "/path/to/grid.nc"
+    >>> data_path = "/path/to/data.nc"
 
     1. Open a grid file with `uxarray.open_grid()`:
 
-    >>> uxgrid = ux.open_grid("filename.g")
+    >>> uxgrid = ux.open_grid(grid_path)
 
     2. Open an unstructured grid dataset file with
-    `uxarray.open_dataset()`, then access `Grid` info:
+    `uxarray.open_dataset()`, then access the ``Grid``.:
 
-    >>> uxds = ux.open_dataset("filename.g")
+    >>> uxds = ux.open_dataset(grid_path, data_path)
+    >>> uxds.uxgrid
     """
 
-    def __init__(self, input_obj, **kwargs):
-        # initialize internal variable names
-        self.__init_grid_var_names__()
-
-        # initialize face_area variable
-        self._face_areas = None
+    def __init__(self,
+                 grid_ds: xr.Dataset,
+                 source_grid_spec: Optional[str] = None,
+                 source_dims_dict: Optional[dict] = {}):
+
+        # check if inputted dataset is a minimum representable 2D UGRID unstructured grid
+        if not _validate_minimum_ugrid(grid_ds):
+            raise ValueError(
+                "Direct use of Grid constructor requires grid_ds to follow the internal unstructured grid definition, "
+                "including variable and dimension names. This grid_ds does not satisfy those requirements. If you are "
+                "not sure about how to do that, using ux.open_grid() or ux.from_dataset() is suggested."
+            )  # TODO: elaborate once we have a formal definition
+
+        # grid spec not provided, check if grid_ds is a minimum representable UGRID dataset
+        if source_grid_spec is None:
+            warn(
+                "Attempting to construct a Grid without passing in source_grid_spec. Direct use of Grid constructor"
+                "is only advised if grid_ds is following the internal unstructured grid definition, including"
+                "variable and dimension names. Using ux.open_grid() or ux.from_dataset() is suggested.",
+                Warning)
+            # TODO: more checks for validate grid (lat/lon coords, etc)
+
+        # mapping of ugrid dimensions and variables to source dataset's conventions
+        self._source_dims_dict = source_dims_dict
+
+        # source grid specification (i.e. UGRID, MPAS, SCRIP, etc.)
+        self.source_grid_spec = source_grid_spec
+
+        # internal xarray dataset for storing grid variables
+        self._ds = grid_ds
 
         # initialize attributes
         self._antimeridian_face_indices = None
+        self._face_areas = None
 
         # initialize cached data structures (visualization)
         self._gdf = None
@@ -96,211 +117,77 @@ class Grid:
         # initialize cached data structures (nearest neighbor operations)
         self._ball_tree = None
 
-        # unpack kwargs with default values set to None
-        kwargs_list = [
-            'gridspec', 'vertices', 'islatlon', 'isconcave', 'source_grid',
-            'use_dual'
-        ]
-        for key in kwargs_list:
-            setattr(self, key, kwargs.get(key, None))
-
-        # check if initializing from verts:
-        if isinstance(input_obj, (list, tuple, np.ndarray)):
-            input_obj = np.asarray(input_obj)
-            self.mesh_type = "From vertices"
-            # grid with multiple faces
-            if input_obj.ndim == 3:
-                self.__from_vert__(input_obj)
-                self.source_grid = "From vertices"
-            # grid with a single face
-            elif input_obj.ndim == 2:
-                input_obj = np.array([input_obj])
-                self.__from_vert__(input_obj)
-                self.source_grid = "From vertices"
-            else:
-                raise RuntimeError(
-                    f"Invalid Input Dimension: {input_obj.ndim}. Expected dimension should be "
-                    f"3: [nMesh2_face, nMesh2_node, Two/Three] or 2 when only "
-                    f"one face is passed in.")
-
-        # 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.__from_ds__(dataset=input_obj)
-        else:
-            raise RuntimeError("Dataset is not a valid input type.")
+    @classmethod
+    def from_dataset(cls,
+                     dataset: xr.Dataset,
+                     use_dual: Optional[bool] = False):
+        """Constructs a ``Grid`` object from an ``xarray.Dataset``.
 
-        # {"Standardized Name" : "Original Name"}
-        self._inverse_grid_var_names = {
-            v: k for k, v in self.grid_var_names.items()
-        }
+        Parameters
+        ----------
+        dataset : xr.Dataset
+            ``xarray.Dataset`` containing unstructured grid coordinates and connectivity variables
+        use_dual : bool, default=False
+            When reading in MPAS formatted datasets, indicates whether to use the Dual Mesh
+        """
+        if not isinstance(dataset, xr.Dataset):
+            raise ValueError("Input must be an xarray.Dataset")
+
+        # determine grid/mesh specification
+        source_grid_spec = _parse_grid_type(dataset)
+
+        if source_grid_spec == "Exodus":
+            grid_ds, source_dims_dict = _read_exodus(dataset)
+        elif source_grid_spec == "Scrip":
+            grid_ds, source_dims_dict = _read_scrip(dataset)
+        elif source_grid_spec == "UGRID":
+            grid_ds, source_dims_dict = _read_ugrid(dataset)
+        elif source_grid_spec == "MPAS":
+            grid_ds, source_dims_dict = _read_mpas(dataset, use_dual=use_dual)
+        elif source_grid_spec == "Shapefile":
+            raise ValueError("Shapefiles not yet supported")
+        else:
+            raise ValueError("Unsupported Grid Format")
 
-    def __init_grid_var_names__(self):
-        """Populates a dictionary for storing uxarray's internal representation
-        of xarray object.
+        return cls(grid_ds, source_grid_spec, source_dims_dict)
 
-        Note ugrid conventions are flexible with names of variables, see:
-        http://ugrid-conventions.github.io/ugrid-conventions/
-        """
-        self.grid_var_names = {
-            "Mesh2": "Mesh2",
-            "Mesh2_node_x": "Mesh2_node_x",
-            "Mesh2_node_y": "Mesh2_node_y",
-            "Mesh2_node_z": "Mesh2_node_z",
-            "Mesh2_face_nodes": "Mesh2_face_nodes",
-            # initialize dims
-            "nMesh2_node": "nMesh2_node",
-            "nMesh2_face": "nMesh2_face",
-            "nMaxMesh2_face_nodes": "nMaxMesh2_face_nodes"
-        }
-
-    def __from_vert__(self, dataset):
-        """Create a grid with faces constructed from vertices specified by the
-        given argument.
+    @classmethod
+    def from_face_vertices(cls,
+                           face_vertices: Union[list, tuple, np.ndarray],
+                           latlon: Optional[bool] = True):
+        """Constructs a ``Grid`` object from user-defined face vertices.
 
         Parameters
         ----------
-        dataset : ndarray, list, tuple, required
-            Input vertex coordinates that form our face(s)
+        face_vertices : list, tuple, np.ndarray
+            array-like input containing the face vertices to construct the grid from
+        latlon : bool, default=True
+            Indicates whether the inputted vertices are in lat/lon, with units in degrees
         """
-        self._ds = xr.Dataset()
-        self._ds["Mesh2"] = xr.DataArray(
-            attrs={
-                "cf_role": "mesh_topology",
-                "long_name": "Topology data of unstructured mesh",
-                "topology_dimension": -1,
-                "node_coordinates": "Mesh2_node_x Mesh2_node_y Mesh2_node_z",
-                "node_dimension": "nMesh2_node",
-                "face_node_connectivity": "Mesh2_face_nodes",
-                "face_dimension": "nMesh2_face"
-            })
-
-        self._ds.Mesh2.attrs['topology_dimension'] = dataset.ndim
-
-        if self.islatlon is not None and self.islatlon is False:
-            x_units = 'm'
-            y_units = 'm'
-            z_units = 'm'
-        else:
-            x_units = "degrees_east"
-            y_units = "degrees_north"
-            z_units = "elevation"
+        if not isinstance(face_vertices, (list, tuple, np.ndarray)):
+            raise ValueError(
+                "Input must be either a list, tuple, or np.ndarray")
 
-        x_coord = dataset[:, :, 0].flatten()
-        y_coord = dataset[:, :, 1].flatten()
+        face_vertices = np.asarray(face_vertices)
 
-        if dataset[0][0].size > 2:
-            z_coord = dataset[:, :, 2].flatten()
-        else:
-            z_coord = x_coord * 0.0
-
-        # Identify unique vertices and their indices
-        unique_verts, indices = np.unique(dataset.reshape(
-            -1, dataset.shape[-1]),
-                                          axis=0,
-                                          return_inverse=True)
-
-        # Nodes index that contain a fill value
-        fill_value_mask = np.logical_or(unique_verts[:, 0] == INT_FILL_VALUE,
-                                        unique_verts[:, 1] == INT_FILL_VALUE)
-        if dataset[0][0].size > 2:
-            fill_value_mask = np.logical_or(
-                unique_verts[:, 0] == INT_FILL_VALUE,
-                unique_verts[:, 1] == INT_FILL_VALUE,
-                unique_verts[:, 2] == INT_FILL_VALUE)
-
-        # Get the indices of all the False values in fill_value_mask
-        false_indices = np.where(fill_value_mask == True)[0]
-
-        # Check if any False values were found
-        indices = indices.astype(INT_DTYPE)
-        if false_indices.size > 0:
-
-            # Remove the rows corresponding to False values in unique_verts
-            unique_verts = np.delete(unique_verts, false_indices, axis=0)
-
-            # Update indices accordingly
-            for i, idx in enumerate(false_indices):
-                indices[indices == idx] = INT_FILL_VALUE
-                indices[(indices > idx) & (indices != INT_FILL_VALUE)] -= 1
-
-        # Create coordinate DataArrays
-        self._ds["Mesh2_node_x"] = xr.DataArray(data=unique_verts[:, 0],
-                                                dims=["nMesh2_node"],
-                                                attrs={"units": x_units})
-        self._ds["Mesh2_node_y"] = xr.DataArray(data=unique_verts[:, 1],
-                                                dims=["nMesh2_node"],
-                                                attrs={"units": y_units})
-        if dataset.shape[-1] > 2:
-            self._ds["Mesh2_node_z"] = xr.DataArray(data=unique_verts[:, 2],
-                                                    dims=["nMesh2_node"],
-                                                    attrs={"units": z_units})
-        else:
-            self._ds["Mesh2_node_z"] = xr.DataArray(data=unique_verts[:, 1] *
-                                                    0.0,
-                                                    dims=["nMesh2_node"],
-                                                    attrs={"units": z_units})
-
-        # Create connectivity array using indices of unique vertices
-        connectivity = indices.reshape(dataset.shape[:-1])
-        self._ds["Mesh2_face_nodes"] = xr.DataArray(
-            data=xr.DataArray(connectivity).astype(INT_DTYPE),
-            dims=["nMesh2_face", "nMaxMesh2_face_nodes"],
-            attrs={
-                "cf_role": "face_node_connectivity",
-                "_FillValue": INT_FILL_VALUE,
-                "start_index": 0
-            })
-
-    # load mesh from a file
-    def __from_ds__(self, dataset):
-        """Loads a mesh dataset."""
-        # call reader as per mesh_type
-        if self.mesh_type == "exo":
-            self._ds = _read_exodus(dataset, self.grid_var_names)
-
-            # Assume Exodus was read as cartesian grid and that coordinates are not set by reader, call the latlon setter
-            _populate_lonlat_coord(self)
+        if face_vertices.ndim == 3:
+            grid_ds = _read_face_vertices(face_vertices, latlon)
+
+        elif face_vertices.ndim == 2:
+            grid_ds = _read_face_vertices(np.array([face_vertices]), latlon)
 
-            # set coordinates
-            # there is leftover cartesian z-coordinate from Exodus mesh
-            # set them to zero
-            self._ds["Mesh2_node_z"] = xr.DataArray(
-                data=np.zeros(self._ds["Mesh2_node_x"].shape),
-                dims=["nMesh2_node"],
-                attrs={
-                    "standard_name": "elevation",
-                    "long_name": "elevation",
-                    "units": "m",
-                })
-
-            ds = self._ds.set_coords(
-                ["Mesh2_node_x", "Mesh2_node_y", "Mesh2_node_z"])
-
-        elif self.mesh_type == "scrip":
-            self._ds = _read_scrip(dataset)
-        elif self.mesh_type == "ugrid":
-            self._ds, self.grid_var_names = _read_ugrid(dataset,
-                                                        self.grid_var_names)
-        elif self.mesh_type == "shp":
-            self._ds = _read_shpfile(dataset)
-        elif self.mesh_type == "mpas":
-            # select whether to use the dual mesh
-            if self.use_dual is not None:
-                self._ds = _read_mpas(dataset, self.use_dual)
-            else:
-                self._ds = _read_mpas(dataset)
         else:
-            raise RuntimeError("unknown mesh type")
+            raise RuntimeError(
+                f"Invalid Input Dimension: {face_vertices.ndim}. Expected dimension should be "
+                f"3: [nMesh2_face, nMesh2_node, Two/Three] or 2 when only "
+                f"one face is passed in.")
 
-        dataset.close()
+        return cls(grid_ds, source_grid_spec="Face Vertices")
 
     def __repr__(self):
         """Constructs a string representation of the contents of a ``Grid``."""
-
         prefix = "<uxarray.Grid>\n"
-        original_grid_str = f"Original Grid Type: {self.mesh_type}\n"
+        original_grid_str = f"Original Grid Type: {self.source_grid_spec}\n"
         dims_heading = "Grid Dimensions:\n"
         dims_str = ""
         # if self.grid_var_names["Mesh2_node_x"] in self._ds:
@@ -310,8 +197,9 @@ class Grid:
         #     dims_str += f"  * nMesh2_face: {self.nMesh2_face}\n"
 
         for key, value in zip(self._ds.dims.keys(), self._ds.dims.values()):
-            if key in self._inverse_grid_var_names:
-                dims_str += f"  * {self._inverse_grid_var_names[key]}: {value}\n"
+            dims_str += f"  * {key}: {value}\n"
+            # if key in self._inverse_grid_var_names:
+            #     dims_str += f"  * {self._inverse_grid_var_names[key]}: {value}\n"
 
         if "nMesh2_edge" in self._ds.dims:
             dims_str += f"  * nMesh2_edge: {self.nMesh2_edge}\n"
@@ -321,7 +209,7 @@ class Grid:
 
         coord_heading = "Grid Coordinate Variables:\n"
         coords_str = ""
-        if self.grid_var_names["Mesh2_node_x"] in self._ds:
+        if "Mesh2_node_x" in self._ds:
             coords_str += f"  * Mesh2_node_x: {self.Mesh2_node_x.shape}\n"
             coords_str += f"  * Mesh2_node_y: {self.Mesh2_node_y.shape}\n"
         if "Mesh2_node_cart_x" in self._ds:
@@ -334,7 +222,7 @@ class Grid:
 
         connectivity_heading = "Grid Connectivity Variables:\n"
         connectivity_str = ""
-        if self.grid_var_names["Mesh2_face_nodes"] in self._ds:
+        if "Mesh2_face_nodes" in self._ds:
             connectivity_str += f"  * Mesh2_face_nodes: {self.Mesh2_face_nodes.shape}\n"
         if "Mesh2_edge_nodes" in self._ds:
             connectivity_str += f"  * Mesh2_edge_nodes: {self.Mesh2_edge_nodes.shape}\n"
@@ -345,7 +233,7 @@ class Grid:
         return prefix + original_grid_str + dims_heading + dims_str + coord_heading + coords_str + \
             connectivity_heading + connectivity_str
 
-    def __eq__(self, other):
+    def __eq__(self, other) -> bool:
         """Two grids are equal if they have matching grid topology variables,
         coordinates, and dims all of which are equal.
 
@@ -358,21 +246,23 @@ class Grid:
         -------
         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:
+
+        if not isinstance(other, Grid):
+            return False
+
+        if self.source_grid_spec != other.source_grid_spec:
+            return False
+
+        if not (self.Mesh2_node_x.equals(other.Mesh2_node_x) or
+                self.Mesh2_node_y.equals(other.Mesh2_node_y)):
+            return False
+
+        if not self.Mesh2_face_nodes.equals(other.Mesh2_face_nodes):
             return False
 
         return True
 
-    def __ne__(self, other):
+    def __ne__(self, other) -> bool:
         """Two grids are not equal if they have differing grid topology
         variables, coordinates, or dims.
 
@@ -388,30 +278,30 @@ class Grid:
         return not self.__eq__(other)
 
     @property
-    def parsed_attrs(self):
+    def parsed_attrs(self) -> dict:
         """Dictionary of parsed attributes from the source grid."""
         return self._ds.attrs
 
     @property
-    def Mesh2(self):
+    def Mesh2(self) -> xr.DataArray:
         """UGRID Attribute ``Mesh2``, which indicates the topology data of a 2D
         unstructured mesh."""
-        return self._ds[self.grid_var_names["Mesh2"]]
+        return self._ds["Mesh2"]
 
     @property
-    def nMesh2_node(self):
+    def nMesh2_node(self) -> int:
         """UGRID Dimension ``nMesh2_node``, which represents the total number
         of nodes."""
-        return self._ds[self.grid_var_names["Mesh2_node_x"]].shape[0]
+        return self._ds.dims["nMesh2_node"]
 
     @property
-    def nMesh2_face(self):
+    def nMesh2_face(self) -> int:
         """UGRID Dimension ``nMesh2_face``, which represents the total number
         of faces."""
-        return self._ds[self.grid_var_names["Mesh2_face_nodes"]].shape[0]
+        return self._ds["Mesh2_face_nodes"].shape[0]
 
     @property
-    def nMesh2_edge(self):
+    def nMesh2_edge(self) -> int:
         """UGRID Dimension ``nMesh2_edge``, which represents the total number
         of edges."""
 
@@ -421,13 +311,13 @@ class Grid:
         return self._ds['Mesh2_edge_nodes'].shape[0]
 
     @property
-    def nMaxMesh2_face_nodes(self):
+    def nMaxMesh2_face_nodes(self) -> int:
         """UGRID Dimension ``nMaxMesh2_face_nodes``, which represents the
         maximum number of nodes that a face may contain."""
         return self.Mesh2_face_nodes.shape[1]
 
     @property
-    def nMaxMesh2_face_edges(self):
+    def nMaxMesh2_face_edges(self) -> xr.DataArray:
         """Dimension ``nMaxMesh2_face_edges``, which represents the maximum
         number of edges per face.
 
@@ -440,7 +330,7 @@ class Grid:
         return self._ds["Mesh2_face_edges"].shape[1]
 
     @property
-    def nNodes_per_face(self):
+    def nNodes_per_face(self) -> xr.DataArray:
         """Dimension Variable ``nNodes_per_face``, which contains the number of
         non-fill-value nodes per face.
 
@@ -451,16 +341,18 @@ class Grid:
         return self._ds["nNodes_per_face"]
 
     @property
-    def Mesh2_node_x(self):
+    def Mesh2_node_x(self) -> xr.DataArray:
         """UGRID Coordinate Variable ``Mesh2_node_x``, which contains the
         longitude of each node in degrees.
 
         Dimensions (``nMesh2_node``)
         """
-        return self._ds[self.grid_var_names["Mesh2_node_x"]]
+        if "Mesh2_node_x" not in self._ds:
+            _populate_lonlat_coord(self)
+        return self._ds["Mesh2_node_x"]
 
     @property
-    def Mesh2_node_cart_x(self):
+    def Mesh2_node_cart_x(self) -> xr.DataArray:
         """Coordinate Variable ``Mesh2_node_cart_x``, which contains the x
         location in meters.
 
@@ -468,10 +360,11 @@ class Grid:
         """
         if "Mesh2_node_cart_x" not in self._ds:
             _populate_cartesian_xyz_coord(self)
+
         return self._ds['Mesh2_node_cart_x']
 
     @property
-    def Mesh2_face_x(self):
+    def Mesh2_face_x(self) -> xr.DataArray:
         """UGRID Coordinate Variable ``Mesh2_face_x``, which contains the
         longitude of each face center.
 
@@ -483,16 +376,19 @@ class Grid:
             return None
 
     @property
-    def Mesh2_node_y(self):
+    def Mesh2_node_y(self) -> xr.DataArray:
         """UGRID Coordinate Variable ``Mesh2_node_y``, which contains the
         latitude of each node.
 
         Dimensions (``nMesh2_node``)
         """
-        return self._ds[self.grid_var_names["Mesh2_node_y"]]
+        if "Mesh2_node_y" not in self._ds:
+            _populate_lonlat_coord(self)
+
+        return self._ds["Mesh2_node_y"]
 
     @property
-    def Mesh2_node_cart_y(self):
+    def Mesh2_node_cart_y(self) -> xr.DataArray:
         """Coordinate Variable ``Mesh2_node_cart_y``, which contains the y
         location in meters.
 
@@ -503,7 +399,7 @@ class Grid:
         return self._ds['Mesh2_node_cart_y']
 
     @property
-    def Mesh2_face_y(self):
+    def Mesh2_face_y(self) -> xr.DataArray:
         """UGRID Coordinate Variable ``Mesh2_face_y``, which contains the
         latitude of each face center.
 
@@ -515,25 +411,7 @@ class Grid:
             return None
 
     @property
-    def _Mesh2_node_z(self):
-        """Coordinate Variable ``_Mesh2_node_z``, which contains the level of
-        each node. It is only a placeholder for now as a protected attribute.
-        UXarray does not support this yet and only handles the 2D flexibile
-        meshes.
-
-        If we introduce handling of 3D meshes in the future, it might be only
-        levels, i.e. the same level(s) for all nodes, instead of separate
-        level for each node that ``_Mesh2_node_z`` suggests.
-
-        Dimensions (``nMesh2_node``)
-        """
-        if self.grid_var_names["Mesh2_node_z"] in self._ds:
-            return self._ds[self.grid_var_names["Mesh2_node_z"]]
-        else:
-            return None
-
-    @property
-    def Mesh2_node_cart_z(self):
+    def Mesh2_node_cart_z(self) -> xr.DataArray:
         """Coordinate Variable ``Mesh2_node_cart_z``, which contains the z
         location in meters.
 
@@ -544,7 +422,7 @@ class Grid:
         return self._ds['Mesh2_node_cart_z']
 
     @property
-    def Mesh2_face_nodes(self):
+    def Mesh2_face_nodes(self) -> xr.DataArray:
         """UGRID Connectivity Variable ``Mesh2_face_nodes``, which maps each
         face to its corner nodes.
 
@@ -557,10 +435,10 @@ class Grid:
         Nodes are in counter-clockwise order.
         """
 
-        return self._ds[self.grid_var_names["Mesh2_face_nodes"]]
+        return self._ds["Mesh2_face_nodes"]
 
     @property
-    def Mesh2_edge_nodes(self):
+    def Mesh2_edge_nodes(self) -> xr.DataArray:
         """UGRID Connectivity Variable ``Mesh2_edge_nodes``, which maps every
         edge to the two nodes that it connects.
 
@@ -575,7 +453,7 @@ class Grid:
         return self._ds['Mesh2_edge_nodes']
 
     @property
-    def Mesh2_face_edges(self):
+    def Mesh2_face_edges(self) -> xr.DataArray:
         """UGRID Connectivity Variable ``Mesh2_face_edges``, which maps every
         face to its edges.
 
@@ -587,17 +465,8 @@ class Grid:
 
         return self._ds["Mesh2_face_edges"]
 
-    # other properties
     @property
-    def antimeridian_face_indices(self):
-        """Index of each face that crosses the antimeridian."""
-        if self._antimeridian_face_indices is None:
-            self._antimeridian_face_indices = _build_antimeridian_face_indices(
-                self)
-        return self._antimeridian_face_indices
-
-    @property
-    def Mesh2_node_faces(self):
+    def Mesh2_node_faces(self) -> xr.DataArray:
         """UGRID Connectivity Variable ``Mesh2_node_faces``, which maps every
         node to its faces.
 
@@ -609,8 +478,17 @@ class Grid:
 
         return self._ds["Mesh2_node_faces"]
 
+    # other properties
+    @property
+    def antimeridian_face_indices(self) -> np.ndarray:
+        """Index of each face that crosses the antimeridian."""
+        if self._antimeridian_face_indices is None:
+            self._antimeridian_face_indices = _build_antimeridian_face_indices(
+                self)
+        return self._antimeridian_face_indices
+
     @property
-    def face_areas(self):
+    def face_areas(self) -> np.ndarray:
         """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:
@@ -646,15 +524,12 @@ class Grid:
 
     def copy(self):
         """Returns a deep copy of this grid."""
-        return Grid(xr.Dataset(self._ds),
-                    gridspec=self.gridspec,
-                    vertices=self.vertices,
-                    islatlon=self.islatlon,
-                    isconcave=self.isconcave,
-                    source_grid=self.source_grid,
-                    use_dual=self.use_dual)
-
-    def encode_as(self, grid_type):
+
+        return Grid(self._ds,
+                    source_grid_spec=self.source_grid_spec,
+                    source_dims_dict=self._source_dims_dict)
+
+    def encode_as(self, grid_type: str) -> xr.Dataset:
         """Encodes the grid as a new `xarray.Dataset` per grid format supplied
         in the `grid_type` argument.
 
@@ -675,20 +550,13 @@ class Grid:
             If provided grid type or file type is unsupported.
         """
 
-        if grid_type == "ugrid":
+        if grid_type == "UGRID":
             out_ds = _encode_ugrid(self._ds)
 
-        elif grid_type == "exodus":
+        elif grid_type == "Exodus":
+            out_ds = _encode_exodus(self._ds)
 
-            # NOTE: We assume that output exodus mesh will be cartesian and coordinate units will be 'm'
-            # If the units are rad or degree, the we must convert to m. Assume unit sphere.
-            if "Mesh2_node_cart_x" not in self._ds.keys():
-                _populate_cartesian_xyz_coord(self)
-
-            # encode to exodus assumes that ds has Mesh2_node_cart_x, Mesh2_node_cart_y, Mesh2_node_cart_z
-            out_ds = _encode_exodus(self._ds, self.grid_var_names)
-
-        elif grid_type == "scrip":
+        elif grid_type == "SCRIP":
             out_ds = _encode_scrip(self.Mesh2_face_nodes, self.Mesh2_node_x,
                                    self.Mesh2_node_y, self.face_areas)
         else:
@@ -696,7 +564,9 @@ class Grid:
 
         return out_ds
 
-    def calculate_total_face_area(self, quadrature_rule="triangular", order=4):
+    def calculate_total_face_area(self,
+                                  quadrature_rule: Optional[str] = "triangular",
+                                  order: Optional[int] = 4) -> float:
         """Function to calculate the total surface area of all the faces in a
         mesh.
 
@@ -717,7 +587,10 @@ class Grid:
 
         return np.sum(face_areas)
 
-    def compute_face_areas(self, quadrature_rule="triangular", order=4):
+    def compute_face_areas(self,
+                           quadrature_rule: Optional[str] = "triangular",
+                           order: Optional[int] = 4,
+                           latlon: Optional[bool] = True):
         """Face areas calculation function for grid class, calculates area of
         all faces in the grid.
 
@@ -747,24 +620,22 @@ class Grid:
         # if self._face_areas is None: # this allows for using the cached result,
         # but is not the expected behavior behavior as we are in need to recompute if this function is called with different quadrature_rule or order
 
-        # area of a face call needs the units for coordinate conversion if spherical grid is used
-        coords_type = "spherical"
-        if not "degree" in self.Mesh2_node_x.units:
+        if latlon:
+            x = self.Mesh2_node_x.data
+            y = self.Mesh2_node_y.data
+            z = np.zeros((self.nMesh2_node))
+            coords_type = "spherical"  # TODO: should really be called latlon?
+        else:
+            x = self.Mesh2_node_cart_x.data
+            y = self.Mesh2_node_cart_y.data
+            z = self.Mesh2_node_cart_z.data
             coords_type = "cartesian"
 
-        face_nodes = self.Mesh2_face_nodes.data
-        nNodes_per_face = self.nNodes_per_face.data
+        # TODO: we dont really need this, but keep for now
         dim = self.Mesh2.attrs['topology_dimension']
 
-        # initialize z
-        z = np.zeros((self.nMesh2_node))
-
-        # call func to cal face area of all nodes
-        x = self.Mesh2_node_x.data
-        y = self.Mesh2_node_y.data
-        # check if z dimension
-        if self.Mesh2.topology_dimension > 2:
-            z = self._Mesh2_node_z.data
+        nNodes_per_face = self.nNodes_per_face.data
+        face_nodes = self.Mesh2_face_nodes.data
 
         # Note: x, y, z are np arrays of type float
         # Using np.issubdtype to check if the type is float
@@ -827,9 +698,9 @@ class Grid:
     #     return integral
 
     def to_geodataframe(self,
-                        override=False,
-                        cache=True,
-                        correct_antimeridian_polygons=True):
+                        override: Optional[bool] = False,
+                        cache: Optional[bool] = True,
+                        correct_antimeridian_polygons: Optional[bool] = True):
         """Constructs a ``spatialpandas.GeoDataFrame`` with a "geometry"
         column, containing a collection of Shapely Polygons or MultiPolygons
         representing the geometry of the unstructured grid. Additionally, any
@@ -864,9 +735,9 @@ class Grid:
         return gdf
 
     def to_polycollection(self,
-                          override=False,
-                          cache=True,
-                          correct_antimeridian_polygons=True):
+                          override: Optional[bool] = False,
+                          cache: Optional[bool] = True,
+                          correct_antimeridian_polygons: Optional[bool] = True):
         """Constructs a ``matplotlib.collections.PolyCollection`` object with
         polygons representing the geometry of the unstructured grid, with
         polygons that cross the antimeridian split.
@@ -899,7 +770,9 @@ class Grid:
 
         return poly_collection, corrected_to_original_faces
 
-    def to_shapely_polygons(self, correct_antimeridian_polygons=True):
+    def to_shapely_polygons(self,
+                            correct_antimeridian_polygons: Optional[bool] = True
+                           ):
         """Constructs an array of Shapely Polygons representing each face, with
         antimeridian polygons split according to the GeoJSON standards.
 
diff --git a/uxarray/io/_exodus.py b/uxarray/io/_exodus.py
index 91daed59..ecc243eb 100644
--- a/uxarray/io/_exodus.py
+++ b/uxarray/io/_exodus.py
@@ -5,16 +5,20 @@ from datetime import datetime
 
 from uxarray.grid.connectivity import _replace_fill_values
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.grid.coordinates import _get_lonlat_from_xyz, _get_xyz_from_lonlat
 
 
 # Exodus Number is one-based.
-def _read_exodus(ext_ds, grid_var_names):
+def _read_exodus(ext_ds):
     """Exodus file reader.
 
     Parameters: xarray.Dataset, required
     Returns: ugrid aware xarray.Dataset
     """
 
+    # TODO: UGRID Variable Mapping
+    source_dims_dict = {}
+
     # Not loading specific variables.
     # as there is no way to know number of face types etc. without loading
     # connect1, connect2, connect3, etc..
@@ -50,22 +54,24 @@ def _read_exodus(ext_ds, grid_var_names):
         elif key == "coord":
             ds.Mesh2.attrs['topology_dimension'] = INT_DTYPE(
                 ext_ds.dims['num_dim'])
-            ds["Mesh2_node_x"] = xr.DataArray(data=ext_ds.coord[0],
-                                              dims=["nMesh2_node"],
-                                              attrs={
-                                                  "standard_name": "x",
-                                                  "long_name": "cartesian x",
-                                                  "units": "m",
-                                              })
-            ds["Mesh2_node_y"] = xr.DataArray(data=ext_ds.coord[1],
-                                              dims=["nMesh2_node"],
-                                              attrs={
-                                                  "standard_name": "y",
-                                                  "long_name": "cartesian y",
-                                                  "units": "m",
-                                              })
+            ds["Mesh2_node_cart_x"] = xr.DataArray(
+                data=ext_ds.coord[0],
+                dims=["nMesh2_node"],
+                attrs={
+                    "standard_name": "x",
+                    "long_name": "cartesian x",
+                    "units": "m",
+                })
+            ds["Mesh2_node_cart_y"] = xr.DataArray(
+                data=ext_ds.coord[1],
+                dims=["nMesh2_node"],
+                attrs={
+                    "standard_name": "y",
+                    "long_name": "cartesian y",
+                    "units": "m",
+                })
             if ext_ds.dims['num_dim'] > 2:
-                ds["Mesh2_node_z"] = xr.DataArray(
+                ds["Mesh2_node_cart_z"] = xr.DataArray(
                     data=ext_ds.coord[2],
                     dims=["nMesh2_node"],
                     attrs={
@@ -74,24 +80,26 @@ def _read_exodus(ext_ds, grid_var_names):
                         "units": "m",
                     })
         elif key == "coordx":
-            ds["Mesh2_node_x"] = xr.DataArray(data=ext_ds.coordx,
-                                              dims=["nMesh2_node"],
-                                              attrs={
-                                                  "standard_name": "x",
-                                                  "long_name": "cartesian x",
-                                                  "units": "m",
-                                              })
+            ds["Mesh2_node_cart_x"] = xr.DataArray(
+                data=ext_ds.coordx,
+                dims=["nMesh2_node"],
+                attrs={
+                    "standard_name": "x",
+                    "long_name": "cartesian x",
+                    "units": "m",
+                })
         elif key == "coordy":
-            ds["Mesh2_node_y"] = xr.DataArray(data=ext_ds.coordx,
-                                              dims=["nMesh2_node"],
-                                              attrs={
-                                                  "standard_name": "y",
-                                                  "long_name": "cartesian y",
-                                                  "units": "m",
-                                              })
+            ds["Mesh2_node_cart_y"] = xr.DataArray(
+                data=ext_ds.coordx,
+                dims=["nMesh2_node"],
+                attrs={
+                    "standard_name": "y",
+                    "long_name": "cartesian y",
+                    "units": "m",
+                })
         elif key == "coordz":
             if ext_ds.dims['num_dim'] > 2:
-                ds["Mesh2_node_z"] = xr.DataArray(
+                ds["Mesh2_node_cart_z"] = xr.DataArray(
                     data=ext_ds.coordx,
                     dims=["nMesh2_node"],
                     attrs={
@@ -137,12 +145,36 @@ def _read_exodus(ext_ds, grid_var_names):
                 INT_DTYPE(
                     0)  # NOTE: This might cause an error if numbering has holes
         })
-    print("Finished reading exodus file.")
 
-    return ds
+    # populate lon/lat coordinates
+    lon, lat = _get_lonlat_from_xyz(ds["Mesh2_node_cart_x"].values,
+                                    ds["Mesh2_node_cart_y"].values,
+                                    ds["Mesh2_node_cart_z"].values)
+
+    # populate dataset
+    ds["Mesh2_node_x"] = xr.DataArray(
+        data=lon,
+        dims=["nMesh2_node"],
+        attrs={
+            "standard_name": "longitude",
+            "long_name": "longitude of mesh nodes",
+            "units": "degrees_east",
+        })
+    ds["Mesh2_node_y"] = xr.DataArray(data=lat,
+                                      dims=["nMesh2_node"],
+                                      attrs={
+                                          "standard_name": "latitude",
+                                          "long_name": "latitude of mesh nodes",
+                                          "units": "degrees_north",
+                                      })
+
+    # set lon/lat coordinates
+    ds = ds.set_coords(["Mesh2_node_x", "Mesh2_node_y"])
+
+    return ds, source_dims_dict
 
 
-def _encode_exodus(ds, grid_var_names, outfile=None):
+def _encode_exodus(ds, outfile=None):
     """Encodes an Exodus file.
 
     Parameters
@@ -203,12 +235,10 @@ def _encode_exodus(ds, grid_var_names, outfile=None):
 
     # Set the dim to 3 as we will always have x/y/z for cartesian grid
     # Note: Don't get orig dimension from Mesh2 attribute topology dimension
-    dim = 3
-
-    if "Mesh2_node_cart_z" not in ds.keys():
-        raise RuntimeError(
-            "The grid must have x/y/z cartesian coordinates to encode to exodus."
-        )
+    if "Mesh2_node_cart_x" not in ds:
+        x, y, z = _get_xyz_from_lonlat(ds["Mesh2_node_x"].values,
+                                       ds["Mesh2_node_y"].values)
+        c_data = xr.DataArray([x, y, z])
     else:
         c_data = xr.DataArray([
             ds["Mesh2_node_cart_x"].data.tolist(),
@@ -220,13 +250,12 @@ def _encode_exodus(ds, grid_var_names, outfile=None):
 
     # process face nodes, this array holds num faces at corresponding location
     # eg num_el_all_blks = [0, 0, 6, 12] signifies 6 TRI and 12 SHELL elements
-    num_el_all_blks = np.zeros(ds[grid_var_names["nMaxMesh2_face_nodes"]].size,
-                               "i8")
+    num_el_all_blks = np.zeros(ds["nMaxMesh2_face_nodes"].size, "i8")
     # this list stores connectivity without filling
     conn_nofill = []
 
     # store the number of faces in an array
-    for row in ds[grid_var_names["Mesh2_face_nodes"]].astype(INT_DTYPE).data:
+    for row in ds["Mesh2_face_nodes"].astype(INT_DTYPE).data:
 
         # find out -1 in each row, this indicates lower than max face nodes
         arr = np.where(row == -1)
@@ -242,8 +271,7 @@ def _encode_exodus(ds, grid_var_names, outfile=None):
             conn_nofill.append(list_node)
         elif arr[0].size == 0:
             # increment the number of faces for this "nMaxMesh2_face_nodes" face
-            num_el_all_blks[ds[grid_var_names["nMaxMesh2_face_nodes"]].size -
-                            1] += 1
+            num_el_all_blks[ds["nMaxMesh2_face_nodes"].size - 1] += 1
             # get integer list nodes
             list_node = list(map(int, row.tolist()))
             conn_nofill.append(list_node)
@@ -321,12 +349,7 @@ def _encode_exodus(ds, grid_var_names, outfile=None):
     eb_names = np.empty(num_blks, dtype='str')
     exo_ds["eb_names"] = xr.DataArray(data=xr.DataArray(eb_names),
                                       dims=["num_el_blk"])
-    if dim == 2:
-        cnames = ["x", "y"]
-    elif dim == 3:
-        cnames = ["x", "y", "z"]
-    else:
-        raise RuntimeError("dim must be 2 or 3")
+    cnames = ["x", "y", "z"]
 
     exo_ds["coor_names"] = xr.DataArray(data=xr.DataArray(
         np.array(cnames, dtype='str')),
diff --git a/uxarray/io/_mpas.py b/uxarray/io/_mpas.py
index fb0cb9cd..b7d5035c 100644
--- a/uxarray/io/_mpas.py
+++ b/uxarray/io/_mpas.py
@@ -17,7 +17,9 @@ def _primal_to_ugrid(in_ds, out_ds):
         conventions
     """
 
-    # set mesh topology
+    source_dims_dict = {}
+
+    # set mesh topologys
     out_ds["Mesh2"] = xr.DataArray(
         attrs={
             "cf_role": "mesh_topology",
@@ -116,6 +118,14 @@ def _primal_to_ugrid(in_ds, out_ds):
     # set global attributes
     _set_global_attrs(in_ds, out_ds)
 
+    # populate source dims
+    source_dims_dict['nVertices'] = 'nMesh2_node'
+    source_dims_dict[in_ds['verticesOnCell'].dims[0]] = 'nMesh2_face'
+    source_dims_dict[in_ds['verticesOnCell'].dims[1]] = 'nMaxMesh2_face_nodes'
+    source_dims_dict[in_ds['verticesOnEdge'].dims[0]] = "nMesh2_edge"
+
+    return source_dims_dict
+
 
 def _dual_to_ugrid(in_ds, out_ds):
     """Encodes the MPAS Dual-Mesh in the UGRID conventions.
@@ -129,6 +139,8 @@ def _dual_to_ugrid(in_ds, out_ds):
         conventions
     """
 
+    source_dims_dict = {}
+
     # set mesh topology
     out_ds["Mesh2"] = xr.DataArray(
         attrs={
@@ -223,6 +235,14 @@ def _dual_to_ugrid(in_ds, out_ds):
     # set global attributes
     _set_global_attrs(in_ds, out_ds)
 
+    # populate source dims
+    source_dims_dict[in_ds['latCell'].dims[0]] = "nMesh2_node"
+    source_dims_dict[in_ds['cellsOnVertex'].dims[0]] = "nMesh2_face"
+    source_dims_dict[in_ds['cellsOnVertex'].dims[1]] = "nMaxMesh2_face_nodes"
+    source_dims_dict[in_ds['cellsOnEdge'].dims[0]] = "nMesh2_edge"
+
+    return source_dims_dict
+
 
 def _set_global_attrs(in_ds, out_ds):
     """Helper to set MPAS global attributes.
@@ -375,9 +395,9 @@ def _read_mpas(ext_ds, use_dual=False):
 
     # convert dual-mesh to UGRID
     if use_dual:
-        _dual_to_ugrid(ext_ds, ds)
+        source_dim_map = _dual_to_ugrid(ext_ds, ds)
     # convert primal-mesh to UGRID
     else:
-        _primal_to_ugrid(ext_ds, ds)
+        source_dim_map = _primal_to_ugrid(ext_ds, ds)
 
-    return ds
+    return ds, source_dim_map
diff --git a/uxarray/io/_scrip.py b/uxarray/io/_scrip.py
index ac3b6a66..ca52c4e7 100644
--- a/uxarray/io/_scrip.py
+++ b/uxarray/io/_scrip.py
@@ -20,6 +20,8 @@ def _to_ugrid(in_ds, out_ds):
         to store reassigned SCRIP variables in UGRID conventions
     """
 
+    source_dims_dict = {}
+
     if in_ds['grid_area'].all():
 
         # Create Mesh2_node_x/y variables from grid_corner_lat/lon
@@ -93,7 +95,10 @@ def _to_ugrid(in_ds, out_ds):
     else:
         raise Exception("Structured scrip files are not yet supported")
 
-    return out_ds
+    # populate source dims
+    source_dims_dict[in_ds['grid_center_lon'].dims[0]] = "nMesh2_face"
+
+    return source_dims_dict
 
 
 def _read_scrip(ext_ds):
@@ -122,7 +127,7 @@ def _read_scrip(ext_ds):
 
     try:
         # If not ugrid compliant, translates scrip to ugrid conventions
-        _to_ugrid(ext_ds, ds)
+        source_dims_dict = _to_ugrid(ext_ds, ds)
 
         # Add necessary UGRID attributes to new dataset
         ds["Mesh2"] = xr.DataArray(
@@ -142,7 +147,7 @@ def _read_scrip(ext_ds):
             "https://earthsystemmodeling.org/docs/release/ESMF_6_2_0/ESMF_refdoc/node3.html#SECTION03024000000000000000",
             "for more information on SCRIP Grid file formatting")
 
-    return ds
+    return ds, source_dims_dict
 
 
 def _encode_scrip(mesh2_face_nodes, mesh2_node_x, mesh2_node_y, face_areas):
diff --git a/uxarray/io/_ugrid.py b/uxarray/io/_ugrid.py
index d0e21845..0b6e67b6 100644
--- a/uxarray/io/_ugrid.py
+++ b/uxarray/io/_ugrid.py
@@ -4,59 +4,59 @@ 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):
+def _read_ugrid(xr_ds):
     """UGRID file reader.
 
     Parameters: xarray.Dataset, required
     Returns: ugrid aware xarray.Dataset
     """
 
+    source_dims_dict = {}
     # TODO: obtain and change to Mesh2 construct, see Issue #27
 
     # get the data variable name that has attribute "cf_role" set to "mesh_topology"
     # this is the base xarray.DataArray name
     base_xarray_var = list(
         xr_ds.filter_by_attrs(cf_role="mesh_topology").keys())[0]
+    # TODO: Allow for parsing datasets with more than just coordinates and face nodes
 
-    var_names_dict["Mesh2"] = base_xarray_var
-
-    coord_names = xr_ds[base_xarray_var].node_coordinates.split()
+    xr_ds = xr_ds.rename({base_xarray_var: "Mesh2"})
 
+    # map and rename coordinates
+    coord_names = xr_ds["Mesh2"].node_coordinates.split()
     if len(coord_names) == 1:
-        var_names_dict["Mesh2_node_x"] = coord_names[0]
+        xr_ds = xr_ds.rename({coord_names[0]: "Mesh2_node_x"})
     elif len(coord_names) == 2:
-        var_names_dict["Mesh2_node_x"] = coord_names[0]
-        var_names_dict["Mesh2_node_y"] = coord_names[1]
-    elif len(coord_names) == 3:
-        var_names_dict["Mesh2_node_x"] = coord_names[0]
-        var_names_dict["Mesh2_node_y"] = coord_names[1]
-        var_names_dict["Mesh2_node_z"] = coord_names[2]
+        xr_ds = xr_ds.rename({
+            coord_names[0]: "Mesh2_node_x",
+            coord_names[1]: "Mesh2_node_y"
+        })
+    # map and rename dimensions
+    coord_dim_name = xr_ds["Mesh2_node_x"].dims
 
-    # set #nodes use x coordinates, y or z will be the same and can also be used
-    coord_dim_name = xr_ds[var_names_dict["Mesh2_node_x"]].dims
-    var_names_dict["nMesh2_node"] = coord_dim_name[0]
+    xr_ds = xr_ds.rename({coord_dim_name[0]: "nMesh2_node"})
 
-    face_node_names = xr_ds[base_xarray_var].face_node_connectivity.split()
+    face_node_names = xr_ds["Mesh2"].face_node_connectivity.split()
 
     face_node_name = face_node_names[0]
-    var_names_dict["Mesh2_face_nodes"] = xr_ds[face_node_name].name
-    var_names_dict["nMesh2_face"] = xr_ds[face_node_name].dims[0]
-    var_names_dict["nMaxMesh2_face_nodes"] = xr_ds[face_node_name].dims[1]
-
-    if len(coord_names) == 2:
-        # set coordinates
-        xr_ds = xr_ds.set_coords(
-            [var_names_dict["Mesh2_node_x"], var_names_dict["Mesh2_node_y"]])
-    else:
-        xr_ds = xr_ds.set_coords([
-            var_names_dict["Mesh2_node_x"], var_names_dict["Mesh2_node_y"],
-            var_names_dict["Mesh2_node_z"]
-        ])
+    xr_ds = xr_ds.rename({xr_ds[face_node_name].name: "Mesh2_face_nodes"})
+
+    xr_ds = xr_ds.rename({
+        xr_ds["Mesh2_face_nodes"].dims[0]: "nMesh2_face",
+        xr_ds["Mesh2_face_nodes"].dims[1]: "nMaxMesh2_face_nodes"
+    })
+
+    xr_ds = xr_ds.set_coords(["Mesh2_node_x", "Mesh2_node_y"])
 
     # standardize fill values and data type for face nodes
-    xr_ds = _standardize_fill_values(xr_ds, var_names_dict)
+    xr_ds = _standardize_fill_values(xr_ds)
 
-    return xr_ds, var_names_dict
+    # populate source dimensions
+    source_dims_dict[coord_dim_name[0]] = "nMesh2_node"
+    source_dims_dict[xr_ds["Mesh2_face_nodes"].dims[0]] = "nMesh2_face"
+    source_dims_dict[xr_ds["Mesh2_face_nodes"].dims[1]] = "nMaxMesh2_face_nodes"
+
+    return xr_ds, source_dims_dict
 
 
 def _encode_ugrid(ds):
@@ -71,15 +71,13 @@ def _encode_ugrid(ds):
     return ds
 
 
-def _standardize_fill_values(ds, var_names_dict):
+def _standardize_fill_values(ds):
     """Standardizes the fill values and data type of index variables.
 
     Parameters
     ----------
     ds : xarray.Dataset
         Input Dataset
-    ds_vars_dict : dict
-        UGRID Variable Name Dictionary
 
     Returns
     ----------
@@ -88,12 +86,12 @@ def _standardize_fill_values(ds, var_names_dict):
     """
 
     # original face nodes
-    face_nodes = ds[var_names_dict['Mesh2_face_nodes']].values
+    face_nodes = ds["Mesh2_face_nodes"].values
 
     # original fill value, if one exists
-    if "_FillValue" in ds[var_names_dict['Mesh2_face_nodes']].attrs:
-        original_fv = ds[var_names_dict['Mesh2_face_nodes']]._FillValue
-    elif np.isnan(ds[var_names_dict['Mesh2_face_nodes']].values).any():
+    if "_FillValue" in ds["Mesh2_face_nodes"].attrs:
+        original_fv = ds["Mesh2_face_nodes"]._FillValue
+    elif np.isnan(ds["Mesh2_face_nodes"].values).any():
         original_fv = np.nan
     else:
         original_fv = None
@@ -106,11 +104,10 @@ def _standardize_fill_values(ds, var_names_dict):
                                               new_fill=INT_FILL_VALUE,
                                               new_dtype=INT_DTYPE)
         # reassign data to use updated face nodes
-        ds[var_names_dict['Mesh2_face_nodes']].data = new_face_nodes
+        ds["Mesh2_face_nodes"].data = new_face_nodes
 
         # use new fill value
-        ds[var_names_dict['Mesh2_face_nodes']].attrs[
-            '_FillValue'] = INT_FILL_VALUE
+        ds["Mesh2_face_nodes"].attrs['_FillValue'] = INT_FILL_VALUE
 
     return ds
 
@@ -129,3 +126,12 @@ def _is_ugrid(ds):
         return True
     else:
         return False
+
+
+def _validate_minimum_ugrid(grid_ds):
+    """Checks whether a given ``grid_ds`` meets the requirements for a minimum
+    unstructured grid encoded in the UGRID conventions, containing a set of (x,
+    y) latlon coordinates and face node connectivity."""
+    return ("Mesh2_node_x" in grid_ds and "Mesh2_node_y" in grid_ds) or (
+        "Mesh2_node_cart_x" in grid_ds and "Mesh2_node_cart_y" in grid_ds and
+        "Mesh2_node_cart_z" in grid_ds) and "Mesh2_face_nodes" in grid_ds
diff --git a/uxarray/io/_vertices.py b/uxarray/io/_vertices.py
new file mode 100644
index 00000000..ff32a989
--- /dev/null
+++ b/uxarray/io/_vertices.py
@@ -0,0 +1,113 @@
+import xarray as xr
+import numpy as np
+
+from uxarray.constants import INT_FILL_VALUE, INT_DTYPE
+
+
+def _read_face_vertices(face_vertices, latlon):
+    """Create a grid with faces constructed from vertices specified by the
+    given argument.
+
+    Parameters
+    ----------
+    dataset : ndarray, list, tuple, required
+        Input vertex coordinates that form our face(s)
+    """
+    grid_ds = xr.Dataset()
+    grid_ds["Mesh2"] = xr.DataArray(
+        attrs={
+            "cf_role": "mesh_topology",
+            "long_name": "Topology data of unstructured mesh",
+            "topology_dimension": -1,
+            "node_coordinates": "Mesh2_node_x Mesh2_node_y Mesh2_node_z",
+            "node_dimension": "nMesh2_node",
+            "face_node_connectivity": "Mesh2_face_nodes",
+            "face_dimension": "nMesh2_face"
+        })
+
+    grid_ds.Mesh2.attrs['topology_dimension'] = face_vertices.ndim
+
+    if latlon:
+        x_units = "degrees_east"
+        y_units = "degrees_north"
+    else:
+        x_units = 'm'
+        y_units = 'm'
+        z_units = 'm'
+
+    x_coord = face_vertices[:, :, 0].flatten()
+    y_coord = face_vertices[:, :, 1].flatten()
+
+    if face_vertices[0][0].size > 2:
+        z_coord = face_vertices[:, :, 2].flatten()
+    else:
+        z_coord = x_coord * 0.0
+
+    # Identify unique vertices and their indices
+    unique_verts, indices = np.unique(face_vertices.reshape(
+        -1, face_vertices.shape[-1]),
+                                      axis=0,
+                                      return_inverse=True)
+
+    # Nodes index that contain a fill value
+    fill_value_mask = np.logical_or(unique_verts[:, 0] == INT_FILL_VALUE,
+                                    unique_verts[:, 1] == INT_FILL_VALUE)
+    if face_vertices[0][0].size > 2:
+        fill_value_mask = np.logical_or(unique_verts[:, 0] == INT_FILL_VALUE,
+                                        unique_verts[:, 1] == INT_FILL_VALUE,
+                                        unique_verts[:, 2] == INT_FILL_VALUE)
+
+    # Get the indices of all the False values in fill_value_mask
+    false_indices = np.where(fill_value_mask == True)[0]
+
+    # Check if any False values were found
+    indices = indices.astype(INT_DTYPE)
+    if false_indices.size > 0:
+
+        # Remove the rows corresponding to False values in unique_verts
+        unique_verts = np.delete(unique_verts, false_indices, axis=0)
+
+        # Update indices accordingly
+        for i, idx in enumerate(false_indices):
+            indices[indices == idx] = INT_FILL_VALUE
+            indices[(indices > idx) & (indices != INT_FILL_VALUE)] -= 1
+
+    if latlon:
+        grid_ds["Mesh2_node_x"] = xr.DataArray(data=unique_verts[:, 0],
+                                               dims=["nMesh2_node"],
+                                               attrs={"units": x_units})
+        grid_ds["Mesh2_node_y"] = xr.DataArray(data=unique_verts[:, 1],
+                                               dims=["nMesh2_node"],
+                                               attrs={"units": y_units})
+    else:
+        grid_ds["Mesh2_node_cart_x"] = xr.DataArray(data=unique_verts[:, 0],
+                                                    dims=["nMesh2_node"],
+                                                    attrs={"units": x_units})
+
+        grid_ds["Mesh2_node_cart_y"] = xr.DataArray(data=unique_verts[:, 1],
+                                                    dims=["nMesh2_node"],
+                                                    attrs={"units": y_units})
+
+        if face_vertices.shape[-1] > 2:
+            grid_ds["Mesh2_node_cart_z"] = xr.DataArray(
+                data=unique_verts[:, 2],
+                dims=["nMesh2_node"],
+                attrs={"units": z_units})
+        else:
+            grid_ds["Mesh2_node_cart_z"] = xr.DataArray(
+                data=unique_verts[:, 1] * 0.0,
+                dims=["nMesh2_node"],
+                attrs={"units": z_units})
+
+    # Create connectivity array using indices of unique vertices
+    connectivity = indices.reshape(face_vertices.shape[:-1])
+    grid_ds["Mesh2_face_nodes"] = xr.DataArray(
+        data=xr.DataArray(connectivity).astype(INT_DTYPE),
+        dims=["nMesh2_face", "nMaxMesh2_face_nodes"],
+        attrs={
+            "cf_role": "face_node_connectivity",
+            "_FillValue": INT_FILL_VALUE,
+            "start_index": 0
+        })
+
+    return grid_ds
diff --git a/uxarray/io/utils.py b/uxarray/io/utils.py
index ce531748..cfd05173 100644
--- a/uxarray/io/utils.py
+++ b/uxarray/io/utils.py
@@ -24,17 +24,17 @@ def _parse_grid_type(dataset):
     """
     # exodus with coord or coordx
     if "coord" in dataset:
-        mesh_type = "exo"
+        mesh_type = "Exodus"
     elif "coordx" in dataset:
-        mesh_type = "exo"
+        mesh_type = "Exodus"
     # scrip with grid_center_lon
     elif "grid_center_lon" in dataset:
-        mesh_type = "scrip"
+        mesh_type = "Scrip"
     # ugrid topology
     elif _is_ugrid(dataset):
-        mesh_type = "ugrid"
+        mesh_type = "UGRID"
     elif "verticesOnCell" in dataset:
-        mesh_type = "mpas"
+        mesh_type = "MPAS"
     else:
         raise RuntimeError(f"Could not recognize dataset format.")
     return mesh_type


Test patch:
diff --git a/test/test_api.py b/test/test_api.py
index 87955107..b9bb6c8d 100644
--- a/test/test_api.py
+++ b/test/test_api.py
@@ -53,8 +53,6 @@ class TestAPI(TestCase):
                         constants.NNODES_outCSne30)
         nt.assert_equal(len(uxds_var2_ne30.uxgrid._ds.data_vars),
                         constants.DATAVARS_outCSne30)
-        nt.assert_equal(uxds_var2_ne30.uxgrid.source_grid,
-                        str(self.gridfile_ne30))
         nt.assert_equal(uxds_var2_ne30.source_datasets,
                         str(self.dsfile_var2_ne30))
 
@@ -69,8 +67,7 @@ class TestAPI(TestCase):
                         constants.NNODES_outCSne30)
         nt.assert_equal(len(uxds_mf_ne30.uxgrid._ds.data_vars),
                         constants.DATAVARS_outCSne30)
-        nt.assert_equal(uxds_mf_ne30.uxgrid.source_grid,
-                        str(self.gridfile_ne30))
+
         nt.assert_equal(uxds_mf_ne30.source_datasets, self.dsfiles_mf_ne30)
 
     def test_open_grid(self):
diff --git a/test/test_dataset.py b/test/test_dataset.py
index a8b79b28..81b43908 100644
--- a/test/test_dataset.py
+++ b/test/test_dataset.py
@@ -2,6 +2,7 @@ import os
 from unittest import TestCase
 from pathlib import Path
 import numpy.testing as nt
+import xarray as xr
 
 import uxarray as ux
 
@@ -20,6 +21,8 @@ dsfiles_mf_ne30 = str(
 gridfile_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
 dsfile_v1_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "v1.nc"
 
+mpas_ds_path = current_path / 'meshfiles' / "mpas" / "QU" / 'mesh.QU.1920km.151026.nc'
+
 
 class TestUxDataset(TestCase):
 
@@ -53,3 +56,13 @@ class TestUxDataset(TestCase):
                 uxds_var2_geoflow.info(show_attrs=True)
             except Exception as exc:
                 assert False, f"'uxds_var2_geoflow.info()' raised an exception: {exc}"
+
+    def test_ugrid_dim_names(self):
+        """Tests the remapping of dimensions to the UGRID conventions."""
+
+        ugrid_dims = ["nMesh2_face", "nMesh2_node", "nMesh2_edge"]
+
+        uxds_remap = ux.open_dataset(mpas_ds_path, mpas_ds_path)
+
+        for dim in ugrid_dims:
+            assert dim in uxds_remap.dims
diff --git a/test/test_exodus.py b/test/test_exodus.py
index a07bbb65..41125867 100644
--- a/test/test_exodus.py
+++ b/test/test_exodus.py
@@ -19,6 +19,7 @@ class TestExodus(TestCase):
         """Read an exodus file and writes a exodus file."""
 
         uxgrid = ux.open_grid(self.exo_filename)
+        pass
 
     def test_init_verts(self):
         """Create a uxarray grid from vertices and saves a 1 face exodus
@@ -35,8 +36,8 @@ class TestExodus(TestCase):
 
         uxgrid = ux.open_grid(self.exo2_filename)
 
-        uxgrid.encode_as("ugrid")
-        uxgrid.encode_as("exodus")
+        uxgrid.encode_as("UGRID")
+        uxgrid.encode_as("Exodus")
 
     def test_standardized_dtype_and_fill(self):
         """Test to see if Mesh2_Face_Nodes uses the expected integer datatype
diff --git a/test/test_geometry.py b/test/test_geometry.py
index 4d289024..1542f088 100644
--- a/test/test_geometry.py
+++ b/test/test_geometry.py
@@ -27,10 +27,7 @@ class TestAntimeridian:
     def test_crossing(self):
         verts = [[[-170, 40], [180, 30], [165, 25], [-170, 20]]]
 
-        uxgrid = ux.open_grid(verts,
-                              vertices=True,
-                              islatlon=True,
-                              isconcave=False)
+        uxgrid = ux.open_grid(verts, latlon=True)
 
         gdf = uxgrid.to_geodataframe()
 
@@ -41,9 +38,6 @@ class TestAntimeridian:
     def test_point_on(self):
         verts = [[[-170, 40], [180, 30], [-170, 20]]]
 
-        uxgrid = ux.open_grid(verts,
-                              vertices=True,
-                              islatlon=True,
-                              isconcave=False)
+        uxgrid = ux.open_grid(verts, latlon=True)
 
         assert len(uxgrid.antimeridian_face_indices) == 1
diff --git a/test/test_grid.py b/test/test_grid.py
index 36cd51b5..ca55b4f1 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -45,13 +45,13 @@ class TestGrid(TestCase):
         """Reads a ugrid file and encodes it as `xarray.Dataset` in various
         types."""
 
-        self.grid_CSne30.encode_as("ugrid")
-        self.grid_RLL1deg.encode_as("ugrid")
-        self.grid_RLL10deg_CSne4.encode_as("ugrid")
+        self.grid_CSne30.encode_as("UGRID")
+        self.grid_RLL1deg.encode_as("UGRID")
+        self.grid_RLL10deg_CSne4.encode_as("UGRID")
 
-        self.grid_CSne30.encode_as("exodus")
-        self.grid_RLL1deg.encode_as("exodus")
-        self.grid_RLL10deg_CSne4.encode_as("exodus")
+        self.grid_CSne30.encode_as("Exodus")
+        self.grid_RLL1deg.encode_as("Exodus")
+        self.grid_RLL10deg_CSne4.encode_as("Exodus")
 
     def test_open_non_mesh2_write_exodus(self):
         """Loads grid files of different formats using uxarray's open_dataset
@@ -59,7 +59,7 @@ class TestGrid(TestCase):
 
         grid_geoflow = ux.open_grid(gridfile_CSne30)
 
-        exods = grid_geoflow.encode_as("exodus")
+        exods = grid_geoflow.encode_as("Exodus")
         # Remove the _FillValue attribute from the variable's attributes
         if '_FillValue' in grid_geoflow._ds['Mesh2_face_nodes'].attrs:
             del grid_geoflow._ds['Mesh2_face_nodes'].attrs['_FillValue']
@@ -121,42 +121,30 @@ class TestGrid(TestCase):
 
         # Now consturct the grid using the faces_coords
         verts_cart = np.array(faces_coords)
-        vgrid = ux.open_grid(verts_cart,
-                             vertices=True,
-                             islatlon=False,
-                             isconcave=False)
+        vgrid = ux.open_grid(verts_cart, latlon=False)
 
-        assert (vgrid.source_grid == "From vertices")
         assert (vgrid.nMesh2_face == 6)
         assert (vgrid.nMesh2_node == 8)
-        vgrid.encode_as("ugrid")
+        vgrid.encode_as("UGRID")
 
         # Test the case when user created a nested one-face grid
         faces_verts_one = np.array([
             np.array([[150, 10], [160, 20], [150, 30], [135, 30], [125, 20],
                       [135, 10]])
         ])
-        vgrid = ux.open_grid(faces_verts_one,
-                             vertices=True,
-                             islatlon=True,
-                             isconcave=False)
-        assert (vgrid.source_grid == "From vertices")
+        vgrid = ux.open_grid(faces_verts_one, latlon=True)
         assert (vgrid.nMesh2_face == 1)
         assert (vgrid.nMesh2_node == 6)
-        vgrid.encode_as("ugrid")
+        vgrid.encode_as("UGRID")
 
         # Test the case when user created a one-face grid
         faces_verts_single_face = np.array([[150, 10], [160, 20], [150, 30],
                                             [135, 30], [125, 20], [135, 10]])
 
-        vgrid = ux.open_grid(faces_verts_single_face,
-                             vertices=True,
-                             islatlon=True,
-                             isconcave=False)
-        assert (vgrid.source_grid == "From vertices")
+        vgrid = ux.open_grid(faces_verts_single_face, latlon=True)
         assert (vgrid.nMesh2_face == 1)
         assert (vgrid.nMesh2_node == 6)
-        vgrid.encode_as("ugrid")
+        vgrid.encode_as("UGRID")
 
     def test_init_verts_different_input_datatype(self):
         """Create a uxarray grid from multiple face vertices with different
@@ -174,15 +162,10 @@ class TestGrid(TestCase):
             np.array([[95, 10], [105, 20], [100, 30], [85, 30], [75, 20],
                       [85, 10]]),
         ])
-        vgrid = ux.open_grid(faces_verts_ndarray,
-                             vertices=True,
-                             islatlon=True,
-                             isconcave=False)
-
-        assert (vgrid.source_grid == "From vertices")
+        vgrid = ux.open_grid(faces_verts_ndarray, latlon=True)
         assert (vgrid.nMesh2_face == 3)
         assert (vgrid.nMesh2_node == 14)
-        vgrid.encode_as("ugrid")
+        vgrid.encode_as("UGRID")
 
         # Test initializing Grid from list
         faces_verts_list = [[[150, 10], [160, 20], [150, 30], [135, 30],
@@ -191,14 +174,10 @@ class TestGrid(TestCase):
                              [100, 30], [105, 20]],
                             [[95, 10], [105, 20], [100, 30], [85, 30], [75, 20],
                              [85, 10]]]
-        vgrid = ux.open_grid(faces_verts_list,
-                             vertices=True,
-                             islatlon=False,
-                             isconcave=False)
-        assert (vgrid.source_grid == "From vertices")
+        vgrid = ux.open_grid(faces_verts_list, latlon=False)
         assert (vgrid.nMesh2_face == 3)
         assert (vgrid.nMesh2_node == 14)
-        vgrid.encode_as("ugrid")
+        vgrid.encode_as("UGRID")
 
         # Test initializing Grid from tuples
         faces_verts_tuples = [
@@ -206,14 +185,10 @@ class TestGrid(TestCase):
             ((125, 20), (135, 30), (125, 60), (110, 60), (100, 30), (105, 20)),
             ((95, 10), (105, 20), (100, 30), (85, 30), (75, 20), (85, 10))
         ]
-        vgrid = ux.open_grid(faces_verts_tuples,
-                             vertices=True,
-                             islatlon=False,
-                             isconcave=False)
-        assert (vgrid.source_grid == "From vertices")
+        vgrid = ux.open_grid(faces_verts_tuples, latlon=False)
         assert (vgrid.nMesh2_face == 3)
         assert (vgrid.nMesh2_node == 14)
-        vgrid.encode_as("ugrid")
+        vgrid.encode_as("UGRID")
 
     def test_init_verts_fill_values(self):
         faces_verts_filled_values = [[[150, 10], [160, 20], [150, 30],
@@ -224,11 +199,10 @@ class TestGrid(TestCase):
                                      [[95, 10], [105, 20], [100, 30], [85, 30],
                                       [ux.INT_FILL_VALUE, ux.INT_FILL_VALUE],
                                       [ux.INT_FILL_VALUE, ux.INT_FILL_VALUE]]]
-        vgrid = ux.open_grid(faces_verts_filled_values,
-                             vertices=True,
-                             islatlon=False,
-                             isconcave=False)
-        assert (vgrid.source_grid == "From vertices")
+        vgrid = ux.open_grid(
+            faces_verts_filled_values,
+            latlon=False,
+        )
         assert (vgrid.nMesh2_face == 3)
         assert (vgrid.nMesh2_node == 12)
 
@@ -238,16 +212,13 @@ class TestGrid(TestCase):
 
         # Dataset with standard UGRID variable names
         # Coordinates
-        xr.testing.assert_equal(
-            self.grid_CSne30.Mesh2_node_x, self.grid_CSne30._ds[
-                self.grid_CSne30.grid_var_names["Mesh2_node_x"]])
-        xr.testing.assert_equal(
-            self.grid_CSne30.Mesh2_node_y, self.grid_CSne30._ds[
-                self.grid_CSne30.grid_var_names["Mesh2_node_y"]])
+        xr.testing.assert_equal(self.grid_CSne30.Mesh2_node_x,
+                                self.grid_CSne30._ds["Mesh2_node_x"])
+        xr.testing.assert_equal(self.grid_CSne30.Mesh2_node_y,
+                                self.grid_CSne30._ds["Mesh2_node_y"])
         # Variables
-        xr.testing.assert_equal(
-            self.grid_CSne30.Mesh2_face_nodes, self.grid_CSne30._ds[
-                self.grid_CSne30.grid_var_names["Mesh2_face_nodes"]])
+        xr.testing.assert_equal(self.grid_CSne30.Mesh2_face_nodes,
+                                self.grid_CSne30._ds["Mesh2_face_nodes"])
 
         # Dimensions
         n_nodes = self.grid_CSne30.Mesh2_node_x.shape[0]
@@ -267,16 +238,13 @@ class TestGrid(TestCase):
         # Dataset with non-standard UGRID variable names
         grid_geoflow = ux.open_grid(gridfile_geoflow)
 
-        xr.testing.assert_equal(
-            grid_geoflow.Mesh2_node_x,
-            grid_geoflow._ds[grid_geoflow.grid_var_names["Mesh2_node_x"]])
-        xr.testing.assert_equal(
-            grid_geoflow.Mesh2_node_y,
-            grid_geoflow._ds[grid_geoflow.grid_var_names["Mesh2_node_y"]])
+        xr.testing.assert_equal(grid_geoflow.Mesh2_node_x,
+                                grid_geoflow._ds["Mesh2_node_x"])
+        xr.testing.assert_equal(grid_geoflow.Mesh2_node_y,
+                                grid_geoflow._ds["Mesh2_node_y"])
         # Variables
-        xr.testing.assert_equal(
-            grid_geoflow.Mesh2_face_nodes,
-            grid_geoflow._ds[grid_geoflow.grid_var_names["Mesh2_face_nodes"]])
+        xr.testing.assert_equal(grid_geoflow.Mesh2_face_nodes,
+                                grid_geoflow._ds["Mesh2_face_nodes"])
         # Dimensions
         n_nodes = grid_geoflow.Mesh2_node_x.shape[0]
         n_faces, n_face_nodes = grid_geoflow.Mesh2_face_nodes.shape
@@ -339,7 +307,8 @@ class TestFaceAreas(TestCase):
     def test_calculate_total_face_area_file(self):
         """Create a uxarray grid from vertices and saves an exodus file."""
 
-        area = self.grid_CSne30.calculate_total_face_area()
+        # = self.grid_CSne30.calculate_total_face_area()
+        area = ux.open_grid(gridfile_CSne30).calculate_total_face_area()
 
         nt.assert_almost_equal(area, constants.MESH30_AREA, decimal=3)
 
@@ -398,10 +367,7 @@ class TestFaceAreas(TestCase):
                       [75, 20, 0], [85, 10, 0]]),
         ])
         # load our vertices into a UXarray Grid object
-        verts_grid = ux.open_grid(faces_verts_ndarray,
-                                  vertices=True,
-                                  islatlon=True,
-                                  isconcave=False)
+        verts_grid = ux.open_grid(faces_verts_ndarray, latlon=True)
 
         face_verts_areas = verts_grid.face_areas
 
@@ -440,17 +406,17 @@ class TestPopulateCoordinates(TestCase):
 
         verts_degree = np.stack((lon_deg, lat_deg), axis=1)
 
-        vgrid = ux.open_grid(verts_degree, islatlon=False)
-        _populate_cartesian_xyz_coord(vgrid)
+        vgrid = ux.open_grid(verts_degree, latlon=True)
+        #_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],
+            nt.assert_almost_equal(vgrid.Mesh2_node_cart_x.values[i],
                                    cart_x[i],
                                    decimal=12)
-            nt.assert_almost_equal(vgrid._ds["Mesh2_node_cart_y"].values[i],
+            nt.assert_almost_equal(vgrid.Mesh2_node_cart_y.values[i],
                                    cart_y[i],
                                    decimal=12)
-            nt.assert_almost_equal(vgrid._ds["Mesh2_node_cart_z"].values[i],
+            nt.assert_almost_equal(vgrid.Mesh2_node_cart_z.values[i],
                                    cart_z[i],
                                    decimal=12)
 
@@ -481,7 +447,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 = ux.open_grid(verts_cart, latlon=False)
         _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))))
@@ -812,10 +778,10 @@ class TestConnectivity(TestCase):
         # Now we don't need the face_nodes_conn_index anymore.
         del face_nodes_conn_index
 
-        vgrid = ux.Grid(face_nodes_conn_lonlat,
-                        vertices=True,
-                        islatlon=True,
-                        concave=False)
+        vgrid = ux.Grid.from_face_vertices(
+            face_nodes_conn_lonlat,
+            latlon=True,
+        )
 
         # We eyeballed the `Grid._face_nodes_connectivity` and wrote the following expected result
         expected = np.array([
@@ -838,7 +804,7 @@ class TestConnectivity(TestCase):
 
         for grid_path in grid_paths:
             grid_xr = xr.open_dataset(grid_path)
-            grid_ux = ux.Grid(grid_xr)
+            grid_ux = ux.Grid.from_dataset(grid_xr)
 
             # use the dictionary method to build the node_face_connectivity
             node_face_connectivity = {}
@@ -866,6 +832,45 @@ class TestConnectivity(TestCase):
                                    face_index_from_dict))
 
 
+class TestClassMethods(TestCase):
+
+    gridfile_ugrid = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
+    gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
+    gridfile_exodus = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g"
+    gridfile_scrip = current_path / "meshfiles" / "scrip" / "outCSne8" / "outCSne8.nc"
+
+    def test_from_dataset(self):
+
+        # UGRID
+        xrds = xr.open_dataset(self.gridfile_ugrid)
+        uxgrid = ux.Grid.from_dataset(xrds)
+
+        # MPAS
+        xrds = xr.open_dataset(self.gridfile_mpas)
+        uxgrid = ux.Grid.from_dataset(xrds, use_dual=False)
+        uxgrid = ux.Grid.from_dataset(xrds, use_dual=True)
+
+        # Exodus
+        xrds = xr.open_dataset(self.gridfile_exodus)
+        uxgrid = ux.Grid.from_dataset(xrds)
+
+        # SCRIP
+        xrds = xr.open_dataset(self.gridfile_scrip)
+        uxgrid = ux.Grid.from_dataset(xrds)
+
+        pass
+
+    def test_from_face_vertices(self):
+        single_face_latlon = [(0.0, 90.0), (-180, 0.0), (0.0, -90)]
+        uxgrid = ux.Grid.from_face_vertices(single_face_latlon, latlon=True)
+
+        multi_face_latlon = [[(0.0, 90.0), (-180, 0.0), (0.0, -90)],
+                             [(0.0, 90.0), (180, 0.0), (0.0, -90)]]
+        uxgrid = ux.Grid.from_face_vertices(multi_face_latlon, latlon=True)
+
+        single_face_cart = [(0.0,)]
+
+
 class TestBallTree(TestCase):
 
     corner_grid_files = [gridfile_CSne30, gridfile_mpas]
diff --git a/test/test_mpas.py b/test/test_mpas.py
index 2c401d6b..cb0b9c79 100644
--- a/test/test_mpas.py
+++ b/test/test_mpas.py
@@ -24,8 +24,8 @@ class TestMPAS(TestCase):
 
     def test_read_mpas(self):
         """Tests execution of _read_mpas()"""
-        mpas_primal_ugrid = _read_mpas(self.mpas_xr_ds, use_dual=False)
-        mpas_dual_ugrid = _read_mpas(self.mpas_xr_ds, use_dual=True)
+        mpas_primal_ugrid, _ = _read_mpas(self.mpas_xr_ds, use_dual=False)
+        mpas_dual_ugrid, _ = _read_mpas(self.mpas_xr_ds, use_dual=True)
 
     def test_mpas_to_grid(self):
         """Tests creation of Grid object from converted MPAS dataset."""
@@ -36,7 +36,7 @@ class TestMPAS(TestCase):
         """Verifies that the Primal-Mesh was converted properly."""
 
         # primal-mesh encoded in the UGRID conventions
-        ds = _read_mpas(self.mpas_xr_ds, use_dual=False)
+        ds, _ = _read_mpas(self.mpas_xr_ds, use_dual=False)
 
         # check for correct dimensions
         expected_ugrid_dims = [
@@ -59,7 +59,7 @@ class TestMPAS(TestCase):
         """Verifies that the Dual-Mesh was converted properly."""
 
         # dual-mesh encoded in the UGRID conventions
-        ds = _read_mpas(self.mpas_xr_ds, use_dual=True)
+        ds, _ = _read_mpas(self.mpas_xr_ds, use_dual=True)
 
         # check for correct dimensions
         expected_ugrid_dims = [
@@ -110,7 +110,7 @@ class TestMPAS(TestCase):
         ]
 
         # included attrs: 'sphere_radius', 'mesh_spec' 'on_a_sphere'
-        ds = _read_mpas(self.mpas_xr_ds)
+        ds, _ = _read_mpas(self.mpas_xr_ds)
 
         # set dummy attrs to test execution
         ds.attrs['mesh_id'] = "12345678"
diff --git a/test/test_scrip.py b/test/test_scrip.py
index ab48296f..d85f3863 100644
--- a/test/test_scrip.py
+++ b/test/test_scrip.py
@@ -63,7 +63,7 @@ class TestUgrid(TestCase):
         """Read an Exodus dataset and encode that as a UGRID format."""
 
         ux_grid = ux.open_grid(gridfile_exo_ne8)
-        ux_grid.encode_as("ugrid")
+        ux_grid.encode_as("UGRID")
 
     def test_standardized_dtype_and_fill(self):
         """Test to see if Mesh2_Face_Nodes uses the expected integer datatype
diff --git a/test/test_ugrid.py b/test/test_ugrid.py
index ab48296f..d85f3863 100644
--- a/test/test_ugrid.py
+++ b/test/test_ugrid.py
@@ -63,7 +63,7 @@ class TestUgrid(TestCase):
         """Read an Exodus dataset and encode that as a UGRID format."""
 
         ux_grid = ux.open_grid(gridfile_exo_ne8)
-        ux_grid.encode_as("ugrid")
+        ux_grid.encode_as("UGRID")
 
     def test_standardized_dtype_and_fill(self):
         """Test to see if Mesh2_Face_Nodes uses the expected integer datatype