Repository: UXARRAY/uxarray

Issue:
Update `.subset.bounding_box()` to use spherical face bounds
Add an option to use the spherical face bounds to subset a grid as opposed to querying the BallTree coordinates.

Discussion:

Solution patch:
diff --git a/docs/api.rst b/docs/api.rst
index 65e7abf5a..010f1c96d 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -325,6 +325,10 @@ Grid
 
    Grid.cross_section
    Grid.cross_section.constant_latitude
+   Grid.cross_section.constant_longitude
+   Grid.cross_section.constant_latitude_interval
+   Grid.cross_section.constant_longitude_interval
+
 
 UxDataArray
 ~~~~~~~~~~~
@@ -335,7 +339,9 @@ UxDataArray
 
    UxDataArray.cross_section
    UxDataArray.cross_section.constant_latitude
-
+   UxDataArray.cross_section.constant_longitude
+   UxDataArray.cross_section.constant_latitude_interval
+   UxDataArray.cross_section.constant_longitude_interval
 Remapping
 ---------
 
diff --git a/docs/user-guide/cross-sections.ipynb b/docs/user-guide/cross-sections.ipynb
index ab408b8b8..da6b5f471 100644
--- a/docs/user-guide/cross-sections.ipynb
+++ b/docs/user-guide/cross-sections.ipynb
@@ -7,21 +7,14 @@
    "source": [
     "# Cross-Sections\n",
     "\n",
-    "This section demonstrates how to extract cross-sections from an unstructured grid using UXarray, which allows the analysis and visualization across slices of grids.\n"
+    "This section demonstrates how to extract cross-sections from an unstructured grid using UXarray, which allows the analysis and visualization across slices of grids. Cross-sections can be performed directly on a `ux.Grid` object or on a `ux.UxDataArray`\n"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "b35ba4a2c30750e4",
-   "metadata": {
-    "ExecuteTime": {
-     "start_time": "2024-10-14T16:39:35.957687Z"
-    },
-    "jupyter": {
-     "is_executing": true
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "import uxarray as ux\n",
@@ -33,24 +26,11 @@
     "projection = ccrs.Robinson()"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "id": "395a3db7-495c-4cff-b733-06bbe522a604",
-   "metadata": {},
-   "source": [
-    "## Data"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "b4160275c09fe6b0",
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2024-10-09T17:50:51.217211Z",
-     "start_time": "2024-10-09T17:50:50.540946Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "base_path = \"../../test/meshfiles/ugrid/outCSne30/\"\n",
@@ -66,16 +46,6 @@
     ")"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "c844a3b5-ef71-41c8-a0ae-02ea937801d6",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "uxds.uxgrid.normalize_cartesian_coordinates()"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "a7a40958-0a4d-47e4-9e38-31925261a892",
@@ -83,7 +53,7 @@
    "source": [
     "## Constant Latitude\n",
     "\n",
-    "Cross-sections along constant latitude lines can be obtained using the ``.cross_section.constant_latitude`` method, available for both ``ux.Grid`` and ``ux.DataArray`` objects. This functionality allows users to extract and analyze slices of data at specified latitudes, providing insights into variations along horizontal sections of the grid.\n"
+    "Cross-sections along constant latitude lines can be obtained by using the  ``.cross_section.constant_latitude(lat)`` method. The sliced grid will be made up of the faces that contain at least one edge that intersects with a line of constant latitude.\n"
    ]
   },
   {
@@ -91,19 +61,14 @@
    "id": "2fbe9f6e5bb59a17",
    "metadata": {},
    "source": [
-    "For example, we can obtain a cross-section at 30 degrees latitude by doing the following:"
+    "For example, we can obtain a cross-section at 0 degrees latitude by doing the following:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "3775daa1-2f1d-4738-bab5-2b69ebd689d9",
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2024-10-09T17:50:53.093314Z",
-     "start_time": "2024-10-09T17:50:53.077719Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "lat = 0\n",
@@ -167,8 +132,7 @@
     "## Constant Longitude\n",
     "\n",
     "\n",
-    "\n",
-    "Cross-sections along constant longitude lines can be obtained using the ``.cross_section.constant_longitude`` method, available for both ``ux.Grid`` and ``ux.DataArray`` objects. \n"
+    "Cross-sections along constant longitude lines can be obtained using the ``.cross_section.constant_longitude(lon)`` method. The sliced grid will be made up of the faces that contain at least one edge that intersects with a line of constant longitude.\n"
    ]
   },
   {
@@ -205,6 +169,94 @@
     ")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "5044b8680d514fdc",
+   "metadata": {},
+   "source": [
+    "## Constant Latitude Interval\n",
+    "\n",
+    "Cross-sections between two lines of latitudes can be obtained using the ``.cross_section.constant_lats_interval(lats)`` method. The sliced grid will contain faces that are strictly between the latitude interval."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fc84e47efe2edf96",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lats = [-20, 20]\n",
+    "\n",
+    "uxda_constant_lat_interval = uxds[\"psi\"].cross_section.constant_latitude_interval(lats)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "60232626ba6c74ad",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "(\n",
+    "    uxda_constant_lat_interval.plot(\n",
+    "        rasterize=False,\n",
+    "        backend=\"bokeh\",\n",
+    "        cmap=\"inferno\",\n",
+    "        projection=projection,\n",
+    "        global_extent=True,\n",
+    "        coastline=True,\n",
+    "        title=f\"Cross Section between {lats[0]} and {lats[1]} degrees latitude\",\n",
+    "        periodic_elements=\"split\",\n",
+    "    )\n",
+    "    * gf.grid(projection=projection)\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4afa3d891a80c597",
+   "metadata": {},
+   "source": [
+    "## Constant Longitude Interval\n",
+    "\n",
+    "Cross-sections between two lines of longitude can be obtained using the ``.cross_section.constant_lons_interval(lons)`` method. The sliced grid will contain faces that are strictly between the longitude interval.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b183d15838aaf6bb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lons = [-25, 25]\n",
+    "\n",
+    "uxda_constant_lon_interval = uxds[\"psi\"].cross_section.constant_longitude_interval(lats)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "22793d56701504ce",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "(\n",
+    "    uxda_constant_lon_interval.plot(\n",
+    "        rasterize=False,\n",
+    "        backend=\"bokeh\",\n",
+    "        cmap=\"inferno\",\n",
+    "        projection=projection,\n",
+    "        global_extent=True,\n",
+    "        coastline=True,\n",
+    "        title=f\"Cross Section between {lons[0]} and {lons[1]} degrees longitude\",\n",
+    "        periodic_elements=\"split\",\n",
+    "    )\n",
+    "    * gf.grid(projection=projection)\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "54d9eff1-67f1-4691-a3b0-1ee0c874c98f",
@@ -219,7 +271,7 @@
    "metadata": {},
    "source": [
     "```{warning}\n",
-    "Arbitrary great circle arc cross sections are not yet supported.\n",
+    "Arbitrary great circle arc cross sections are not yet implemented.\n",
     "```"
    ]
   }
@@ -240,7 +292,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.8"
+   "version": "3.12.3"
   }
  },
  "nbformat": 4,
diff --git a/docs/user-guide/subset.ipynb b/docs/user-guide/subset.ipynb
index c1e7a97df..3ec66f3f2 100644
--- a/docs/user-guide/subset.ipynb
+++ b/docs/user-guide/subset.ipynb
@@ -256,20 +256,6 @@
     "lat_bounds = (41.8781 - 2, 41.8781 + 2)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "collapsed": false,
-    "jupyter": {
-     "outputs_hidden": false
-    }
-   },
-   "source": [
-    "Our bounding box ensures that the coordinates of our select element (`nodes`, `edge_centers`, or `face_centers`) are within the defined bounding box range.\n",
-    "\n",
-    "Below is an example using the corner nodes for our subset."
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -282,7 +268,8 @@
    "outputs": [],
    "source": [
     "bbox_subset_nodes = uxds[\"relhum_200hPa\"][0].subset.bounding_box(\n",
-    "    lon_bounds, lat_bounds, element=\"nodes\"\n",
+    "    lon_bounds,\n",
+    "    lat_bounds,\n",
     ")\n",
     "\n",
     "bbox_subset_nodes.plot(\n",
@@ -294,58 +281,6 @@
     ") * features"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "collapsed": false,
-    "jupyter": {
-     "outputs_hidden": false
-    }
-   },
-   "source": [
-    "And similarly using the face center coordinates."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": false,
-    "jupyter": {
-     "outputs_hidden": false
-    }
-   },
-   "outputs": [],
-   "source": [
-    "bbox_subset_faces = uxds[\"relhum_200hPa\"][0].subset.bounding_box(\n",
-    "    lon_bounds, lat_bounds, element=\"face centers\"\n",
-    ")\n",
-    "\n",
-    "bbox_subset_faces.plot(\n",
-    "    rasterize=True,\n",
-    "    periodic_elements=\"exclude\",\n",
-    "    clim=clim,\n",
-    "    title=\"Bounding Box Subset (Face Center Query)\",\n",
-    "    **plot_opts,\n",
-    ") * features"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "collapsed": false,
-    "jupyter": {
-     "outputs_hidden": false
-    }
-   },
-   "source": [
-    "While the bounding box is generally the same, you will notice differences along the border depending on which element is used to query.\n",
-    "\n",
-    ":::{note}\n",
-    "Specifying which element to query (i.e. ``nodes``, ``edgecenters``, or ``face centers``) is supported by all subsetting methods.\n",
-    ":::"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {
diff --git a/uxarray/cross_sections/dataarray_accessor.py b/uxarray/cross_sections/dataarray_accessor.py
index 52599c7a4..43365b79f 100644
--- a/uxarray/cross_sections/dataarray_accessor.py
+++ b/uxarray/cross_sections/dataarray_accessor.py
@@ -1,7 +1,7 @@
 from __future__ import annotations
 
 
-from typing import TYPE_CHECKING, Union, List, Set
+from typing import TYPE_CHECKING, Union, List, Set, Tuple
 
 if TYPE_CHECKING:
     pass
@@ -17,8 +17,10 @@ def __repr__(self):
         prefix = "<uxarray.UxDataArray.cross_section>\n"
         methods_heading = "Supported Methods:\n"
 
-        methods_heading += "  * constant_latitude(lat)\n"
-        methods_heading += "  * constant_longitude(lon)\n"
+        methods_heading += "  * constant_latitude(lat, inverse_indices)\n"
+        methods_heading += "  * constant_longitude(lon, inverse_indices)\n"
+        methods_heading += "  * constant_latitude_interval(lats, inverse_indices)\n"
+        methods_heading += "  * constant_longitude_interval(lons, inverse_indices)\n"
 
         return prefix + methods_heading
 
@@ -34,8 +36,10 @@ def constant_latitude(
             The latitude at which to extract the cross-section, in degrees.
             Must be between -90.0 and 90.0
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
 
         Returns
         -------
@@ -79,8 +83,10 @@ def constant_longitude(
             The latitude at which to extract the cross-section, in degrees.
             Must be between -180.0 and 180.0
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
 
         Returns
         -------
@@ -114,13 +120,103 @@ def constant_longitude(
 
         return self.uxda.isel(n_face=faces, inverse_indices=inverse_indices)
 
-    def gca(self, *args, **kwargs):
-        raise NotImplementedError
+    def constant_latitude_interval(
+        self,
+        lats: Tuple[float, float],
+        inverse_indices: Union[List[str], Set[str], bool] = False,
+    ):
+        """Extracts a cross-section of data by selecting all faces that
+        are within a specified latitude interval.
 
-    def bounded_latitude(self, *args, **kwargs):
-        raise NotImplementedError
+        Parameters
+        ----------
+        lats : Tuple[float, float]
+            The latitude interval (min_lat, max_lat) at which to extract the cross-section,
+            in degrees. Values must be between -90.0 and 90.0
+        return_face_indices : bool, optional
+            If True, also returns the indices of the faces that intersect with the
+            latitude interval.
+        inverse_indices : Union[List[str], Set[str], bool], optional
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
+
+        Returns
+        -------
+        uxarray.UxDataArray
+            A subset of the original data array containing only the faces that are within a specified latitude interval.
+
+        Raises
+        ------
+        ValueError
+            If no faces are found within the specified latitude interval.
+
+        Examples
+        --------
+        >>> # Extract data between 30°S and 30°N latitude
+        >>> cross_section = uxda.cross_section.constant_latitude_interval(
+        ...     lats=(-30.0, 30.0)
+        ... )
+
+
+        Notes
+        -----
+        The initial execution time may be significantly longer than subsequent runs
+        due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching.
+        """
+        faces = self.uxda.uxgrid.get_faces_between_latitudes(lats)
+
+        return self.uxda.isel(n_face=faces, inverse_indices=inverse_indices)
 
-    def bounded_longitude(self, *args, **kwargs):
+    def constant_longitude_interval(
+        self,
+        lons: Tuple[float, float],
+        inverse_indices: Union[List[str], Set[str], bool] = False,
+    ):
+        """Extracts a cross-section of data by selecting all faces are within a specifed longitude interval.
+
+        Parameters
+        ----------
+        lons : Tuple[float, float]
+            The longitude interval (min_lon, max_lon) at which to extract the cross-section,
+            in degrees. Values must be between -180.0 and 180.0
+        return_face_indices : bool, optional
+            If True, also returns the indices of the faces that intersect are within a specifed longitude interval.
+        inverse_indices : Union[List[str], Set[str], bool], optional
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
+
+        Returns
+        -------
+        uxarray.UxDataArray
+            A subset of the original data array containing only the faces that intersect
+            with the specified longitude interval.
+
+        Raises
+        ------
+        ValueError
+            If no faces are found within the specified longitude interval.
+
+        Examples
+        --------
+        >>> # Extract data between 0° and 45° longitude
+        >>> cross_section = uxda.cross_section.constant_longitude_interval(
+        ...     lons=(0.0, 45.0)
+        ... )
+
+        Notes
+        -----
+        The initial execution time may be significantly longer than subsequent runs
+        due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching.
+        """
+        faces = self.uxda.uxgrid.get_faces_between_longitudes(lons)
+
+        return self.uxda.isel(n_face=faces, inverse_indices=inverse_indices)
+
+    def gca(self, *args, **kwargs):
         raise NotImplementedError
 
     def gca_gca(self, *args, **kwargs):
diff --git a/uxarray/cross_sections/grid_accessor.py b/uxarray/cross_sections/grid_accessor.py
index 76485fbda..921a6060a 100644
--- a/uxarray/cross_sections/grid_accessor.py
+++ b/uxarray/cross_sections/grid_accessor.py
@@ -1,7 +1,6 @@
 from __future__ import annotations
 
-
-from typing import TYPE_CHECKING, Union, List, Set
+from typing import TYPE_CHECKING, Union, List, Set, Tuple
 
 if TYPE_CHECKING:
     from uxarray.grid import Grid
@@ -17,8 +16,14 @@ def __repr__(self):
         prefix = "<uxarray.Grid.cross_section>\n"
         methods_heading = "Supported Methods:\n"
 
-        methods_heading += "  * constant_latitude(lat, return_face_indices)\n"
-        methods_heading += "  * constant_longitude(lon, return_face_indices)\n"
+        methods_heading += (
+            "  * constant_latitude(lat, return_face_indices, inverse_indices)\n"
+        )
+        methods_heading += (
+            "  * constant_longitude(lon, return_face_indices, inverse_indices)\n"
+        )
+        methods_heading += "  * constant_latitude_interval(lats, return_face_indices, inverse_indices)\n"
+        methods_heading += "  * constant_longitude_interval(lons, return_face_indices, inverse_indices)\n"
         return prefix + methods_heading
 
     def constant_latitude(
@@ -32,20 +37,25 @@ def constant_latitude(
 
         Parameters
         ----------
-        lat : float
-            The latitude at which to extract the cross-section, in degrees.
+        lon : float
+            The longitude at which to extract the cross-section, in degrees.
             Must be between -90.0 and 90.0
         return_face_indices : bool, optional
-            If True, also returns the indices of the faces that intersect with the line of constant latitude.
+            If True, also returns the indices of the faces that intersect with the
+            line of constant latitude.
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face  indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
 
         Returns
         -------
         uxarray.Grid
             A subset of the original grid containing only the faces that intersect
             with the specified latitude.
+        Tuple[uxarray.Grid, numpy.ndarray], optional
+            If return_face_indices=True, returns a tuple of (grid_subset, face_indices)
 
         Raises
         ------
@@ -54,8 +64,12 @@ def constant_latitude(
 
         Examples
         --------
-        >>> # Extract data at 15.5°S latitude
-        >>> cross_section = grid.cross_section.constant_latitude(lat=-15.5)
+        >>> # Extract grid at 25° latitude
+        >>> cross_section = grid.cross_section.constant_latitude(lat=25.0)
+        >>> # With face indices
+        >>> cross_section, faces = grid.cross_section.constant_latitude(
+        ...     lat=25.0, return_face_indices=True
+        ... )
 
         Notes
         -----
@@ -92,18 +106,23 @@ def constant_longitude(
         ----------
         lon : float
             The longitude at which to extract the cross-section, in degrees.
-            Must be between -90.0 and 90.0
+            Must be between -180.0 and 180.0
         return_face_indices : bool, optional
-            If True, also returns the indices of the faces that intersect with the line of constant longitude.
+            If True, also returns the indices of the faces that intersect with the
+            line of constant longitude.
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face  indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
 
         Returns
         -------
         uxarray.Grid
             A subset of the original grid containing only the faces that intersect
             with the specified longitude.
+        Tuple[uxarray.Grid, numpy.ndarray], optional
+            If return_face_indices=True, returns a tuple of (grid_subset, face_indices)
 
         Raises
         ------
@@ -112,8 +131,12 @@ def constant_longitude(
 
         Examples
         --------
-        >>> # Extract data at 0° longitude
-        >>> cross_section = grid.cross_section.constant_latitude(lon=0.0)
+        >>> # Extract grid at 0° longitude (Prime Meridian)
+        >>> cross_section = grid.cross_section.constant_longitude(lon=0.0)
+        >>> # With face indices
+        >>> cross_section, faces = grid.cross_section.constant_longitude(
+        ...     lon=0.0, return_face_indices=True
+        ... )
 
         Notes
         -----
@@ -139,11 +162,128 @@ def constant_longitude(
     def gca(self, *args, **kwargs):
         raise NotImplementedError
 
-    def bounded_latitude(self, *args, **kwargs):
-        raise NotImplementedError
+    def constant_latitude_interval(
+        self,
+        lats: Tuple[float, float],
+        return_face_indices: bool = False,
+        inverse_indices: Union[List[str], Set[str], bool] = False,
+    ):
+        """Extracts a cross-section of the grid by selecting all faces that
+        are within a specified latitude interval.
 
-    def bounded_longitude(self, *args, **kwargs):
-        raise NotImplementedError
+        Parameters
+        ----------
+        lats : Tuple[float, float]
+            The latitude interval (min_lat, max_lat) at which to extract the cross-section,
+            in degrees. Values must be between -90.0 and 90.0
+        return_face_indices : bool, optional
+            If True, also returns the indices of the faces that intersect with the
+            latitude interval.
+        inverse_indices : Union[List[str], Set[str], bool], optional
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
+
+        Returns
+        -------
+        uxarray.Grid
+            A subset of the original grid containing only the faces that are within a specified latitude interval.
+        Tuple[uxarray.Grid, numpy.ndarray], optional
+            If return_face_indices=True, returns a tuple of (grid_subset, face_indices)
+
+        Raises
+        ------
+        ValueError
+            If no faces are found within the specified latitude interval.
+
+        Examples
+        --------
+        >>> # Extract grid between 30°S and 30°N latitude
+        >>> cross_section = grid.cross_section.constant_latitude_interval(
+        ...     lats=(-30.0, 30.0)
+        ... )
+        >>> # With face indices
+        >>> cross_section, faces = grid.cross_section.constant_latitude_interval(
+        ...     lats=(-30.0, 30.0), return_face_indices=True
+        ... )
+
+        Notes
+        -----
+        The initial execution time may be significantly longer than subsequent runs
+        due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching.
+        """
+        faces = self.uxgrid.get_faces_between_latitudes(lats)
+
+        grid_between_lats = self.uxgrid.isel(
+            n_face=faces, inverse_indices=inverse_indices
+        )
+
+        if return_face_indices:
+            return grid_between_lats, faces
+        else:
+            return grid_between_lats
+
+    def constant_longitude_interval(
+        self,
+        lons: Tuple[float, float],
+        return_face_indices: bool = False,
+        inverse_indices: Union[List[str], Set[str], bool] = False,
+    ):
+        """Extracts a cross-section of the grid by selecting all faces are within a specifed longitude interval.
+
+        Parameters
+        ----------
+        lons : Tuple[float, float]
+            The longitude interval (min_lon, max_lon) at which to extract the cross-section,
+            in degrees. Values must be between -180.0 and 180.0
+        return_face_indices : bool, optional
+            If True, also returns the indices of the faces that intersect are within a specifed longitude interval.
+        inverse_indices : Union[List[str], Set[str], bool], optional
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
+
+        Returns
+        -------
+        uxarray.Grid
+            A subset of the original grid containing only the faces that intersect
+            with the specified longitude interval.
+        Tuple[uxarray.Grid, numpy.ndarray], optional
+            If return_face_indices=True, returns a tuple of (grid_subset, face_indices)
+
+        Raises
+        ------
+        ValueError
+            If no faces are found within the specified longitude interval.
+
+        Examples
+        --------
+        >>> # Extract grid between 0° and 45° longitude
+        >>> cross_section = grid.cross_section.constant_longitude_interval(
+        ...     lons=(0.0, 45.0)
+        ... )
+        >>> # With face indices
+        >>> cross_section, faces = grid.cross_section.constant_longitude_interval(
+        ...     lons=(0.0, 45.0), return_face_indices=True
+        ... )
+
+        Notes
+        -----
+        The initial execution time may be significantly longer than subsequent runs
+        due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching.
+        """
+        faces = self.uxgrid.get_faces_between_longitudes(lons)
+
+        grid_between_lons = self.uxgrid.isel(
+            n_face=faces, inverse_indices=inverse_indices
+        )
+
+        if return_face_indices:
+            return grid_between_lons, faces
+        else:
+            return grid_between_lons
 
     def gca_gca(self, *args, **kwargs):
         raise NotImplementedError
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
index 5c11c8281..45e00ba42 100644
--- a/uxarray/grid/coordinates.py
+++ b/uxarray/grid/coordinates.py
@@ -5,7 +5,7 @@
 
 from uxarray.conventions import ugrid
 
-from numba import njit
+from numba import njit, prange
 from uxarray.constants import ERROR_TOLERANCE
 from typing import Union
 
@@ -305,6 +305,7 @@ def _populate_face_centroids(grid, repopulate=False):
         )
 
 
+@njit(cache=True, parallel=True)
 def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face):
     """Constructs the xyz centroid coordinate for each face using Cartesian
     Averaging.
@@ -330,14 +331,21 @@ def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_fa
     centroid_x = np.zeros((face_nodes.shape[0]), dtype=np.float64)
     centroid_y = np.zeros((face_nodes.shape[0]), dtype=np.float64)
     centroid_z = np.zeros((face_nodes.shape[0]), dtype=np.float64)
+    n_face = n_nodes_per_face.shape[0]
 
-    for face_idx, n_max_nodes in enumerate(n_nodes_per_face):
-        # Compute Cartesian Average
-        centroid_x[face_idx] = np.mean(node_x[face_nodes[face_idx, 0:n_max_nodes]])
-        centroid_y[face_idx] = np.mean(node_y[face_nodes[face_idx, 0:n_max_nodes]])
-        centroid_z[face_idx] = np.mean(node_z[face_nodes[face_idx, 0:n_max_nodes]])
+    for i_face in prange(n_face):
+        n_max_nodes = n_nodes_per_face[i_face]
 
-    return _normalize_xyz(centroid_x, centroid_y, centroid_z)
+        x = np.mean(node_x[face_nodes[i_face, 0:n_max_nodes]])
+        y = np.mean(node_y[face_nodes[i_face, 0:n_max_nodes]])
+        z = np.mean(node_z[face_nodes[i_face, 0:n_max_nodes]])
+
+        x, y, z = _normalize_xyz_scalar(x, y, z)
+
+        centroid_x[i_face] = x
+        centroid_y[i_face] = y
+        centroid_z[i_face] = z
+    return centroid_x, centroid_y, centroid_z
 
 
 def _welzl_recursive(points, boundary, R):
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 7f51c47c5..aa93192d2 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -11,6 +11,7 @@
     Union,
     List,
     Set,
+    Tuple,
 )
 
 # reader and writer imports
@@ -81,6 +82,8 @@
     constant_lon_intersections_no_extreme,
     constant_lat_intersections_face_bounds,
     constant_lon_intersections_face_bounds,
+    faces_within_lon_bounds,
+    faces_within_lat_bounds,
 )
 
 from spatialpandas import GeoDataFrame
@@ -203,8 +206,7 @@ def __init__(
         self._ds.assign_attrs({"source_grid_spec": self.source_grid_spec})
         self._is_subset = is_subset
 
-        if inverse_indices is not None:
-            self._inverse_indices = inverse_indices
+        self._inverse_indices = inverse_indices
 
         # cached parameters for GeoDataFrame conversions
         self._gdf_cached_parameters = {
@@ -2415,3 +2417,35 @@ def get_faces_at_constant_longitude(self, lon: float):
 
         faces = constant_lon_intersections_face_bounds(lon, self.face_bounds_lon.values)
         return faces
+
+    def get_faces_between_longitudes(self, lons: Tuple[float, float]):
+        """Identifies the indices of faces that are strictly between two lines of constant longitude.
+
+        Parameters
+        ----------
+        lons: Tuple[float, float]
+            A tuple of longitudes that define that minimum and maximum longitude.
+
+        Returns
+        -------
+        faces : numpy.ndarray
+            An array of face indices that are strictly between two lines of constant longitude.
+
+        """
+        return faces_within_lon_bounds(lons, self.face_bounds_lon.values)
+
+    def get_faces_between_latitudes(self, lats: Tuple[float, float]):
+        """Identifies the indices of faces that are strictly between two lines of constant latitude.
+
+        Parameters
+        ----------
+        lats: Tuple[float, float
+            A tuple of latitudes that define that minimum and maximum latitudes.
+
+        Returns
+        -------
+        faces : numpy.ndarray
+            An array of face indices that are strictly between two lines of constant latitude.
+
+        """
+        return faces_within_lat_bounds(lats, self.face_bounds_lat.values)
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 4550598bb..cd9b702ed 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -129,28 +129,26 @@ def constant_lon_intersections_no_extreme(lon, edge_node_x, edge_node_y, n_edge)
     return np.unique(intersecting_edges)
 
 
-# @njit(cache=True)
-def constant_lat_intersections_face_bounds(lat, face_bounds_lat):
-    """Identifies the candidate faces on a grid that intersect with a given
-    constant latitude.
-
-    This function checks whether the specified latitude, `lat`, in degrees lies within
-    the latitude bounds of grid faces, defined by `face_min_lat_rad` and `face_max_lat_rad`,
-    which are given in radians. The function returns the indices of the faces where the
-    latitude is within these bounds.
+@njit(cache=True)
+def constant_lat_intersections_face_bounds(lat: float, face_bounds_lat: np.ndarray):
+    """
+    Identify candidate faces that intersect with a given constant latitude line.
 
     Parameters
     ----------
     lat : float
         The latitude in degrees for which to find intersecting faces.
-    TODO:
+    face_bounds_lat : numpy.ndarray
+        A 2D array of shape (n_faces, 2), where each row represents the latitude
+        bounds of a face. The first element of each row is the minimum latitude
+        and the second element is the maximum latitude of the face.
 
     Returns
     -------
     candidate_faces : numpy.ndarray
-        A 1D array containing the indices of the faces that intersect with the given latitude.
+        A 1D array of integers containing the indices of the faces that intersect
+        the given latitude.
     """
-
     face_bounds_lat_min = face_bounds_lat[:, 0]
     face_bounds_lat_max = face_bounds_lat[:, 1]
 
@@ -160,27 +158,25 @@ def constant_lat_intersections_face_bounds(lat, face_bounds_lat):
 
 
 @njit(cache=True)
-def constant_lon_intersections_face_bounds(lon, face_bounds_lon):
-    """Identifies the candidate faces on a grid that intersect with a given
-    constant longitude.
-
-    This function checks whether the specified longitude, `lon`, in degrees lies within
-    the longitude bounds of grid faces, defined by `face_min_lon_rad` and `face_max_lon_rad`,
-    which are given in radians. The function returns the indices of the faces where the
-    longitude is within these bounds.
+def constant_lon_intersections_face_bounds(lon: float, face_bounds_lon: np.ndarray):
+    """
+    Identify candidate faces that intersect with a given constant longitude line.
 
     Parameters
     ----------
     lon : float
         The longitude in degrees for which to find intersecting faces.
-    TODO:
+    face_bounds_lon : numpy.ndarray
+        A 2D array of shape (n_faces, 2), where each row represents the longitude
+        bounds of a face. The first element of each row is the minimum longitude
+        and the second element is the maximum longitude of the face.
 
     Returns
     -------
     candidate_faces : numpy.ndarray
-        A 1D array containing the indices of the faces that intersect with the given longitude.
+        A 1D array of integers containing the indices of the faces that intersect
+        the given longitude.
     """
-
     face_bounds_lon_min = face_bounds_lon[:, 0]
     face_bounds_lon_max = face_bounds_lon[:, 1]
     n_face = face_bounds_lon.shape[0]
@@ -191,16 +187,116 @@ def constant_lon_intersections_face_bounds(lon, face_bounds_lon):
         cur_face_bounds_lon_max = face_bounds_lon_max[i]
 
         if cur_face_bounds_lon_min < cur_face_bounds_lon_max:
-            if (lon >= cur_face_bounds_lon_min) & (lon <= cur_face_bounds_lon_max):
+            if (lon >= cur_face_bounds_lon_min) and (lon <= cur_face_bounds_lon_max):
                 candidate_faces.append(i)
         else:
             # antimeridian case
-            if (lon >= cur_face_bounds_lon_min) | (lon <= cur_face_bounds_lon_max):
+            if (lon >= cur_face_bounds_lon_min) or (lon <= cur_face_bounds_lon_max):
                 candidate_faces.append(i)
 
     return np.array(candidate_faces, dtype=INT_DTYPE)
 
 
+@njit(cache=True)
+def faces_within_lon_bounds(lons, face_bounds_lon):
+    """
+    Identify candidate faces that lie within a specified longitudinal interval.
+
+    Parameters
+    ----------
+    lons : tuple or list of length 2
+        A pair (min_lon, max_lon) specifying the query interval. If `min_lon <= max_lon`,
+        the interval is [min_lon, max_lon]. If `min_lon > max_lon`, the interval
+        crosses the antimeridian and should be interpreted as [min_lon, 180] U [-180, max_lon].
+    face_bounds_lon : numpy.ndarray
+        A 2D array of shape (n_faces, 2), where each row represents the longitude bounds
+        of a face. The first element is the minimum longitude and the second is the maximum
+        longitude for that face. Bounds may cross the antimeridian.
+
+    Returns
+    -------
+    candidate_faces : numpy.ndarray
+        A 1D array of integers containing the indices of the faces whose longitude bounds
+        overlap with the specified interval.
+    """
+    face_bounds_lon_min = face_bounds_lon[:, 0]
+    face_bounds_lon_max = face_bounds_lon[:, 1]
+    n_face = face_bounds_lon.shape[0]
+
+    min_lon, max_lon = lons
+
+    # For example, a query of (160, -160) would cross the antimeridian
+    antimeridian = min_lon > max_lon
+
+    candidate_faces = []
+    for i in range(n_face):
+        cur_face_min = face_bounds_lon_min[i]
+        cur_face_max = face_bounds_lon_max[i]
+
+        # Check if the face itself crosses the antimeridian
+        face_crosses_antimeridian = cur_face_min > cur_face_max
+
+        if not antimeridian:
+            # Normal case: min_lon <= max_lon
+            # Face must be strictly contained within [min_lon, max_lon]
+            if not face_crosses_antimeridian:
+                if (cur_face_min >= min_lon) and (cur_face_max <= max_lon):
+                    candidate_faces.append(i)
+            else:
+                # If face crosses antimeridian, it cannot be strictly contained
+                # in a non-antimeridian query interval
+                continue
+        else:
+            # Antimeridian case: interval crosses the -180/180 boundary
+            # The query interval is effectively [min_lon, 180] U [-180, max_lon]
+
+            if face_crosses_antimeridian:
+                # If face crosses antimeridian, check if it's contained in the full query range
+                if (cur_face_min >= min_lon) and (cur_face_max <= max_lon):
+                    candidate_faces.append(i)
+            else:
+                # For non-crossing faces, check if they're strictly contained in either part
+                contained_part1 = (cur_face_min >= min_lon) and (cur_face_max <= 180)
+                contained_part2 = (cur_face_min >= -180) and (cur_face_max <= max_lon)
+
+                if contained_part1 or contained_part2:
+                    candidate_faces.append(i)
+
+    return np.array(candidate_faces, dtype=INT_DTYPE)
+
+
+@njit(cache=True)
+def faces_within_lat_bounds(lats, face_bounds_lat):
+    """
+    Identify candidate faces that lie within a specified latitudinal interval.
+
+    Parameters
+    ----------
+    lats : tuple or list of length 2
+        A pair (min_lat, max_lat) specifying the query interval. All returned faces
+        must be fully contained within this interval.
+    face_bounds_lat : numpy.ndarray
+        A 2D array of shape (n_faces, 2), where each row represents the latitude
+        bounds of a face. The first element is the minimum latitude and the second
+        is the maximum latitude for that face.
+
+    Returns
+    -------
+    candidate_faces : numpy.ndarray
+        A 1D array of integers containing the indices of the faces whose latitude
+        bounds lie completely within the specified interval.
+    """
+
+    min_lat, max_lat = lats
+
+    face_bounds_lat_min = face_bounds_lat[:, 0]
+    face_bounds_lat_max = face_bounds_lat[:, 1]
+
+    within_bounds = (face_bounds_lat_max <= max_lat) & (face_bounds_lat_min >= min_lat)
+    candidate_faces = np.where(within_bounds)[0]
+    return candidate_faces
+
+
 def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz):
     gca_a_xyz = np.asarray(gca_a_xyz)
     gca_b_xyz = np.asarray(gca_b_xyz)
@@ -321,7 +417,6 @@ def gca_const_lat_intersection(gca_cart, constZ, fma_disabled=True, verbose=Fals
         return res
 
     # If the constant latitude is not the same as the GCA endpoints, calculate the intersection point
-    # TODO:
     lat_min = _extreme_gca_latitude_cartesian(gca_cart, extreme_type="min")
     lat_max = _extreme_gca_latitude_cartesian(gca_cart, extreme_type="max")
 
diff --git a/uxarray/subset/dataarray_accessor.py b/uxarray/subset/dataarray_accessor.py
index 2d966c587..1775baa68 100644
--- a/uxarray/subset/dataarray_accessor.py
+++ b/uxarray/subset/dataarray_accessor.py
@@ -19,10 +19,12 @@ def __repr__(self):
         prefix = "<uxarray.UxDataArray.subset>\n"
         methods_heading = "Supported Methods:\n"
 
-        methods_heading += "  * nearest_neighbor(center_coord, k, element, **kwargs)\n"
-        methods_heading += "  * bounding_circle(center_coord, r, element, **kwargs)\n"
+        methods_heading += "  * nearest_neighbor(center_coord, k, element, inverse_indices, **kwargs)\n"
         methods_heading += (
-            "  * bounding_box(lon_bounds, lat_bounds, element, method, **kwargs)\n"
+            "  * bounding_circle(center_coord, r, element, inverse_indices, **kwargs)\n"
+        )
+        methods_heading += (
+            "  * bounding_box(lon_bounds, lat_bounds, inverse_indices, **kwargs)\n"
         )
 
         return prefix + methods_heading
@@ -31,10 +33,7 @@ def bounding_box(
         self,
         lon_bounds: Union[Tuple, List, np.ndarray],
         lat_bounds: Union[Tuple, List, np.ndarray],
-        element: Optional[str] = "nodes",
-        method: Optional[str] = "coords",
         inverse_indices: Union[List[str], Set[str], bool] = False,
-        **kwargs,
     ):
         """Subsets an unstructured grid between two latitude and longitude
         points which form a bounding box.
@@ -55,14 +54,14 @@ def bounding_box(
         element: str
             Element for use with `coords` comparison, one of `nodes`, `face centers`, or `edge centers`
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
         """
-        grid = self.uxda.uxgrid.subset.bounding_box(
-            lon_bounds, lat_bounds, element, method, inverse_indices=inverse_indices
-        )
+        grid = self.uxda.uxgrid.subset.bounding_box(lon_bounds, lat_bounds)
 
-        return self.uxda._slice_from_grid(grid)
+        return self.uxda._slice_from_grid(grid, inverse_indices=inverse_indices)
 
     def bounding_circle(
         self,
@@ -84,8 +83,10 @@ def bounding_circle(
         element: str
             Element for use with `coords` comparison, one of `nodes`, `face centers`, or `edge centers`
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
         """
         grid = self.uxda.uxgrid.subset.bounding_circle(
             center_coord, r, element, inverse_indices=inverse_indices, **kwargs
@@ -112,8 +113,10 @@ def nearest_neighbor(
         element: str
             Element for use with `coords` comparison, one of `nodes`, `face centers`, or `edge centers`
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
         """
 
         grid = self.uxda.uxgrid.subset.nearest_neighbor(
diff --git a/uxarray/subset/grid_accessor.py b/uxarray/subset/grid_accessor.py
index a504179f1..e42720db7 100644
--- a/uxarray/subset/grid_accessor.py
+++ b/uxarray/subset/grid_accessor.py
@@ -19,22 +19,19 @@ def __repr__(self):
         prefix = "<uxarray.Grid.subset>\n"
         methods_heading = "Supported Methods:\n"
 
-        methods_heading += "  * nearest_neighbor(center_coord, k, element, **kwargs)\n"
-        methods_heading += "  * bounding_circle(center_coord, r, element, **kwargs)\n"
+        methods_heading += "  * nearest_neighbor(center_coord, k, element, inverse_indices, **kwargs)\n"
         methods_heading += (
-            "  * bounding_box(lon_bounds, lat_bounds, element, method, **kwargs)\n"
+            "  * bounding_circle(center_coord, r, element, inverse_indices, **kwargs)\n"
         )
+        methods_heading += "  * bounding_box(lon_bounds, lat_bounds, inverse_indices)\n"
 
         return prefix + methods_heading
 
     def bounding_box(
         self,
-        lon_bounds: Union[Tuple, List, np.ndarray],
-        lat_bounds: Union[Tuple, List, np.ndarray],
-        element: Optional[str] = "nodes",
-        method: Optional[str] = "coords",
+        lon_bounds: Tuple[float, float],
+        lat_bounds: Tuple[float, float],
         inverse_indices: Union[List[str], Set[str], bool] = False,
-        **kwargs,
     ):
         """Subsets an unstructured grid between two latitude and longitude
         points which form a bounding box.
@@ -55,64 +52,18 @@ def bounding_box(
         element: str
             Element for use with `coords` comparison, one of `nodes`, `face centers`, or `edge centers`
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
         """
 
-        if method == "coords":
-            # coordinate comparison bounding box
-
-            if element == "nodes":
-                lat, lon = self.uxgrid.node_lat.values, self.uxgrid.node_lon.values
-            elif element == "face centers":
-                lat, lon = self.uxgrid.face_lat.values, self.uxgrid.face_lon.values
-            elif element == "edge centers":
-                lat, lon = self.uxgrid.edge_lat.values, self.uxgrid.edge_lon.values
-            else:
-                raise ValueError("TODO")
-
-            if lon_bounds[0] > lon_bounds[1]:
-                # split across antimeridian
-
-                lon_indices_lhs = np.argwhere(
-                    np.logical_and(lon >= -180, lon < lon_bounds[1])
-                )
-
-                lon_indices_rhs = np.argwhere(
-                    np.logical_and(lon >= lon_bounds[0], lon < 180)
-                )
-
-                lon_indices = np.union1d(
-                    lon_indices_lhs.squeeze(), lon_indices_rhs.squeeze()
-                )
-            else:
-                # continuous bound
-
-                lon_indices = np.argwhere(
-                    np.logical_and(lon > lon_bounds[0], lon < lon_bounds[1])
-                )
-
-            lat_indices = np.argwhere(
-                np.logical_and(lat > lat_bounds[0], lat < lat_bounds[1])
-            )
-
-            # treat both indices as a set, find the intersection of both
-            indices = np.intersect1d(lat_indices, lon_indices)
-
-            if len(indices) == 0:
-                raise ValueError(
-                    f"No elements founding within the bounding box when querying {element}"
-                )
+        faces_between_lons = self.uxgrid.get_faces_between_longitudes(lon_bounds)
+        face_between_lats = self.uxgrid.get_faces_between_latitudes(lat_bounds)
 
-            if element == "nodes":
-                return self.uxgrid.isel(inverse_indices, n_node=indices)
-            elif element == "face centers":
-                return self.uxgrid.isel(inverse_indices, n_face=indices)
-            elif element == "edge centers":
-                return self.uxgrid.isel(inverse_indices, n_edge=indices)
+        faces = np.intersect1d(faces_between_lons, face_between_lats)
 
-        else:
-            raise ValueError(f"Method '{method}' not supported.")
+        return self.uxgrid.isel(n_face=faces, inverse_indices=inverse_indices)
 
     def bounding_circle(
         self,
@@ -134,8 +85,10 @@ def bounding_circle(
         element: str
             Element for use with `coords` comparison, one of `nodes`, `face centers`, or `edge centers`
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
         """
 
         coords = np.asarray(center_coord)
@@ -171,8 +124,10 @@ def nearest_neighbor(
         element: str
             Element for use with `coords` comparison, one of `nodes`, `face centers`, or `edge centers`
         inverse_indices : Union[List[str], Set[str], bool], optional
-            Indicates whether to store the original grids indices. Passing `True` stores the original face centers,
-            other reverse indices can be stored by passing any or all of the following: (["face", "edge", "node"], True)
+            Controls storage of original grid indices. Options:
+            - True: Stores original face indices
+            - List/Set of strings: Stores specified index types (valid values: "face", "edge", "node")
+            - False: No index storage (default)
         """
 
         coords = np.asarray(center_coord)


Test patch:
diff --git a/test/test_cross_sections.py b/test/test_cross_sections.py
index 0e249c57a..8eaf035da 100644
--- a/test/test_cross_sections.py
+++ b/test/test_cross_sections.py
@@ -1,8 +1,10 @@
-import uxarray as ux
-import pytest
-import numpy as np
 from pathlib import Path
 
+import numpy as np
+import pytest
+
+import uxarray as ux
+
 # Define the current path and file paths for grid and data
 current_path = Path(__file__).resolve().parent
 quad_hex_grid_path = current_path / 'meshfiles' / "ugrid" / "quad-hexagon" / 'grid.nc'
@@ -40,7 +42,6 @@ def test_constant_lon_cross_section_grid():
 
 def test_constant_lat_cross_section_uxds():
     uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path)
-    uxds.uxgrid.normalize_cartesian_coordinates()
 
     da_top_two = uxds['t2m'].cross_section.constant_latitude(lat=0.1)
     np.testing.assert_array_equal(da_top_two.data, uxds['t2m'].isel(n_face=[1, 2]).data)
@@ -56,7 +57,6 @@ def test_constant_lat_cross_section_uxds():
 
 def test_constant_lon_cross_section_uxds():
     uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path)
-    uxds.uxgrid.normalize_cartesian_coordinates()
 
     da_left_two = uxds['t2m'].cross_section.constant_longitude(lon=-0.1)
     np.testing.assert_array_equal(da_left_two.data, uxds['t2m'].isel(n_face=[0, 2]).data)
@@ -115,3 +115,56 @@ def test_constant_lat_out_of_bounds():
     )
 
     assert len(candidate_faces) == 0
+
+
+
+def test_constant_longitude_interval_grid():
+    uxgrid = ux.open_grid(quad_hex_grid_path)
+
+    grid_all_four = uxgrid.cross_section.constant_longitude_interval(lons=(-1, 1))
+
+
+    assert grid_all_four.n_face == 4
+
+    grid_left_two = uxgrid.cross_section.constant_longitude_interval(lons=(-1, 0.2))
+
+    assert grid_left_two.n_face == 2
+
+    grid_right_two = uxgrid.cross_section.constant_longitude_interval(lons=(-0.1, 1.0))
+
+    assert grid_right_two.n_face == 2
+
+def test_constant_latitude_interval_grid():
+    uxgrid = ux.open_grid(quad_hex_grid_path)
+
+    grid_all_four = uxgrid.cross_section.constant_latitude_interval(lats=(-0.3, 0.3))
+
+    assert grid_all_four.n_face == 4
+
+    grid_top_two = uxgrid.cross_section.constant_latitude_interval(lats=(-0.1, 0.3))
+    assert grid_top_two.n_face == 2
+
+    grid_bottom_two = uxgrid.cross_section.constant_latitude_interval(lats=(-0.3, 0.1))
+
+    assert grid_bottom_two.n_face == 2
+
+def test_constant_latitude_interval_uxds():
+    uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path)
+
+    da_top_two = uxds['t2m'].cross_section.constant_latitude_interval(lats=(-0.1, 0.3))
+    np.testing.assert_array_equal(da_top_two.data, uxds['t2m'].isel(n_face=[1, 2]).data)
+
+    da_bottom_two = uxds['t2m'].cross_section.constant_latitude_interval(lats=(-0.3, 0.1))
+    np.testing.assert_array_equal(da_bottom_two.data, uxds['t2m'].isel(n_face=[0, 3]).data)
+
+    da_all_four = uxds['t2m'].cross_section.constant_latitude_interval(lats=(-0.3, 0.3))
+    np.testing.assert_array_equal(da_all_four.data, uxds['t2m'].data)
+
+def test_constant_longitude_interval_uxds():
+    uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path)
+
+    da_left_two = uxds['t2m'].cross_section.constant_longitude_interval(lons=(-1, 0.2))
+    np.testing.assert_array_equal(da_left_two.data, uxds['t2m'].isel(n_face=[0, 2]).data)
+
+    da_right_two = uxds['t2m'].cross_section.constant_longitude_interval(lons=(-0.1, 1.0))
+    np.testing.assert_array_equal(da_right_two.data, uxds['t2m'].isel(n_face=[1, 3]).data)
diff --git a/test/test_subset.py b/test/test_subset.py
index f13ae4919..252719500 100644
--- a/test/test_subset.py
+++ b/test/test_subset.py
@@ -97,11 +97,10 @@ def test_grid_bounding_box_subset():
             grid = ux.open_grid(grid_path)
 
             grid_subset = grid.subset.bounding_box(bbox[0],
-                                                   bbox[1],
-                                                   element=element)
+                                                   bbox[1],)
 
             grid_subset_antimeridian = grid.subset.bounding_box(
-                bbox_antimeridian[0], bbox_antimeridian[1], element=element)
+                bbox_antimeridian[0], bbox_antimeridian[1])
 
 
 def test_uxda_isel():
@@ -132,7 +131,7 @@ def test_inverse_indices():
 
     # Test bounding box subsetting
     box = [(-10, 10), (-10, 10)]
-    subset = grid.subset.bounding_box(box[0], box[1], element="face centers", inverse_indices=True)
+    subset = grid.subset.bounding_box(box[0], box[1], inverse_indices=True)
 
     assert subset.inverse_indices is not None