Merge pull request #535 from hongyuchen1030/Maximum_GCA_Latitude

Maximum Latitude of a Great Circle Arc

docs/internal_api/index.rst

@@ -151,6 +151,8 @@ Arcs
    :toctree: _autosummary
 
    grid.arcs._angle_of_2_vectors
+   grid.arcs._angle_of_2_vectors
+
 
 Utils
 -----

docs/user_api/index.rst

@@ -278,6 +278,7 @@ Arcs
 
    grid.arcs.in_between
    grid.arcs.point_within_gca
+   grid.arcs.extreme_gca_latitude
 
 Intersections
 -------------

test/test_geometry.py

@@ -7,6 +7,7 @@
 from pathlib import Path
 
 import uxarray as ux
+from uxarray.constants import ERROR_TOLERANCE
 
 from spatialpandas.geometry import MultiPolygon
 
@@ -153,3 +154,231 @@ def test_pole_point_inside_polygon_from_vertice_cross(self):
         result = ux.grid.geometry._pole_point_inside_polygon(
             'North', face_edge_cart)
         self.assertTrue(result, "North pole should be inside the polygon")
+
+
+class TestLatlonBound(TestCase):
+
+    def _max_latitude_rad_iterative(self, gca_cart):
+        """Calculate the maximum latitude of a great circle arc defined by two
+        points.
+
+        Parameters
+        ----------
+        gca_cart : numpy.ndarray
+            An array containing two 3D vectors that define a great circle arc.
+
+        Returns
+        -------
+        float
+            The maximum latitude of the great circle arc in radians.
+
+        Raises
+        ------
+        ValueError
+            If the input vectors are not valid 2-element lists or arrays.
+
+        Notes
+        -----
+        The method divides the great circle arc into subsections, iteratively refining the subsection of interest
+        until the maximum latitude is found within a specified tolerance.
+        """
+
+        # Convert input vectors to radians and Cartesian coordinates
+
+        v1_cart, v2_cart = gca_cart
+        b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v1_cart.tolist())
+        c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v2_cart.tolist())
+
+        # Initialize variables for the iterative process
+        v_temp = ux.grid.utils.cross_fma(v1_cart, v2_cart)
+        v0 = ux.grid.utils.cross_fma(v_temp, v1_cart)
+        v0 = ux.grid.coordinates.normalize_in_place(v0.tolist())
+        max_section = [v1_cart, v2_cart]
+
+        # Iteratively find the maximum latitude
+        while np.abs(b_lonlat[1] - c_lonlat[1]) >= ERROR_TOLERANCE or np.abs(
+                b_lonlat[0] - c_lonlat[0]) >= ERROR_TOLERANCE:
+            max_lat = -np.pi
+            v_b, v_c = max_section
+            angle_v1_v2_rad = ux.grid.arcs._angle_of_2_vectors(v_b, v_c)
+            v0 = ux.grid.utils.cross_fma(v_temp, v_b)
+            v0 = ux.grid.coordinates.normalize_in_place(v0.tolist())
+            avg_angle_rad = angle_v1_v2_rad / 10.0
+
+            for i in range(10):
+                angle_rad_prev = avg_angle_rad * i
+                angle_rad_next = angle_rad_prev + avg_angle_rad if i < 9 else angle_v1_v2_rad
+                w1_new = np.cos(angle_rad_prev) * v_b + np.sin(
+                    angle_rad_prev) * np.array(v0)
+                w2_new = np.cos(angle_rad_next) * v_b + np.sin(
+                    angle_rad_next) * np.array(v0)
+                w1_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
+                    w1_new.tolist())
+                w2_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
+                    w2_new.tolist())
+
+                # Adjust latitude boundaries to avoid error accumulation
+                if i == 0:
+                    w1_lonlat[1] = b_lonlat[1]
+                elif i >= 9:
+                    w2_lonlat[1] = c_lonlat[1]
+
+                # Update maximum latitude and section if needed
+                max_lat = max(max_lat, w1_lonlat[1], w2_lonlat[1])
+                if np.abs(w2_lonlat[1] -
+                          w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[
+                              1] == max_lat == w2_lonlat[1]:
+                    max_section = [w1_new, w2_new]
+                    break
+                if np.abs(max_lat - w1_lonlat[1]) <= ERROR_TOLERANCE:
+                    max_section = [w1_new, w2_new] if i != 0 else [v_b, w2_new]
+                elif np.abs(max_lat - w2_lonlat[1]) <= ERROR_TOLERANCE:
+                    max_section = [w1_new, w2_new] if i != 9 else [w1_new, v_c]
+
+            # Update longitude and latitude for the next iteration
+            b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
+                max_section[0].tolist())
+            c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
+                max_section[1].tolist())
+
+        return np.average([b_lonlat[1], c_lonlat[1]])
+
+    def _min_latitude_rad_iterative(self, gca_cart):
+        """Calculate the minimum latitude of a great circle arc defined by two
+        points.
+
+        Parameters
+        ----------
+        gca_cart : numpy.ndarray
+            An array containing two 3D vectors that define a great circle arc.
+
+        Returns
+        -------
+        float
+            The minimum latitude of the great circle arc in radians.
+
+        Raises
+        ------
+        ValueError
+            If the input vectors are not valid 2-element lists or arrays.
+
+        Notes
+        -----
+        The method divides the great circle arc into subsections, iteratively refining the subsection of interest
+        until the minimum latitude is found within a specified tolerance.
+        """
+
+        # Convert input vectors to radians and Cartesian coordinates
+        v1_cart, v2_cart = gca_cart
+        b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v1_cart.tolist())
+        c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v2_cart.tolist())
+
+        # Initialize variables for the iterative process
+        v_temp = ux.grid.utils.cross_fma(v1_cart, v2_cart)
+        v0 = ux.grid.utils.cross_fma(v_temp, v1_cart)
+        v0 = np.array(ux.grid.coordinates.normalize_in_place(v0.tolist()))
+        min_section = [v1_cart, v2_cart]
+
+        # Iteratively find the minimum latitude
+        while np.abs(b_lonlat[1] - c_lonlat[1]) >= ERROR_TOLERANCE or np.abs(
+                b_lonlat[0] - c_lonlat[0]) >= ERROR_TOLERANCE:
+            min_lat = np.pi
+            v_b, v_c = min_section
+            angle_v1_v2_rad = ux.grid.arcs._angle_of_2_vectors(v_b, v_c)
+            v0 = ux.grid.utils.cross_fma(v_temp, v_b)
+            v0 = np.array(ux.grid.coordinates.normalize_in_place(v0.tolist()))
+            avg_angle_rad = angle_v1_v2_rad / 10.0
+
+            for i in range(10):
+                angle_rad_prev = avg_angle_rad * i
+                angle_rad_next = angle_rad_prev + avg_angle_rad if i < 9 else angle_v1_v2_rad
+                w1_new = np.cos(angle_rad_prev) * v_b + np.sin(
+                    angle_rad_prev) * v0
+                w2_new = np.cos(angle_rad_next) * v_b + np.sin(
+                    angle_rad_next) * v0
+                w1_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
+                    w1_new.tolist())
+                w2_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
+                    w2_new.tolist())
+
+                # Adjust latitude boundaries to avoid error accumulation
+                if i == 0:
+                    w1_lonlat[1] = b_lonlat[1]
+                elif i >= 9:
+                    w2_lonlat[1] = c_lonlat[1]
+
+                # Update minimum latitude and section if needed
+                min_lat = min(min_lat, w1_lonlat[1], w2_lonlat[1])
+                if np.abs(w2_lonlat[1] -
+                          w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[
+                              1] == min_lat == w2_lonlat[1]:
+                    min_section = [w1_new, w2_new]
+                    break
+                if np.abs(min_lat - w1_lonlat[1]) <= ERROR_TOLERANCE:
+                    min_section = [w1_new, w2_new] if i != 0 else [v_b, w2_new]
+                elif np.abs(min_lat - w2_lonlat[1]) <= ERROR_TOLERANCE:
+                    min_section = [w1_new, w2_new] if i != 9 else [w1_new, v_c]
+
+            # Update longitude and latitude for the next iteration
+            b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
+                min_section[0].tolist())
+            c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
+                min_section[1].tolist())
+
+        return np.average([b_lonlat[1], c_lonlat[1]])
+
+    def test_extreme_gca_latitude_max(self):
+        # Define a great circle arc that is symmetrical around 0 degrees longitude
+        gca_cart = np.array([
+            ux.grid.coordinates.normalize_in_place([0.5, 0.5, 0.5]),
+            ux.grid.coordinates.normalize_in_place([-0.5, 0.5, 0.5])
+        ])
+
+        # Calculate the maximum latitude
+        max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max')
+
+        # Check if the maximum latitude is correct
+        expected_max_latitude = self._max_latitude_rad_iterative(gca_cart)
+        self.assertAlmostEqual(max_latitude,
+                               expected_max_latitude,
+                               delta=ERROR_TOLERANCE)
+
+        # Define a great circle arc in 3D space
+        gca_cart = np.array([[0.0, 0.0, 1.0], [1.0, 0.0, 0.0]])
+
+        # Calculate the maximum latitude
+        max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max')
+
+        # Check if the maximum latitude is correct
+        expected_max_latitude = np.pi / 2  # 90 degrees in radians
+        self.assertAlmostEqual(max_latitude,
+                               expected_max_latitude,
+                               delta=ERROR_TOLERANCE)
+
+    def test_extreme_gca_latitude_min(self):
+        # Define a great circle arc that is symmetrical around 0 degrees longitude
+        gca_cart = np.array([
+            ux.grid.coordinates.normalize_in_place([0.5, 0.5, -0.5]),
+            ux.grid.coordinates.normalize_in_place([-0.5, 0.5, -0.5])
+        ])
+
+        # Calculate the minimum latitude
+        min_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'min')
+
+        # Check if the minimum latitude is correct
+        expected_min_latitude = self._min_latitude_rad_iterative(gca_cart)
+        self.assertAlmostEqual(min_latitude,
+                               expected_min_latitude,
+                               delta=ERROR_TOLERANCE)
+
+        # Define a great circle arc in 3D space
+        gca_cart = np.array([[0.0, 0.0, -1.0], [1.0, 0.0, 0.0]])
+
+        # Calculate the minimum latitude
+        min_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'min')
+
+        # Check if the minimum latitude is correct
+        expected_min_latitude = -np.pi / 2  # 90 degrees in radians
+        self.assertAlmostEqual(min_latitude,
+                               expected_min_latitude,
+                               delta=ERROR_TOLERANCE)

