DRAFT: Accurate Spherical Operations User Guide

docs/user-guide/accurate-helpers.ipynb

@@ -0,0 +1,229 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "3e7251cd7ba2b5ba",
+   "metadata": {},
+   "source": [
+    "# Accurate Spherical Operations\n",
+    "\n",
+    "UXarray provides a suite of accurate operators specifically tailored for calculations of spherical surfaces. These operators significantly enhance the precision of computations involving normal vectors, intersections of geodesic arcs (GCAs), and intersections between geodesic arcs and constant latitude lines."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "4029f4f13fa139ee",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-08-12T21:57:34.063375Z",
+     "start_time": "2024-08-12T21:57:34.060889Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import uxarray as ux\n",
+    "import numpy as np\n",
+    "\n",
+    "from uxarray.grid.coordinates import _lonlat_rad_to_xyz"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4931b26b33df7af0",
+   "metadata": {},
+   "source": "## Normal Vector Calculation"
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6287148a4fcec66a",
+   "metadata": {},
+   "source": "## GCA - GCA Intersection"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "bf794159b317cbb2",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-08-12T21:57:34.485559Z",
+     "start_time": "2024-08-12T21:57:34.483936Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from uxarray.grid.intersections import gca_gca_intersection"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9fd090adfb976922",
+   "metadata": {},
+   "source": "### Parallel"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "64351482be087d3a",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-08-12T21:57:34.997613Z",
+     "start_time": "2024-08-12T21:57:34.988864Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([], dtype=float64)"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "gca_a = np.array(\n",
+    "    [_lonlat_rad_to_xyz(0.3 * np.pi, 0.0), _lonlat_rad_to_xyz(0.5 * np.pi, 0.0)]\n",
+    ")\n",
+    "gca_b = np.array(\n",
+    "    [_lonlat_rad_to_xyz(0.5 * np.pi, 0.0), _lonlat_rad_to_xyz(-0.5 * np.pi - 0.01, 0.0)]\n",
+    ")\n",
+    "\n",
+    "gca_gca_intersection(gca_a, gca_b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "245ad090216e4773",
+   "metadata": {},
+   "source": "### Perpendicular"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "80f55804148fbadb",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-08-12T21:57:36.157666Z",
+     "start_time": "2024-08-12T21:57:36.150888Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([-0.98480775,  0.17364818, -0.        ])"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "gca_a = np.array(\n",
+    "    [\n",
+    "        _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(0.0)),\n",
+    "        _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(10.0)),\n",
+    "    ]\n",
+    ")\n",
+    "gca_b = np.array(\n",
+    "    [\n",
+    "        _lonlat_rad_to_xyz(*[0.5 * np.pi, 0.0]),\n",
+    "        _lonlat_rad_to_xyz(*[-0.5 * np.pi - 0.01, 0.0]),\n",
+    "    ]\n",
+    ")\n",
+    "gca_gca_intersection(gca_a, gca_b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "68ece4b2313bbce0",
+   "metadata": {},
+   "source": "## GCA  - Constant Latitude Intersection"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "f2a9ce9d59c2fbdf",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-08-12T21:57:37.606276Z",
+     "start_time": "2024-08-12T21:57:37.603195Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from uxarray.grid.intersections import gca_const_lat_intersection"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ddba01e5f16ff101",
+   "metadata": {},
+   "source": "## Point Within GCA"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "cb34238a0783d88",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-08-12T21:57:38.090111Z",
+     "start_time": "2024-08-12T21:57:38.086859Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from uxarray.grid.arcs import point_within_gca"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7009ba8e123ed53d",
+   "metadata": {},
+   "source": "## Spherical Bounding Box"
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7425fb621947a531",
+   "metadata": {},
+   "source": [
+    "## Takeaways\n",
+    "\n",
+    "These advancements address long-standing geoscience issues, including the degeneration of closely positioned points, inconsistent intersection points, and other geometry problems caused by error propagation. As a result, our work significantly enhances the reliability of geospatial analysis. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3c9a45b50bdae834",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/userguide.rst

@@ -58,6 +58,9 @@ These user guides provide detailed explanations of the core functionality in UXa
 `Face Area Calculations <user-guide/area_calc.ipynb>`_
  Methods for computing the area of each face
 
+`Accurate Spherical Operators <user-guide/accurate-helpers.ipynb>`_
+ SHORT DESCRIPTION TODO
+
 Supplementary Guides
 --------------------
 
@@ -80,7 +83,8 @@ These user guides provide additional detail about specific features in UXarray.
    user-guide/advanced-plotting.ipynb
    user-guide/subset.ipynb
    user-guide/topological-aggregations.ipynb
-   user-guide/area_calc.ipynb
    user-guide/holoviz.ipynb
    user-guide/remapping.ipynb
    user-guide/tree_structures.ipynb
+   user-guide/area_calc.ipynb
+   user-guide/accurate-helpers.ipynb

test/test_intersections.py

@@ -6,7 +6,7 @@
 # from uxarray.grid.coordinates import node_lonlat_rad_to_xyz, node_xyz_to_lonlat_rad
 
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad
-from uxarray.grid.intersections import gca_gca_intersection, gca_constLat_intersection
+from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection
 
 
 class TestGCAGCAIntersection(TestCase):
@@ -100,33 +100,33 @@ def test_get_GCA_GCA_intersections_perpendicular(self):
 
 class TestGCAconstLatIntersection(TestCase):
 
-    def test_GCA_constLat_intersections_antimeridian(self):
+    def test_gca_const_lat_intersections_antimeridian(self):
         GCR1_cart = np.array([
             _lonlat_rad_to_xyz(np.deg2rad(170.0),
                               np.deg2rad(89.99)),
             _lonlat_rad_to_xyz(np.deg2rad(170.0),
                               np.deg2rad(10.0))
         ])
 
-        res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(60.0)), verbose=True)
+        res = gca_const_lat_intersection(GCR1_cart, np.sin(np.deg2rad(60.0)), verbose=True)
         res_lonlat_rad = _xyz_to_lonlat_rad(*(res[0].tolist()))
         self.assertTrue(
             np.allclose(res_lonlat_rad,
                         np.array([np.deg2rad(170.0),
                                   np.deg2rad(60.0)])))
 
