Point In Face & Get Faces Containing Point

benchmarks/mpas_ocean.py

@@ -194,6 +194,28 @@ def time_const_lat(self, resolution, lat_step):
             self.uxgrid.cross_section.constant_latitude(lat)
 
 
+class PointInPolygon:
+    param_names = ['resolution']
+    params = ['480km', '120km']
+
+    def setup(self, resolution):
+        self.uxgrid = ux.open_grid(file_path_dict[resolution][0])
+        self.uxgrid.normalize_cartesian_coordinates()
+
+        _ = self.uxgrid.face_edge_connectivity
+
+        point = np.array([0.0, 0.0, 1.0])
+        res = self.uxgrid.get_faces_containing_point(point)
+
+    def teardown(self, resolution):
+        del self.uxgrid
+
+    def time_face_search(self, resolution):
+        point_xyz = np.array([self.uxgrid.face_x[0].values, self.uxgrid.face_y[0].values, self.uxgrid.face_z[0].values])
+
+        self.uxgrid.get_faces_containing_point(point=point_xyz)
+
+
 class ZonalAverage(DatasetBenchmark):
     def setup(self, resolution, *args, **kwargs):
         self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])

docs/api.rst

@@ -139,6 +139,7 @@ Descriptors
    Grid.partial_sphere_coverage
    Grid.global_sphere_coverage
    Grid.triangular
+   Grid.max_face_radius
 
 Attributes
 ~~~~~~~~~~
@@ -159,6 +160,7 @@ Methods
    Grid.calculate_total_face_area
    Grid.normalize_cartesian_coordinates
    Grid.construct_face_centers
+   Grid.get_faces_containing_point
 
 Inheritance of Xarray Functionality
 -----------------------------------

docs/user-guide/point_in_face.ipynb

@@ -0,0 +1,170 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "336fe698-c5a1-481e-9aea-9de158174a96",
+   "metadata": {},
+   "source": [
+    "# Point In Face"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a5dbe362-99b5-4c2d-a148-cd3816cb3334",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import uxarray as ux\n",
+    "import numpy as np\n",
+    "from uxarray.grid.utils import _get_cartesian_face_edge_nodes\n",
+    "from uxarray.grid.geometry import point_in_face\n",
+    "\n",
+    "grid_path = \"../../test/meshfiles/ugrid/quad-hexagon/grid.nc\"\n",
+    "grid = ux.open_grid(grid_path)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "30c4cff2-f760-46d3-9ff9-11f48a51cea8",
+   "metadata": {},
+   "source": [
+    "The function `point_in_face` determines whether or not a given face contains a point. The function parameters are as follows:\n",
+    "\n",
+    "* `edges_xyz` are the supplied edges of the face, with shape [n_edges, 2]\n",
+    "* `point_xyz` the cartesian point which is being checked. The dtype should be `np.float64`\n",
+    "* `inclusive` bool flag indicating whether to classify a point that lies on either the edges or nodes as within the face"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5e70681b-e7c3-4817-bcdc-3241c0f3e46e",
+   "metadata": {},
+   "source": [
+    "We can take a face and check to see if it's face center lies within the face."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c8c6f19a-38f2-4211-8a7b-d5f383eb5dc5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get the face edge arrays of each face in the grid\n",
+    "faces_edge_cartesian = _get_cartesian_face_edge_nodes(\n",
+    "    grid.face_node_connectivity.values,\n",
+    "    grid.n_face,\n",
+    "    grid.n_max_face_nodes,\n",
+    "    grid.node_x.values,\n",
+    "    grid.node_y.values,\n",
+    "    grid.node_z.values,\n",
+    ")\n",
+    "\n",
+    "face_center = np.array(\n",
+    "    [grid.face_x.values[0], grid.face_y.values[0], grid.face_z.values[0]],\n",
+    "    dtype=np.float64,\n",
+    ")\n",
+    "point_in_face(faces_edge_cartesian[0], face_center, inclusive=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5a656ef1-3898-4586-8ca7-3d43f7517adc",
+   "metadata": {},
+   "source": [
+    "We can do something similar but with the nodes of the face, to showcase the inclusive flag use."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7e686262-da40-4b86-84af-9ac2a10887ca",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "face_node = np.array(\n",
+    "    [grid.node_x.values[0], grid.node_y.values[0], grid.node_z.values[0]],\n",
+    "    dtype=np.float64,\n",
+    ")\n",
+    "inclusive = point_in_face(faces_edge_cartesian[0], face_node, inclusive=True)\n",
+    "non_inclusive = point_in_face(faces_edge_cartesian[0], face_node, inclusive=False)\n",
+    "\n",
+    "print(\n",
+    "    f\"Node is in the face when inclusive is set to True: {inclusive}\\nNode is not in the face when inclusive is set to False: {non_inclusive}\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c3e3ead1-40cf-44c7-84b2-ca7b4c4bbdac",
+   "metadata": {},
+   "source": [
+    "# Get Faces Containing Point"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "632c7222-154e-459a-b063-98b4260729e8",
+   "metadata": {},
+   "source": [
+    "Inside UXarray we can also retrieve any faces that contain a given point. By simply passing in a point, the function will return an array that contains the indices of the faces that the point is contained within. A point can be contained within many faces if the point is on a node or an edge. This is because many different faces may use the same node, and an edge will likely have two faces saddling it. This function considers these points as part of the faces that contain that edge or node. If an empty array is returned, this is likely because the grid is a partial and the point is within a region that has no faces. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3f2277e6-0073-4432-a424-e86efbaee1b8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Face center\n",
+    "grid.get_faces_containing_point(face_center)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "39ee9098-d639-4605-b15b-9d73085498a9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Node which is apart of 3 different faces\n",
+    "grid.get_faces_containing_point(face_node)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "daa70b10-514c-4942-ba7b-1b76bef626f9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Pole point which this partial grid does not cover\n",
+    "random_xyz = np.array([0, 0, 1], dtype=np.float64)\n",
+    "grid.get_faces_containing_point(random_xyz)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/userguide.rst

@@ -82,6 +82,8 @@ These user guides provide detailed explanations of the core functionality in UXa
 `Dual Mesh Construction <user-guide/dual-mesh.ipynb>`_
  Construct the Dual Mesh of an unstructured grid
 
+`Point in Face <user-guide/point_in_face.ipynb>`_
+ Check if a point is within a face and find the faces containing a point
 
 Supplementary Guides
 --------------------
@@ -123,3 +125,4 @@ These user guides provide additional details about specific features in UXarray.
    user-guide/healpix.ipynb
    user-guide/holoviz.ipynb
    user-guide/from_file.ipynb
+   user-guide/point_in_face.ipynb