uxarray/grid/arcs.py

@@ -1,7 +1,7 @@
 import numpy as np
 
-from uxarray.grid.coordinates import node_xyz_to_lonlat_rad
-
+from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place
+from enum import Enum
 from uxarray.constants import ERROR_TOLERANCE
 
 
@@ -150,3 +150,52 @@ def _angle_of_2_vectors(u, v):
     angle_u_v_rad = 2 * np.arctan2(np.linalg.norm(vec_minus),
                                    np.linalg.norm(vec_sum))
     return angle_u_v_rad
+
+
+def extreme_gca_latitude(gca_cart, extreme_type):
+    """Calculate the maximum or minimum latitude of a great circle arc defined
+    by two 3D points.
+
+    Parameters
+    ----------
+    gca_cart : numpy.ndarray
+        An array containing two 3D vectors that define a great circle arc.
+
+    extreme_type : str
+        The type of extreme latitude to calculate. Must be either 'max' or 'min'.
+
+    Returns
+    -------
+    float
+        The maximum or minimum latitude of the great circle arc in radians.
+
+    Raises
+    ------
+    ValueError
+        If `extreme_type` is not 'max' or 'min'.
+    """
+    extreme_type = extreme_type.lower()
+
+    if extreme_type not in ('max', 'min'):
+        raise ValueError("extreme_type must be either 'max' or 'min'")
+
+    n1, n2 = gca_cart
+    dot_n1_n2 = np.dot(n1, n2)
+    denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0)
+    d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom
+
+    d_a_max = np.clip(d_a_max, 0, 1) if np.isclose(
+        d_a_max, [0, 1], atol=ERROR_TOLERANCE).any() else d_a_max
+    lat_n1, lat_n2 = node_xyz_to_lonlat_rad(
+        n1.tolist())[1], node_xyz_to_lonlat_rad(n2.tolist())[1]
+
+    if 0 < d_a_max < 1:
+        node3 = (1 - d_a_max) * n1 + d_a_max * n2
+        node3 = np.array(normalize_in_place(node3.tolist()))
+        d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1))
+
+        return max(d_lat_rad, lat_n1, lat_n2) if extreme_type == 'max' else min(
+            d_lat_rad, lat_n1, lat_n2)
+    else:
+        return max(lat_n1, lat_n2) if extreme_type == 'max' else min(
+            lat_n1, lat_n2)