Area correction for when an edge is the line of constant latitude (#1120)

* o First draft area correction

* o Add test and notebook to showcase area correction

* o Add more to notebook, mathematical formulation of correction and plot of original triangle

* o fix notebook

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* o Cleanup notebook

* o Fix spaces

* Use $ for formula in notebook, mathjax seems to be configured correctly

* Attempted fix for notebook error

* notebook fix

* Changed markdown style of notebook

* Update area_calc.ipynb

* Update area_calc.ipynb

* o Remove bad file

* o checked in files my mistake

* o Remove unnecessary files

* o Complete TODO fix area correct sign based on edge lon and hemisphere. Add notebook example. TODO: Add test and corner cases

* o Fix tests

* o change grid coords to match new constants

* o Correct math formula display for doc page html, use latex in notebook

* Fix tests and notebook, add plot, cleanup code, getting ready to get out of draft mode

* Remove all o/p from notebook

* o Delete bad files

* o Use new mathjax

* o Fix tests and notebook, also use MathJAX3 and 2 is no longer supported

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* Update uxarray/grid/area.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* o Add docstring for correct_area

* o notebook clear output

* o Add sep for _fesom, removes warnings, remove one more comment

* o ruff fix

* o Use new simple formatting for doc math formulae

* o Add docstring for correct_area in grid.py

* o Add another example for face area calculation that motivates the usage of correction for zonal averaging calculation etc.

* o Remove hemisphere string and add check for skipping when lon_diff of edge is more than 180 deg

* o extra normalization and simplify calc

* o Use pythonic variable names for Gaussian and triangular quadrature points for area calc.

* o Use latitude_adjusted_area as a flag to indicate area correction requirement when users wants to specify  that all edges in a mesh that have the same latitude endpoints are lines of constant latitude, Add TODO for const lattitude edge flag and spherical changes

---------

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>
Co-authored-by: Aaron Zedwick <aazedwick@gmail.com>
Co-authored-by: Aaron Zedwick <95507181+aaronzedwick@users.noreply.github.com>

Author: 3 people

docs/user-guide/area_calc.ipynb

@@ -6,7 +6,8 @@
 NNODES_outCSne30 = 5402
 NNODES_outRLL1deg = 64442
 DATAVARS_outCSne30 = 4
-TRI_AREA = 1.047
+TRI_AREA = 0.02216612469199045
+CORRECTED_TRI_AREA = 0.02244844510268421
 # 4*Pi is 12.56
 MESH30_AREA = 12.566
 PSI_INTG = 12.566

test/constants.py

@@ -265,24 +265,27 @@ def test_operators_ne():
 
 def test_face_areas_calculate_total_face_area_triangle():
     """Create a uxarray grid from vertices and saves an exodus file."""
-    verts = [[[0.57735027, -5.77350269e-01, -0.57735027],
-              [0.57735027, 5.77350269e-01, -0.57735027],
-              [-0.57735027, 5.77350269e-01, -0.57735027]]]
+    verts = [
+    [[0.02974582, -0.74469018, 0.66674712],
+    [0.1534193, -0.88744577, 0.43462917],
+    [0.18363692, -0.72230586, 0.66674712]]
+    ]
 
     grid_verts = ux.open_grid(verts, latlon=False)
 
     # validate the grid
     assert grid_verts.validate()
 
-    # calculate area
-    area_gaussian = grid_verts.calculate_total_face_area(
-        quadrature_rule="gaussian", order=5)
-    nt.assert_almost_equal(area_gaussian, constants.TRI_AREA, decimal=3)
-
+    # calculate area without correction
     area_triangular = grid_verts.calculate_total_face_area(
         quadrature_rule="triangular", order=4)
     nt.assert_almost_equal(area_triangular, constants.TRI_AREA, decimal=1)
 
+    # calculate area
+    area_gaussian = grid_verts.calculate_total_face_area(
+        quadrature_rule="gaussian", order=5, latitude_adjusted_area=True)
+    nt.assert_almost_equal(area_gaussian, constants.CORRECTED_TRI_AREA, decimal=3)
+
 
 def test_face_areas_calculate_total_face_area_file():
     """Create a uxarray grid from vertices and saves an exodus file."""

test/test_grid.py

@@ -32,41 +32,62 @@
 def test_face_area_coords():
     """Test function for helper function get_all_face_area_from_coords."""
     # Note: currently only testing one face, but this can be used to get area of multiple faces
-    x = np.array([0.57735027, 0.57735027, -0.57735027])
-    y = np.array([-5.77350269e-01, 5.77350269e-01, 5.77350269e-01])
-    z = np.array([-0.57735027, -0.57735027, -0.57735027])
-
-    face_nodes = np.array([[0, 1, 2]]).astype(INT_DTYPE)
+    # Cartesian coordinates (x, y, z) for each city
+    # Index 0: Chicago, Index 1: Miami, Index 2: Newburgh, New York, USA.
+    x = np.array([0.02974582, 0.1534193, 0.18363692])
+    y = np.array([-0.74469018, -0.88744577, -0.72230586])
+    z = np.array([0.66674712, 0.43462917, 0.66674712])
+    face_nodes = np.array([[0, 1, 2]])
     face_dimension = np.array([3], dtype=INT_DTYPE)
 
-    area, jacobian = ux.grid.area.get_all_face_area_from_coords(
+    area, _ = ux.grid.area.get_all_face_area_from_coords(
         x, y, z, face_nodes, face_dimension, 3, coords_type="cartesian")
+    nt.assert_almost_equal(area, constants.TRI_AREA, decimal=5)
 
-    nt.assert_almost_equal(area, constants.TRI_AREA, decimal=1)
 
 def test_calculate_face_area():
     """Test function for helper function calculate_face_area - only one face."""
     # Note: currently only testing one face, but this can be used to get area of multiple faces
     # Also note, this does not need face_nodes, assumes nodes are in counterclockwise orientation
-    x = np.array([0.57735027, 0.57735027, -0.57735027])
-    y = np.array([-5.77350269e-01, 5.77350269e-01, 5.77350269e-01])
-    z = np.array([-0.57735027, -0.57735027, -0.57735027])
+    x = np.array([0.02974582, 0.1534193, 0.18363692])
+    y = np.array([-0.74469018, -0.88744577, -0.72230586])
+    z = np.array([0.66674712, 0.43462917, 0.66674712])
+
+    area, _ = ux.grid.area.calculate_face_area(
+        x, y, z, "gaussian", 5, "cartesian", latitude_adjusted_area=False)
+
+    nt.assert_almost_equal(area, constants.TRI_AREA, decimal=5)
+
+    area_corrected, _ = ux.grid.area.calculate_face_area(
+        x, y, z, "gaussian", 5, "cartesian", latitude_adjusted_area=True)
 
-    area, jacobian = ux.grid.area.calculate_face_area(
-        x, y, z, "gaussian", 5, "cartesian")
+    nt.assert_almost_equal(area_corrected, constants.CORRECTED_TRI_AREA, decimal=5)
+
+    # Make the same grid using lon/lat check area = constants.TRI_AREA
+    lon = np.array([-87.7126, -80.1918, -75.7355])
+    lat = np.array([41.8165, 25.7617, 41.8165])
+    face_nodes = np.array([[0, 1, 2]])
+
+    grid = ux.Grid.from_topology(
+        node_lon=lon,
+        node_lat=lat,
+        face_node_connectivity=face_nodes,
+        fill_value=-1,
+    )
 
-    nt.assert_almost_equal(area, constants.TRI_AREA, decimal=3)
+    area, _ = grid.compute_face_areas()
+    nt.assert_almost_equal(area, constants.TRI_AREA, decimal=5)
 
 def test_quadrature():
     order = 1
-    dG, dW = ux.grid.area.get_tri_quadratureDG(order)
+    dG, dW = ux.grid.area.get_tri_quadrature_dg(order)
     G = np.array([[0.33333333, 0.33333333, 0.33333333]])
     W = np.array([1.0])
 
     np.testing.assert_array_almost_equal(G, dG)
     np.testing.assert_array_almost_equal(W, dW)
 
-    dG, dW = ux.grid.area.get_gauss_quadratureDG(order)
+    dG, dW = ux.grid.area.get_gauss_quadrature_dg(order)
 
     G = np.array([[0.5]])
     W = np.array([1.0])