diff --git a/docs/user-guide/accurate-helpers.ipynb b/docs/user-guide/accurate-helpers.ipynb new file mode 100644 index 000000000..3c31d2026 --- /dev/null +++ b/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 +} diff --git a/docs/userguide.rst b/docs/userguide.rst index 8d1f598ba..f533c4ec2 100644 --- a/docs/userguide.rst +++ b/docs/userguide.rst @@ -58,6 +58,9 @@ These user guides provide detailed explanations of the core functionality in UXa `Face Area Calculations `_ Methods for computing the area of each face +`Accurate Spherical Operators `_ + 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 diff --git a/test/test_intersections.py b/test/test_intersections.py index b134d9029..f348fb421 100644 --- a/test/test_intersections.py +++ b/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,7 +100,7 @@ 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)), @@ -108,14 +108,14 @@ def test_GCA_constLat_intersections_antimeridian(self): 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)), @@ -123,10 +123,10 @@ def test_GCA_constLat_intersections_empty(self): 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) diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index 16cc13194..60ddf3737 100644 --- a/uxarray/grid/integrate.py +++ b/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: diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index d57969db5..a0707310a 100644 --- a/uxarray/grid/intersections.py +++ b/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