Working with MPAS Grids#
Authors: Philip Chmielowiec, Ian Franda
This notebook showcases how to work with datasets from the Model for Prediction Across Scales (MPAS).
MPAS Grid Overview#
The defining feature of MPAS when compared to other models is that its unstructured grid is composed of Voronoi Meshes with a C-grid staggering. This means that the grid can be broken down into two meshes: Primal and Dual. The Primal Mesh is composed of Voronoi regions and the Dual Mesh is composed of Delaunay Triangles. The figure below showcases this relationship, with the dashed triangles being the dual of the Voronoi mesh.
Since the Primal Mesh is predominantly made up of hexagons, with pentagons and heptagons occasionally being present, and the Dual Mesh being strictly triangular, this notebook will showcase how we can represent both of these meshes in the UGRID conventions using UXarray.
Dataset Overview#
Before diving straight into using UXarray, it is important to first investigate how MPAS datasets are represented.
As mentioned in earlier notebooks, the grid definition and data variables are typically stored in separate files. However, in this example, our dataset will contain both within the same file, which is often the case when working with smaller datasets. Additionally, even when working with separate Grid and Data files in MPAS, the definition of the Primal and Dual mesh are still stored under the same Grid file.
Below we can take a quick look into the dataset by opening it with Xarray.
# Dataset Path
mpas_root_filepath = "../../test/meshfiles/mpas/"
mpas_dataset_filepath = mpas_root_filepath + "QU/mesh.QU.1920km.151026.nc"
import xarray as xr
xrds_mpas = xr.open_dataset(mpas_dataset_filepath)
xrds_mpas
<xarray.Dataset> Dimensions: (nCells: 162, nEdges: 480, nVertices: 320, maxEdges: 6, maxEdges2: 12, TWO: 2, vertexDegree: 3) Dimensions without coordinates: nCells, nEdges, nVertices, maxEdges, maxEdges2, TWO, vertexDegree Data variables: (12/42) latCell (nCells) float64 ... lonCell (nCells) float64 ... xCell (nCells) float64 ... yCell (nCells) float64 ... zCell (nCells) float64 ... indexToCellID (nCells) int32 ... ... ... cellQuality (nCells) float64 ... gridSpacing (nCells) float64 ... triangleQuality (nVertices) float64 ... triangleAngleQuality (nVertices) float64 ... obtuseTriangle (nVertices) int32 ... meshDensity (nCells) float64 ... Attributes: on_a_sphere: YES sphere_radius: 1.0 is_periodic: NO history: MpasMeshConverter.x base_grids/x1.162.grid.nc base_meshes... mesh_spec: 1.0 Conventions: MPAS source: MpasMeshConverter.x file_id: rku96q0z66
Here we opened up the dataset to get an overview of the full set of grid variables needed to describe an MPAS grid as outlined in the MPAS Specification Document [2]. Below is a list of the key grid variables that are used for representing and constructing the Primal and Dual meshes.
Primal Mesh#
lonVertex, latVertex: Corner Vertices of Primal Mesh cells
lonCell, latCell: Center Vertices of Primal Mesh cells
verticesOnCell: Vertex indices that surround each Primal Mesh cell
verticesOnEdge: Vertex indices that saddle a given edge
nEdgesOnCell: Maximum number of edges that can surround a cell
Dual Mesh#
lonCell, latCell: Corner Vertices of Dual Mesh cells
lonVertex, latVertex: Center Vertices of Dual Mesh cells
cellsOnVertex: Vertex indices that surround each Dual Mesh cell
cellsOnEdge: Vertex indices that saddle a given edge
Constructing a Grid Object#
Note
Since we only have a Grid file and no Data file in this example, we will be working exclusively with the Grid
class to investigate the grid topology and not with UxDataset
or UxDataArray
data structures.
The xarray.Dataset
that we opened in the previous section stores the coordinates and connectivity variables according to the MPAS specification standards for both the Primal and Dual meshes together in a single dataset. Here, instead of opening up the dataset using Xarray, we can pass through the path into our open_grid
method to construct an instance of a Grid
class. This Grid
can take in a use_dual
parameter to select whether to construct the Primal or Dual mesh, parsing and encoding the appropriate variables in the UGRID conventions.
import uxarray as ux
primal_mesh = ux.open_grid(mpas_dataset_filepath, use_dual=False)
dual_mesh = ux.open_grid(mpas_dataset_filepath, use_dual=True)
primal_mesh
/home/docs/checkouts/readthedocs.org/user_builds/uxarray/conda/v2023.11.0/lib/python3.11/site-packages/uxarray/grid/grid.py:141: FutureWarning: 'Mesh2' prefix used in dimension, coordinate, and connectivity attributes (i.e. Mesh2_face_nodes) will be dropped in a future release.
warn(
<uxarray.Grid>
Original Grid Type: MPAS
Grid Dimensions:
* nMesh2_node: 320
* nMesh2_face: 162
* nMesh2_edge: 480
* nMaxMesh2_face_nodes: 6
* nMaxNumFacesPerNode: 3
* Two: 2
Grid Coordinates (Latitude & Longitude):
* Mesh2_node_x: (320,)
* Mesh2_node_y: (320,)
* Mesh2_edge_x: (480,)
* Mesh2_edge_y: (480,)
* Mesh2_face_x: (162,)
* Mesh2_face_y: (162,)
Grid Coordinates (Cartesian):
* Mesh2_node_cart_x: (320,)
* Mesh2_node_cart_y: (320,)
* Mesh2_node_cart_z: (320,)
* Mesh2_edge_cart_x: (480,)
* Mesh2_edge_cart_y: (480,)
* Mesh2_edge_cart_z: (480,)
* Mesh2_face_cart_x: (162,)
* Mesh2_face_cart_y: (162,)
* Mesh2_face_cart_z: (162,)
Grid Connectivity Variables:
* Mesh2_node_faces: (320, 3)
* Mesh2_edge_nodes: (480, 2)
* Mesh2_face_nodes: (162, 6)
* Mesh2_face_edges: (162, 6)
* nNodes_per_face: (162,)
dual_mesh
<uxarray.Grid>
Original Grid Type: MPAS
Grid Dimensions:
* nMesh2_node: 162
* nMesh2_face: 320
* nMesh2_edge: 480
* nMaxMesh2_face_nodes: 3
* nMaxNumFacesPerNode: 6
* Two: 2
Grid Coordinates (Latitude & Longitude):
* Mesh2_node_x: (162,)
* Mesh2_node_y: (162,)
* Mesh2_edge_x: (480,)
* Mesh2_edge_y: (480,)
* Mesh2_face_x: (320,)
* Mesh2_face_y: (320,)
Grid Coordinates (Cartesian):
* Mesh2_node_cart_x: (162,)
* Mesh2_node_cart_y: (162,)
* Mesh2_node_cart_z: (162,)
* Mesh2_edge_cart_x: (480,)
* Mesh2_edge_cart_y: (480,)
* Mesh2_edge_cart_z: (480,)
* Mesh2_face_cart_x: (320,)
* Mesh2_face_cart_y: (320,)
* Mesh2_face_cart_z: (320,)
Grid Connectivity Variables:
* Mesh2_node_faces: (162, 6)
* Mesh2_edge_nodes: (480, 2)
* Mesh2_face_nodes: (320, 3)
* Mesh2_face_edges: (320, 3)
* nNodes_per_face: (320,)
Relationship between MPAS and UGRID#
In the previous two sections, we outlined the set of grid variables used to describe the Primal and Dual meshes and how to open an MPAS grid in UXarray. Here, we provide an overview of how we represent both meshes in the UGRID conventions and how the original grid variables were modified to meet these conventions.
Grid Variables (Primal Mesh)#
Mesh2_node_x
& Mesh2_node_y
Longitude and Latitude coordinates of the Primal Mesh corner nodes
Derived from
lonVertex
&latVertex
Converted from Radians to Degrees
Mesh2_face_x
& Mesh2_face_y
Longitude and Latitude coordinates of the Primal Mesh center nodes
Derived from
lonCell
&latCell
Converted from Radians to Degrees
Mesh2_face_nodes
Connectivity array describing which nodes make up a face
Derived from
verticesOnCell
Padding is replaced with
INT_FILL_VALUE
Missing Values (zeros) replaced with
INT_FILL_VALUE
Converted to zero-index
Mesh2_edge_nodes
Connectivity array describing which nodes link to form each edge
Derived from
verticesOnEdge
Padding is replaced with
INT_FILL_VALUE
Missing Values (zeros) replaced with
INT_FILL_VALUE
Converted to zero-index
Grid Variables (Dual Mesh)#
Mesh2_node_x
& Mesh2_node_y
Longitude and Latitude coordinates of the Dual Mesh vertices
Derived from
lonCell
&latCell
, the centers of the Primal MeshConverted from Radians to Degrees
Mesh2_face_x
& Mesh2_face_y
Longitude and Latitude coordinates of the Dual Mesh centers
Derived from
lonVertex
&latVertex
, the vertices of the Primal MeshConverted from Radians to Degrees
Mesh2_face_nodes
Connectivity array describing which nodes make up a face
Derived from
verticesOnCell
Padding is replaced with
INT_FILL_VALUE
Missing Values (zeros) replaced with
INT_FILL_VALUE
Converted to zero-index
Mesh2_edge_nodes
Connectivity array describing which nodes link to form each edge
Derived from
verticesOnEdge
Padding is replaced with
INT_FILL_VALUE
Missing Values (zeros) replaced with
INT_FILL_VALUE
Converted to zero-index
Functionality#
Face Area Calculation#
Using our parsed attributes, we can determine whether our unstructured grid lies on the surface of a sphere by accessing the on_a_sphere
attribute.
primal_mesh.parsed_attrs['on_a_sphere']
'YES'
Simiarly, we can access the sphere_radius
attribute.
primal_mesh.parsed_attrs['sphere_radius']
1.0
Since our mesh lies on a sphere, we would expect our total surface area to follow the equation
\(4{\pi}{r^2}\)
We can use the value of the sphere_radius
attribute to calculate the expected total surface area.
import numpy as np
sphere_r = primal_mesh.parsed_attrs['sphere_radius']
expected_area = 4 * np.pi * (sphere_r)**2
expected_area
12.566370614359172
UXarray
can be used to compute the face area of each face on our grid.
primal_mesh_face_areas = primal_mesh.face_areas
primal_mesh_face_areas
array([0.06733676, 0.06733676, 0.06733676, 0.06733676, 0.06733676,
0.06733676, 0.06733676, 0.06733676, 0.06733676, 0.06733676,
0.06733676, 0.06733676, 0.07879559, 0.07879561, 0.07879561,
0.07879559, 0.07879559, 0.07879559, 0.07879559, 0.07879559,
0.07879559, 0.07879561, 0.07879561, 0.07879559, 0.07879561,
0.07879561, 0.07879559, 0.07879559, 0.07879561, 0.07879559,
0.07879559, 0.07879561, 0.07879559, 0.07879559, 0.07879561,
0.07879559, 0.07879559, 0.07879559, 0.07879559, 0.07879559,
0.07879559, 0.07879561, 0.07631252, 0.07631253, 0.07631253,
0.07631252, 0.07631254, 0.07631254, 0.07631252, 0.07631253,
0.07631254, 0.07631252, 0.07631253, 0.07631254, 0.07631252,
0.07631254, 0.07631253, 0.07631254, 0.07631254, 0.07631253,
0.07631253, 0.07631252, 0.07631254, 0.07631253, 0.07631253,
0.07631252, 0.07631254, 0.07631253, 0.07631253, 0.07631252,
0.07631254, 0.07631253, 0.07631253, 0.07631253, 0.07631254,
0.07631253, 0.07631254, 0.07631253, 0.07631252, 0.07631252,
0.07631254, 0.07631253, 0.07631253, 0.07631252, 0.07631254,
0.07631254, 0.07631252, 0.07631254, 0.07631253, 0.07631253,
0.07631253, 0.07631253, 0.07631252, 0.07631254, 0.07631254,
0.07631252, 0.07631253, 0.07631252, 0.07631254, 0.07631254,
0.07631252, 0.07631254, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192, 0.08026192, 0.08026192, 0.08026192,
0.08026192, 0.08026192])
The total face (surface) area can be computed by summing over each value.
primal_mesh_face_areas.sum()
12.5663761450242
We can then compute the absolute error of our calculation.
abs(expected_area - primal_mesh_face_areas.sum())
5.530665028175008e-06
The same can be done for the Dual Mesh.
dual_mesh_face_areas = dual_mesh.face_areas
dual_mesh_face_areas.sum()
12.566372994960696
abs(expected_area - dual_mesh_face_areas.sum())
2.3806015239102862e-06
We can see that the total face area of both the Primal and Dual meshes is within 1e-6 of the expected area. For a more detailed explanation of the face area calculation and ways to obtain more precision, check out our other notebooks.
Visualization#
To visually confirm that the Primal and Dual meshes have been correctly translated to the UGRID conventions, the following showcases a basic visualization of the mesh structure.
To display the mesh, the vertices are used to construct a PolyCollection
object, which represents a group of polygons on a plane. The collection is then fed into Matplotlib to be visualized.
Note
The following visualizations showcase the mesh structure excluding any faces that are located on the boundary between 0 and 360 degrees. There is work currently being done to handle these cases and to provide a more sophisticated visualization API in UXarray. Stay tuned!
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection
from uxarray.grid.connectivity import close_face_nodes
Helper Functions#
In order to make our data compatible wth the PolyCollection
class, we need to first convert it into an appropriate representation. The following helper function takes in a Grid
object (either Primal or Dual mesh) and returns the vertices of each face that can be interpreted by the PolyCollection
.
def get_verts(grid):
nodes_per_face = grid.nNodes_per_face.values
x = grid.Mesh2_node_x.values
y = grid.Mesh2_node_y.values
closed_face_nodes = close_face_nodes(grid.Mesh2_face_nodes.values,
grid.nMesh2_face,
grid.nMaxMesh2_face_nodes)
verts = []
for face_node, n_nodes in zip(closed_face_nodes, nodes_per_face):
polygon_x = x[face_node[0:n_nodes]]
polygon_y = y[face_node[0:n_nodes]]
# exclude faces near the boundary
if np.any(np.abs(polygon_x - 180) > 170):
continue
vert = np.array([polygon_x, polygon_y])
verts.append(vert.T)
return verts
Primal Mesh#
The following displays the Primal Mesh, which is composed mostly of hexagons representing the Voronoi regions,
primal_vertices = get_verts(primal_mesh)
primal_polygons = PolyCollection(primal_vertices, edgecolors='Black', alpha=0.5)
fig, ax = plt.subplots()
ax.add_collection(primal_polygons, autolim=False)
ax.set_xlim(0, 360)
ax.set_ylim(-80, 80)
ax.set_title("Primal Mesh")
Text(0.5, 1.0, 'Primal Mesh')
Dual Mesh#
The follwing displays the Dual Mesh, which is composed entirely of Delaunay triangles.
dual_vertices = get_verts(dual_mesh)
dual_polygons = PolyCollection(dual_vertices, edgecolors='Black', alpha=0.5)
fig, ax = plt.subplots()
ax.add_collection(dual_polygons, autolim=False)
ax.set_xlim(0, 360)
ax.set_ylim(-80, 80)
ax.set_title("Dual Mesh")
Text(0.5, 1.0, 'Dual Mesh')
Primal and Dual Meshes Overlayed#
As mentioned earlier, the Primal and Dual meshes are related to one another. The vertices of the Primal mesh are each at the center of a Delaunay triangle in the Dual mesh. Likewise, the vertices of the Dual mesh are each at the center of a Voronoi region in the Primal mesh. This relationship can be seen in the visualization below
dual_vertices = get_verts(dual_mesh)
dual_polygons = PolyCollection(dual_vertices, edgecolors='Black', alpha=0.5)
primal_vertices = get_verts(primal_mesh)
primal_polygons = PolyCollection(primal_vertices, edgecolors='Black', alpha=0.5)
fig, ax = plt.subplots()
ax.add_collection(primal_polygons, autolim=False)
ax.add_collection(dual_polygons, autolim=False)
ax.set_xlim(0, 360)
ax.set_ylim(-80, 80)
ax.set_title("Dual & Primal Mesh")
Text(0.5, 1.0, 'Dual & Primal Mesh')
References#
[1] https://mpas-dev.github.io/
[2] https://mpas-dev.github.io/files/documents/MPAS-MeshSpec.pdf