# 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'.

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)
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'.

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)
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'.

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

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(
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
```