diff --git a/docs/user-guide/area_calc.ipynb b/docs/user-guide/area_calc.ipynb index 49e4d1ace..c9d9b1d86 100644 --- a/docs/user-guide/area_calc.ipynb +++ b/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", + "\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, diff --git a/test/constants.py b/test/constants.py index 4a54e1b1a..eeb23ecf3 100644 --- a/test/constants.py +++ b/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 diff --git a/test/test_grid.py b/test/test_grid.py index 4c8d98c76..8300b7590 100644 --- a/test/test_grid.py +++ b/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() diff --git a/uxarray/grid/area.py b/uxarray/grid/area.py index b785f87f8..ccbb18c36 100644 --- a/uxarray/grid/area.py +++ b/uxarray/grid/area.py @@ -7,7 +7,13 @@ @njit(cache=True) def calculate_face_area( - x, y, z, quadrature_rule="gaussian", order=4, coords_type="spherical" + x, + y, + z, + quadrature_rule="gaussian", + order=4, + coords_type="spherical", + correct_area=False, ): """Calculate area of a face on sphere. @@ -56,6 +62,22 @@ def calculate_face_area( # num triangles is two less than the total number of nodes num_triangles = num_nodes - 2 + if coords_type == "spherical": + # Preallocate arrays for Cartesian coordinates + n_points = len(x) + x_cartesian = np.empty(n_points) + y_cartesian = np.empty(n_points) + z_cartesian = np.empty(n_points) + + # Convert all points to Cartesian coordinates using an explicit loop + for i in range(n_points): + lon_rad = np.deg2rad(x[i]) + lat_rad = np.deg2rad(y[i]) + cartesian = _lonlat_rad_to_xyz(lon_rad, lat_rad) + x_cartesian[i], y_cartesian[i], z_cartesian[i] = cartesian + + x, y, z = x_cartesian, y_cartesian, z_cartesian + # Using tempestremap GridElements: https://github.com/ClimateGlobalChange/tempestremap/blob/master/src/GridElements.cpp # loop through all sub-triangles of face for j in range(0, num_triangles): @@ -63,16 +85,6 @@ def calculate_face_area( node2 = np.array([x[j + 1], y[j + 1], z[j + 1]], dtype=x.dtype) node3 = np.array([x[j + 2], y[j + 2], z[j + 2]], dtype=x.dtype) - if coords_type == "spherical": - node1 = _lonlat_rad_to_xyz(np.deg2rad(x[0]), np.deg2rad(y[0])) - node1 = np.asarray(node1) - - node2 = _lonlat_rad_to_xyz(np.deg2rad(x[j + 1]), np.deg2rad(y[j + 1])) - node2 = np.asarray(node2) - - node3 = _lonlat_rad_to_xyz(np.deg2rad(x[j + 2]), np.deg2rad(y[j + 2])) - node3 = np.asarray(node3) - for p in range(len(dW)): if quadrature_rule == "gaussian": for q in range(len(dW)): @@ -92,9 +104,85 @@ def calculate_face_area( area += dW[p] * jacobian jacobian += jacobian + # check if the any edge is on the line of constant latitude + # which means we need to check edges for same z-coordinates and call area correction routine + correction = 0.0 + if correct_area: + for i in range(num_nodes): + node1 = np.array([x[i], y[i], z[i]], dtype=x.dtype) + node2 = np.array( + [ + x[(i + 1) % num_nodes], + y[(i + 1) % num_nodes], + z[(i + 1) % num_nodes], + ], + dtype=x.dtype, + ) + # Check if z-coordinates are approximately equal + if np.isclose(node1[2], node2[2]): + # Check if the edge passes through a pole + passes_through_pole = edge_passes_through_pole(node1, node2) + print("Check if edge passes through pole: ", passes_through_pole) + + if passes_through_pole: + # Skip the edge if it passes through a pole + continue + else: + # Calculate the correction term + correction = area_correction(node1, node2) + print( + "For Node 1 ", + node1, + "\n and Node 2", + node2, + "\nCORRECTION", + correction, + ) + correction += correction + + if correction != 0.0: + print("AREA Before Correction", area) + + # TODO: Fix sign of the calculated correction? + area += correction + + if correction != 0.0: + print("AREA After Correction", area) + return area, jacobian +@njit(cache=True) +def edge_passes_through_pole(node1, node2): + """ + Check if the edge passes through a pole. + + Parameters: + - node1: first node of the edge (normalized). + - node2: second node of the edge (normalized). + + Returns: + - bool: True if the edge passes through a pole, False otherwise. + """ + # Calculate the normal vector to the plane defined by the origin, node1, and node2 + n = np.cross(node1, node2) + + # Check for numerical stability issues with the normal vector + if np.allclose(n, 0): + # Handle cases where the cross product is near zero, such as when nodes are nearly identical or opposite + return False + + # Normalize the normal vector + n = n / np.linalg.norm(n) + + # North and South Pole vectors + p_north = np.array([0.0, 0.0, 1.0]) + p_south = np.array([0.0, 0.0, -1.0]) + + # Check if the normal vector is orthogonal to either pole + return np.isclose(np.dot(n, p_north), 0) or np.isclose(np.dot(n, p_south), 0) + + @njit(cache=True) def get_all_face_area_from_coords( x, @@ -106,6 +194,7 @@ def get_all_face_area_from_coords( quadrature_rule="triangular", order=4, coords_type="spherical", + correct_area=False, ): """Given coords, connectivity and other area calculation params, this routine loop over all faces and return an numpy array with areas of each @@ -141,6 +230,10 @@ def get_all_face_area_from_coords( ------- area of all faces : ndarray """ + # this casting helps to prevent the type mismatch + x = np.asarray(x, dtype=np.float64) + y = np.asarray(y, dtype=np.float64) + z = np.asarray(z, dtype=np.float64) n_face, n_max_face_nodes = face_nodes.shape @@ -161,7 +254,7 @@ def get_all_face_area_from_coords( # After getting all the nodes of a face assembled call the cal. face area routine face_area, face_jacobian = calculate_face_area( - face_x, face_y, face_z, quadrature_rule, order, coords_type + face_x, face_y, face_z, quadrature_rule, order, coords_type, correct_area ) # store current face area area[face_idx] = face_area @@ -170,6 +263,55 @@ def get_all_face_area_from_coords( return area, jacobian +@njit(cache=True) +def area_correction(node1, node2): + """ + Calculate the area correction A using the given formula. + + Parameters: + - node1: first node of the edge (normalized). + - node2: second node of the edge (normalized). + - z: z-coordinate (shared by both points and part of the formula, normalized). + + Returns: + - A: correction term of the area, when one of the edges is a line of constant latitude + """ + x1 = node1[0] + y1 = node1[1] + x2 = node2[0] + y2 = node2[1] + z = node1[2] + + # Calculate terms + term1 = x1 * y2 - x2 * y1 + den1 = x1**2 + y1**2 + x1 * x2 + y1 * y2 + den2 = x1 * x2 + y1 * y2 + + # Helper function to handle arctan quadrants + def arctan_quad(y, x): + if x > 0: + return np.arctan(y / x) + elif x < 0 and y >= 0: + return np.arctan(y / x) + np.pi + elif x < 0 and y < 0: + return np.arctan(y / x) - np.pi + elif x == 0 and y > 0: + return np.pi / 2 + elif x == 0 and y < 0: + return -np.pi / 2 + else: + return 0 # x == 0 and y == 0 case + + # Compute angles using arctan + angle1 = arctan_quad(z * term1, den1) + angle2 = arctan_quad(term1, den2) + + # Compute A + A = abs(2 * angle1 - z * angle2) + print(x1, y1, x2, y2, z, "correction:", A) + return A + + @njit(cache=True) def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB): """Calculate Jacobian of a spherical triangle. This is a helper function diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index aa93192d2..1b11d7d05 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -1766,7 +1766,10 @@ def encode_as(self, grid_type: str) -> xr.Dataset: return out_ds def calculate_total_face_area( - self, quadrature_rule: Optional[str] = "triangular", order: Optional[int] = 4 + self, + quadrature_rule: Optional[str] = "triangular", + order: Optional[int] = 4, + correct_area: Optional[bool] = False, ) -> float: """Function to calculate the total surface area of all the faces in a mesh. @@ -1784,7 +1787,9 @@ def calculate_total_face_area( """ # call function to get area of all the faces as a np array - face_areas, face_jacobian = self.compute_face_areas(quadrature_rule, order) + face_areas, face_jacobian = self.compute_face_areas( + quadrature_rule, order, correct_area=correct_area + ) return np.sum(face_areas) @@ -1793,6 +1798,7 @@ def compute_face_areas( quadrature_rule: Optional[str] = "triangular", order: Optional[int] = 4, latlon: Optional[bool] = True, + correct_area: Optional[bool] = False, ): """Face areas calculation function for grid class, calculates area of all faces in the grid. @@ -1860,6 +1866,7 @@ def compute_face_areas( quadrature_rule, order, coords_type, + correct_area, ) min_jacobian = np.min(self._face_jacobian)