Face Area Calculations#

Authors: Rajeev Jain

Overview#

This notebook will showcase the different area calculation options provided by uxarray

For more details on how to load in data, check out our previous usage example

This notebook has the following sections:

  1. Calculate Total Face Area

  2. Options for Grid.calculate_total_face_area Function

  3. Getting Area of Individual Faces

  4. Calculate Area of a Single Triangle in Cartesian Coordinates

  5. Calculate Area from Multiple Faces in Spherical Coordinates

  6. Area Calculation without Grid Object

We will be using the outCSne30.ug grid file, which is encoded in the UGRID convention.

import uxarray as ux
import numpy as np
import xarray as xr

# Base data path
base_path = "../../test/meshfiles/"

# Path to Grid files
ugrid_path = base_path + "/ugrid/outCSne30/outCSne30.ug"

# Load grid files and create UXarray Grid objects
ugrid_ds = xr.open_dataset(ugrid_path)

ugrid = ux.open_grid(ugrid_ds)
ugrid

<uxarray.Grid>
Original Grid Type: UGRID
Grid Dimensions:
  * nMesh2_face: 5400
  * nMaxMesh2_face_nodes: 4
  * nMesh2_node: 5402
Grid Coordinate Variables:
  * Mesh2_node_x: (5402,)
  * Mesh2_node_y: (5402,)
Grid Connectivity Variables:
  * Mesh2_face_nodes: (5400, 4)
  * nNodes_per_face: (5400,)

1. Calculate Total Face Area#

We can calculate the total face area by calling the function Grid.calculate_total_face_area(). Since our dataset lies on the unit sphere, our expected area is 4*pi, which is approximately 12.56

t4_area = ugrid.calculate_total_face_area()
t4_area

12.566370614678554

2. Options for Grid.calculate_total_face_area Function#

By default, Grid.calculate_total_face_area uses a Triangular Quadrature Rule and an Order of 4. However, you can specify the Quadrature Rule and Order as follows:

Order:

   1 to 10              for gaussian

   1, 4, 8, 10 and 12   for triangular

t1_area = ugrid.calculate_total_face_area(quadrature_rule="triangular", order=1)
t1_area

12.571403993719983

For the result above, notice that the area is slightly different than the first calculation we made.

Now we use Triangular Quadrature Rule and Order 1

Using a lower order is faster, but at the sacrifice of accuracy.

Generally, gaussian quadrature rule is more accurate than the triangular quadrature rule. Additionally, a higher order comes at the cost of computation time, but produces a more accurate result. See uxarray/get_quadratureDG.py file for details on quadrature points and weights.

3. Getting Area of Individual Faces#

We can access the Grid attribute Grid.face_area to access the individual face areas. If you have not run a face area calculation yet, it will run the Grid.compute_face_areas and cache the value.

ugrid.face_areas

array([0.00211238, 0.00211285, 0.00210788, ..., 0.00210788, 0.00211285,
       0.00211238])

Now calculate the area again with the Grid.compute_face_areas function using arguments: quadrature_rule “gaussian” and order 4

all_face_areas = ugrid.compute_face_areas(quadrature_rule="gaussian", order=4)
g4_area = all_face_areas.sum()
g4_area

12.566370614359112

Now we compare the values with actual know value and report error for each of the three cases above.

actual_area = 4 * np.pi
diff_t4_area = np.abs(t4_area - actual_area)
diff_t1_area = np.abs(t1_area - actual_area)
diff_g4_area = np.abs(g4_area - actual_area)

diff_t1_area, diff_t4_area, diff_g4_area

(0.005033379360810386, 3.1938185429680743e-10, 6.039613253960852e-14)

As we can see, it is clear that the Gaussian Quadrature Rule with Order 4 is the most accurate, and the Triangular Quadrature Rule with Order 1 is the least accurate.

4. Calculate Area of a Single Triangle in Cartesian Coordinates#

For this section, we create a single triangle face with 3 vertices. By default, in uxarray, we assume that the coordinate system is spherical (lat / lon), however if you want to use cartesian coordinates, you must pass through islatlon = False into the Grid constructor.

Assume the units in meters - this is a big triangle!

verts = [[[0.57735027, -5.77350269e-01, -0.57735027],
         [0.57735027, 5.77350269e-01, -0.57735027],
         [-0.57735027, 5.77350269e-01, -0.57735027]]]

# load our vertices into a UXarray Grid object
vgrid = ux.open_grid(verts, islatlon=False)

vgrid

<uxarray.Grid>
Original Grid Type: Face Vertices
Grid Dimensions:
  * nMesh2_node: 3
  * nMesh2_face: 1
  * nMaxMesh2_face_nodes: 3
Grid Coordinate Variables:
  * Mesh2_node_cart_x: (3,)
  * Mesh2_node_cart_y: (3,)
  * Mesh2_node_cart_z: (3,)
