DRAFT: area correction for when an edge is the line of constant lattitude

docs/user-guide/area_calc.ipynb

@@ -22,7 +22,8 @@
     "2. Options for `Grid.calculate_total_face_area` Function\n",
     "3. Getting Area of Individual Faces\n",
     "4. Calculate Area of a Single Triangle in Cartesian Coordinates\n",
-    "5. Calculate Area from Multiple Faces in Spherical Coordinates"
+    "5. Calculate Area from Multiple Faces in Spherical Coordinates\n",
+    "6. Demonstrate Area Correction at Line of Constant Lattitude"
    ]
   },
   {
@@ -370,6 +371,164 @@
     "area, jacobian = verts_grid.compute_face_areas()\n",
     "area"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 6. Area Correction\n",
+    "\n",
+    "The correction, \\(A\\), is calculated using the following formula:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "```{math}\n",
+    "A = 2arctan(z(x_{1}y_{2} - x_{2} y_{1}) / {x_{1}^2 + y_{1}^2 + x_{1} x_{2} + y_{1} y_{2}) - z (x_{1} y_{2} - x_{2} y_{1} / x_{1} x_{2} + y_{1} y_{2}\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Where:\n",
+    "<ul>\n",
+    "\n",
+    "<li>(x<sub>1</sub>, y<sub>1</sub>, z) are the Cartesian coordinates of the first node.</li>\n",
+    "\n",
+    "<li>(x<sub>2</sub>, y<sub>2</sub>, z) are the Cartesian coordinates of the second node (note that the z coordinate is the same for both nodes).</li>\n",
+    "\n",
+    "</ul>\n",
+    "\n",
+    "### Assumptions:\n",
+    "- This formula assumes that the input coordinates \\((x_1, y_1)\\) and \\((x_2, y_2)\\) are normalized (i.e., they lie on the unit sphere).\n",
+    "\n",
+    "\n",
+    "For the same large triangle used in **Section 4**, we calculate the area correction term when an edge lies along the line of constant latitude.\n",
+    "\n",
+    "### The following code:\n",
+    "- Plots the triangle.\n",
+    "- Marks the edges with different colors.\n",
+    "- Highlights the edge along the line of constant latitude.\n",
+    "- Marks the coordinates of the vertices."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline \n",
+    "import matplotlib.pyplot as plt\n",
+    "import cartopy.crs as ccrs\n",
+    "import cartopy.feature as cfeature\n",
+    "\n",
+    "# Convert the points to latitude and longitude\n",
+    "points = [\n",
+    "    (point[0] * 180 / 3.14159, point[1] * 180 / 3.14159)\n",
+    "    for vert in verts\n",
+    "    for point in vert\n",
+    "]\n",
+    "\n",
+    "# Node names\n",
+    "node_names = [\"Node 1\", \"Node 2\", \"Node 3\"]\n",
+    "\n",
+    "# Edge names and colors\n",
+    "edge_names = [\"Edge 1\", \"Edge 2\", \"Edge 3\"]\n",
+    "edge_colors = [\"red\", \"green\", \"blue\"]\n",
+    "\n",
+    "# Create a new figure\n",
+    "fig = plt.figure(figsize=(10, 5))\n",
+    "\n",
+    "# Create a GeoAxes in the Orthographic projection\n",
+    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(central_longitude=0))\n",
+    "\n",
+    "# Add coastlines and gridlines for reference\n",
+    "ax.add_feature(cfeature.COASTLINE)\n",
+    "ax.gridlines(draw_labels=True)\n",
+    "\n",
+    "# Plot the points, edges, and labels\n",
+    "for i in range(len(points)):\n",
+    "    lon1, lat1 = points[i]\n",
+    "    lon2, lat2 = points[(i + 1) % len(points)]\n",
+    "\n",
+    "    # Plot the edge\n",
+    "    ax.plot(\n",
+    "        [lon1, lon2],\n",
+    "        [lat1, lat2],\n",
+    "        color=edge_colors[i],\n",
+    "        transform=ccrs.Geodetic(),\n",
+    "        label=edge_names[i],\n",
+    "    )\n",
+    "\n",
+    "    # Add edge label (adjust the offset as needed)\n",
+    "    ax.text(\n",
+    "        (lon1 + lon2) / 2,\n",
+    "        (lat1 + lat2) / 2,\n",
+    "        edge_names[i],\n",
+    "        transform=ccrs.Geodetic(),\n",
+    "        color=edge_colors[i],\n",
+    "    )\n",
+    "\n",
+    "    # Add node label with Cartesian coordinates\n",
+    "    cartesian_coords = verts[0][i]  # Get Cartesian coordinates from verts\n",
+    "    label = f\"{node_names[i]}\\n[{cartesian_coords[0]:.2f}, {cartesian_coords[1]:.2f}, {cartesian_coords[2]:.2f}]\"\n",
+    "    ax.text(lon1, lat1, label, transform=ccrs.Geodetic(), fontsize=8, color=\"black\")\n",
+    "\n",
+    "# Show the full globe\n",
+    "ax.set_global()\n",
+    "\n",
+    "# Add a legend for the edges\n",
+    "plt.legend()\n",
+    "\n",
+    "# Display the plot\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The default for correcttion is False, so we set it to True to calculate the correction term."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Area of the triangle with correction\n",
+    "area, jacobian = vgrid.compute_face_areas(quadrature_rule=\"gaussian\", order=5)\n",
+    "area"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "- **Edge 3** (blue) and **Edge 2** (green) of the triangle lie along lines of constant latitude.\n",
+    "- Without correction, the area is calculated as `1.047`. With the correction, the area becomes `1.047 + 2 * 0.140`, which gives the correct area of the triangle.\n",
+    "\n",
+    "\n",
+    "- **Edge 3** (red) also lies along a line of constant latitude but passes through the poles. Therefore, it is not considered for the correction term.\n",
+    "\n",
+    "To compute the area with the correction term applied, pass the argument ``correct_area=True`` to the function call. This enables automatic adjustment for edges along lines of constant latitude."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Continuing with our example from above and now we will correct the area\n",
+    "area = vgrid.compute_face_areas(quadrature_rule=\"gaussian\", order=5, correct_area=True)"
+   ]
   }
  ],
  "metadata": {
@@ -388,7 +547,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.9"
+   "version": "3.11.11"
   }
  },
  "nbformat": 4,

test/constants.py

@@ -7,6 +7,7 @@
 NNODES_outRLL1deg = 64442
 DATAVARS_outCSne30 = 4
 TRI_AREA = 1.047
+CORRECTED_TRI_AREA = 1.3281
 # 4*Pi is 12.56
 MESH30_AREA = 12.566
 PSI_INTG = 12.566

test/test_grid.py

@@ -266,15 +266,16 @@ def test_face_areas_calculate_total_face_area_triangle():
     # 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, correct_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."""
     area = ux.open_grid(gridfile_CSne30).calculate_total_face_area()