-    def test_GCA_constLat_intersections_empty(self):
+    def test_gca_const_lat_intersections_empty(self):
         GCR1_cart = np.array([
             _lonlat_rad_to_xyz(np.deg2rad(170.0),
                                     np.deg2rad(89.99)),
             _lonlat_rad_to_xyz(np.deg2rad(170.0),
                                     np.deg2rad(10.0))
         ])
 
-        res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(-10.0)), verbose=False)
+        res = gca_const_lat_intersection(GCR1_cart, np.sin(np.deg2rad(-10.0)), verbose=False)
         self.assertTrue(res.size == 0)
 
-    def test_GCA_constLat_intersections_two_pts(self):
+    def test_gca_const_lat_intersections_two_pts(self):
         GCR1_cart = np.array([
             _lonlat_rad_to_xyz(np.deg2rad(10.0),
                                     np.deg2rad(10)),
@@ -137,5 +137,5 @@ def test_GCA_constLat_intersections_two_pts(self):
 
         query_lat = (np.deg2rad(10.0) + max_lat) / 2.0
 
-        res = gca_constLat_intersection(GCR1_cart, np.sin(query_lat), verbose=False)
+        res = gca_const_lat_intersection(GCR1_cart, np.sin(query_lat), verbose=False)
         self.assertTrue(res.shape[0] == 2)

uxarray/grid/integrate.py

@@ -1,6 +1,6 @@
 import numpy as np
 from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
-from uxarray.grid.intersections import gca_constLat_intersection
+from uxarray.grid.intersections import gca_const_lat_intersection
 from uxarray.grid.coordinates import _xyz_to_lonlat_rad
 import pandas as pd
 
@@ -178,7 +178,7 @@ def _get_faces_constLat_intersection_info(
         # Calculate intersections (assuming a batch-capable intersection function)
         for idx, edge in enumerate(valid_edges):
             if is_GCA[idx]:
-                intersections = gca_constLat_intersection(
+                intersections = gca_const_lat_intersection(
                     edge, latitude_cart, is_directed=is_directed
                 )
                 if intersections.size == 0:

uxarray/grid/intersections.py

@@ -118,7 +118,7 @@ def gca_gca_intersection(gca1_cart, gca2_cart):
     return res
 
 
-def gca_constLat_intersection(
+def gca_const_lat_intersection(
     gca_cart, constZ, fma_disabled=False, verbose=False, is_directed=False
 ):
     """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a