Grid Connectivity Variables:
  * Mesh2_face_nodes: (1, 3)
  * nNodes_per_face: (1,)

vgrid.calculate_total_face_area()

/home/docs/checkouts/readthedocs.org/user_builds/uxarray/conda/latest/lib/python3.11/site-packages/uxarray/grid/coordinates.py:63: NumbaPendingDeprecationWarning: 
Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'node' of function 'normalize_in_place'.

For more information visit https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types

File "../../../../conda/latest/lib/python3.11/site-packages/uxarray/grid/coordinates.py", line 85:
@njit(cache=ENABLE_JIT_CACHE)
def normalize_in_place(node):
^

  [dx, dy, dz] = normalize_in_place(node_coord)
/home/docs/checkouts/readthedocs.org/user_builds/uxarray/conda/latest/lib/python3.11/site-packages/uxarray/grid/coordinates.py:63: NumbaPendingDeprecationWarning: 
Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'object' of function 'impl_np_array.<locals>.impl'.

For more information visit https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types

File "../../../../conda/latest/lib/python3.11/site-packages/numba/np/arrayobj.py", line 5252:

    def impl(object, dtype=None):
    ^

  [dx, dy, dz] = normalize_in_place(node_coord)
/home/docs/checkouts/readthedocs.org/user_builds/uxarray/conda/latest/lib/python3.11/site-packages/numba/core/ir_utils.py:2149: NumbaPendingDeprecationWarning: 
Encountered the use of a type that is scheduled for deprecation: type 'reflected list' found for argument 'node_coord' of function 'node_xyz_to_lonlat_rad'.

For more information visit https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-reflection-for-list-and-set-types

File "../../../../conda/latest/lib/python3.11/site-packages/uxarray/grid/coordinates.py", line 41:
@njit(cache=ENABLE_JIT_CACHE)
def node_xyz_to_lonlat_rad(node_coord):
^

  warnings.warn(NumbaPendingDeprecationWarning(msg, loc=loc))

1.071941993854821

Additionally, if you are using a unit other than meters, you can update the units as follows

# # Set correct units for the x and y coordinates
# vgrid.Mesh2_node_x.attrs["units"] = "km"
# vgrid.Mesh2_node_y.attrs["units"] = "km"
# vgrid.Mesh2_node_z.attrs["units"] = "km"  # This is just a placeholder, UXarray does not support 3D meshes

# Calculate the area of the triangle
area_gaussian = vgrid.calculate_total_face_area(
            quadrature_rule="gaussian", order=5)
area_gaussian

1.0475702709086987

5. Calculate Area from Multiple Faces in Spherical Coordinates#

Similar to above, we can construct a Grid object with multiple faces by passing through a set of vertices. Here we define 3 six-sided faces.

faces_verts_ndarray = np.array([
    np.array([[150, 10, 0], [160, 20, 0], [150, 30, 0], [135, 30, 0], [125, 20, 0],
              [135, 10, 0]]),
    np.array([[125, 20, 0], [135, 30, 0], [125, 60, 0], [110, 60, 0], [100, 30, 0],
              [105, 20, 0]]),
    np.array([[95, 10, 0], [105, 20, 0], [100, 30, 0], [85, 30, 0], [75, 20, 0],
              [85, 10, 0]]),
])

We want our units to be spherical, so we pass through islatlon=True. Additionally, if islatlon is not passed through, it will default to spherical coordinates.

verts_grid = ux.open_grid(faces_verts_ndarray, islatlon=False)

verts_grid

<uxarray.Grid>
Original Grid Type: Face Vertices
Grid Dimensions:
  * nMesh2_node: 14
  * nMesh2_face: 3
  * nMaxMesh2_face_nodes: 6
Grid Coordinate Variables:
  * Mesh2_node_cart_x: (14,)
  * Mesh2_node_cart_y: (14,)
  * Mesh2_node_cart_z: (14,)
Grid Connectivity Variables:
  * Mesh2_face_nodes: (3, 6)
  * nNodes_per_face: (3,)

verts_grid.compute_face_areas(latlon=False)

array([1.52225680e-19, 1.77829591e-18, 6.46116331e-19])

6. Area Calculation without Grid Object#

If you want to compute the face area without using the Grid object, many of the functions for computing the face are can be accessed through uxarray.helpers

from uxarray.grid.area import calculate_face_area

Grid.calculate_face_area takes in three coordinate variables (x, y, z) in the form of numpy arrays and the coordinate type (either spherical or artesian) and computes the face area from the set of points

cart_x = np.array([
   0.577340924821405, 0.577340924821405, 0.577340924821405,
   0.577340924821405
])
cart_y = np.array([
   0.577343045516932, 0.577343045516932, -0.577343045516932,
   -0.577343045516932
])
cart_z = np.array([
   0.577366836872017, -0.577366836872017, 0.577366836872017,
   -0.577366836872017
])

calculate_face_area(cart_x, cart_y, cart_z, coords_type="cartesian")

2.0901881354848064