Bugfix pt_on_gca (#698)

* Initial Draft Algorithm * Revert "Initial Draft Algorithm" This reverts commit 99573b8. * Initial commit * Revert "Initial commit" This reverts commit 1ba6880. * Initial commit * Fix complete --------- Co-authored-by: Philip Chmielowiec <[email protected]>
UXARRAY · Feb 13, 2024 · d3a4c5f · d3a4c5f
1 parent 871e270
commit d3a4c5f
Show file tree

Hide file tree

Showing 4 changed files with 229 additions and 103 deletions.
diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
@@ -161,7 +161,7 @@ Arcs
    :toctree: generated/
 
    grid.arcs._angle_of_2_vectors
-   grid.arcs._angle_of_2_vectors
+   grid.arcs._decide_pole_latitude
 
 
 Utils

diff --git a/test/test_arcs.py b/test/test_arcs.py
@@ -0,0 +1,149 @@
+import os
+import numpy as np
+import numpy.testing as nt
+import random
+import xarray as xr
+
+from unittest import TestCase
+from pathlib import Path
+
+import uxarray as ux
+
+from uxarray.grid.coordinates import node_lonlat_rad_to_xyz
+from uxarray.grid.arcs import point_within_gca
+
+try:
+    import constants
+except ImportError:
+    from . import constants
+
+# Data files
+current_path = Path(os.path.dirname(os.path.realpath(__file__)))
+
+gridfile_exo_CSne8 = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g"
+gridfile_scrip_CSne8 = current_path / 'meshfiles' / "scrip" / "outCSne8" / 'outCSne8.nc'
+gridfile_geoflowsmall_grid = current_path / 'meshfiles' / "ugrid" / "geoflow-small" / 'grid.nc'
+gridfile_geoflowsmall_var = current_path / 'meshfiles' / "ugrid" / "geoflow-small" / 'v1.nc'
+
+
+class TestIntersectionPoint(TestCase):
+
+    def test_pt_within_gcr(self):
+        # The GCR that's eexactly 180 degrees will have Value Error raised
+        gcr_180degree_cart = [
+            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0]),
+            ux.grid.coordinates.node_lonlat_rad_to_xyz([np.pi, 0.0])
+        ]
+        pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0])
+        with self.assertRaises(ValueError):
+            point_within_gca(pt_same_lon_in, gcr_180degree_cart)
+
+        gcr_180degree_cart = [
+            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, np.pi / 2.0]),
+            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, -np.pi / 2.0])
+        ]
+
+        pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0])
+        with self.assertRaises(ValueError):
+            point_within_gca(pt_same_lon_in, gcr_180degree_cart)
+
+        # Test when the point and the GCA all have the same longitude
+        gcr_same_lon_cart = [
+            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 1.5]),
+            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, -1.5])
+        ]
+        pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0])
+        self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart))
+
+        pt_same_lon_out = ux.grid.coordinates.node_lonlat_rad_to_xyz(
+            [0.0, 1.500000000000001])
+        res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart)
+        self.assertFalse(res)
+
+        pt_same_lon_out_2 = ux.grid.coordinates.node_lonlat_rad_to_xyz(
+            [0.1, 1.0])
+        res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart)
+        self.assertFalse(res)
+
+        # And if we increase the digital place by one, it should be true again
+        pt_same_lon_out_add_one_place = ux.grid.coordinates.node_lonlat_rad_to_xyz(
+            [0.0, 1.5000000000000001])
+        res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
+        self.assertTrue(res)
+
+        # Normal case
+        # GCR vertex0 in radian : [1.3003315590159483, -0.007004587172323237],
+        # GCR vertex1 in radian : [3.5997458123873827, -1.4893379576608758]
+        # Point in radian : [1.3005410084914981, -0.010444274637648326]
+        gcr_cart_2 = np.array([[0.267, 0.963, -0.007], [-0.073, -0.036,
+                                                        -0.997]])
+        pt_cart_within = np.array(
+            [0.25616109352676675, 0.9246590335292105, -0.010021496695000144])
+        self.assertTrue(point_within_gca(pt_cart_within, gcr_cart_2, True))
+
+        # Test other more complicate cases : The anti-meridian case
+
+    def test_pt_within_gcr_antimeridian(self):
+        # GCR vertex0 in radian : [5.163808182822441, 0.6351384888657234],
+        # GCR vertex1 in radian : [0.8280410325693055, 0.42237025187091526]
+        # Point in radian : [0.12574759138415173, 0.770098701904903]
+        gcr_cart = np.array([[0.351, -0.724, 0.593], [0.617, 0.672, 0.410]])
+        pt_cart = np.array(
+            [0.9438777657502077, 0.1193199333436068, 0.922714737029319])
+        self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=True))
+        # If we swap the gcr, it should throw a value error since it's larger than 180 degree
+        gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
+                                                          0.593]])
+        with self.assertRaises(ValueError):
+            point_within_gca(pt_cart, gcr_cart_flip, is_directed=True)
+
+        # If we flip the gcr in the undirected mode, it should still work
+        self.assertTrue(
+            point_within_gca(pt_cart, gcr_cart_flip, is_directed=False))
+
+        # 2nd anti-meridian case
+        # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
+        # GCR vertex1 in radian : [2.4269979227622533, -0.007003212877856825]
+        # Point in radian : [0.43400375562899113, -0.49554509841586936]
+        gcr_cart_1 = np.array([[-0.491, -0.706, 0.510], [-0.755, 0.655,
+                                                         -0.007]])
+        pt_cart_within = np.array(
+            [0.6136726305712109, 0.28442243941920053, -0.365605190899831])
+        self.assertFalse(
+            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=True))
+        self.assertFalse(
+            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=False))
+
+        # The first case should not work and the second should work
+        v1_rad = [0.1, 0.0]
+        v2_rad = [2 * np.pi - 0.1, 0.0]
+        v1_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz(v1_rad)
+        v2_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz(v2_rad)
+        gcr_cart = np.array([v1_cart, v2_cart])
+        pt_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.01, 0.0])
+        with self.assertRaises(ValueError):
+            point_within_gca(pt_cart, gcr_cart, is_directed=True)
+        gcr_car_flipped = np.array([v2_cart, v1_cart])
+        self.assertTrue(
+            point_within_gca(pt_cart, gcr_car_flipped, is_directed=True))
+
+    def test_pt_within_gcr_cross_pole(self):
+        gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, 0.3]])
+        pt_cart = np.array(
+            [0.10, 0.0, 0.8])
+
+        # Normalize the point abd the GCA
+        pt_cart = pt_cart / np.linalg.norm(pt_cart)
+        gcr_cart = np.array([x / np.linalg.norm(x) for x in gcr_cart])
+        self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=False))
+
+        gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, -0.6]])
+        pt_cart = np.array(
+            [0.10, 0.0, 0.8])
+
+        # When the point is not within the GCA
+        pt_cart = pt_cart / np.linalg.norm(pt_cart)
+        gcr_cart = np.array([x / np.linalg.norm(x) for x in gcr_cart])
+        self.assertFalse(point_within_gca(pt_cart, gcr_cart, is_directed=False))
+        with self.assertRaises(ValueError):
+            point_within_gca(pt_cart, gcr_cart, is_directed=True)
diff --git a/test/test_helpers.py b/test/test_helpers.py
@@ -220,107 +220,6 @@ def test_convert_face_node_conn_to_sparse_matrix(self):
         nt.assert_array_equal(nodes_indices, expected_nodes_indices)
 
 
-class TestIntersectionPoint(TestCase):
-
-    def test_pt_within_gcr(self):
-        # The GCR that's eexactly 180 degrees will have Value Error raised
-        gcr_180degree_cart = [
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0]),
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([np.pi, 0.0])
-        ]
-        pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0])
-        with self.assertRaises(ValueError):
-            point_within_gca(pt_same_lon_in, gcr_180degree_cart)
-
-        gcr_180degree_cart = [
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, np.pi / 2.0]),
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, -np.pi / 2.0])
-        ]
-
-        pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0])
-        with self.assertRaises(ValueError):
-            point_within_gca(pt_same_lon_in, gcr_180degree_cart)
-
-        # Test when the point and the GCA all have the same longitude
-        gcr_same_lon_cart = [
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 1.5]),
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, -1.5])
-        ]
-        pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0])
-        self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart))
-
-        pt_same_lon_out = ux.grid.coordinates.node_lonlat_rad_to_xyz(
-            [0.0, 1.500000000000001])
-        res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart)
-        self.assertFalse(res)
-
-        pt_same_lon_out_2 = ux.grid.coordinates.node_lonlat_rad_to_xyz(
-            [0.1, 1.0])
-        res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart)
-        self.assertFalse(res)
-
-        # And if we increase the digital place by one, it should be true again
-        pt_same_lon_out_add_one_place = ux.grid.coordinates.node_lonlat_rad_to_xyz(
-            [0.0, 1.5000000000000001])
-        res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
-        self.assertTrue(res)
-
-        # Normal case
-        # GCR vertex0 in radian : [1.3003315590159483, -0.007004587172323237],
-        # GCR vertex1 in radian : [3.5997458123873827, -1.4893379576608758]
-        # Point in radian : [1.3005410084914981, -0.010444274637648326]
-        gcr_cart_2 = np.array([[0.267, 0.963, -0.007], [-0.073, -0.036,
-                                                        -0.997]])
-        pt_cart_within = np.array(
-            [0.25616109352676675, 0.9246590335292105, -0.010021496695000144])
-        self.assertTrue(point_within_gca(pt_cart_within, gcr_cart_2, True))
-
-        # Test other more complicate cases : The anti-meridian case
-
-        # GCR vertex0 in radian : [5.163808182822441, 0.6351384888657234],
-        # GCR vertex1 in radian : [0.8280410325693055, 0.42237025187091526]
-        # Point in radian : [0.12574759138415173, 0.770098701904903]
-        gcr_cart = np.array([[0.351, -0.724, 0.593], [0.617, 0.672, 0.410]])
-        pt_cart = np.array(
-            [0.9438777657502077, 0.1193199333436068, 0.922714737029319])
-        self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=True))
-        # If we swap the gcr, it should throw a value error since it's larger than 180 degree
-        gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
-                                                          0.593]])
-        with self.assertRaises(ValueError):
-            point_within_gca(pt_cart, gcr_cart_flip, is_directed=True)
-
-        # If we flip the gcr in the undirected mode, it should still work
-        self.assertTrue(
-            point_within_gca(pt_cart, gcr_cart_flip, is_directed=False))
-
-        # 2nd anti-meridian case
-        # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
-        # GCR vertex1 in radian : [2.4269979227622533, -0.007003212877856825]
-        # Point in radian : [0.43400375562899113, -0.49554509841586936]
-        gcr_cart_1 = np.array([[-0.491, -0.706, 0.510], [-0.755, 0.655,
-                                                         -0.007]])
-        pt_cart_within = np.array(
-            [0.6136726305712109, 0.28442243941920053, -0.365605190899831])
-        self.assertFalse(
-            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=True))
-        self.assertFalse(
-            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=False))
-
-        # The first case should not work and the second should work
-        v1_rad = [0.1, 0.0]
-        v2_rad = [2 * np.pi - 0.1, 0.0]
-        v1_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz(v1_rad)
-        v2_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz(v2_rad)
-        gcr_cart = np.array([v1_cart, v2_cart])
-        pt_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.01, 0.0])
-        with self.assertRaises(ValueError):
-            point_within_gca(pt_cart, gcr_cart, is_directed=True)
-        gcr_car_flipped = np.array([v2_cart, v1_cart])
-        self.assertTrue(
-            point_within_gca(pt_cart, gcr_car_flipped, is_directed=True))
-
-
 class TestOperators(TestCase):
 
     def test_in_between(self):

diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
@@ -28,6 +28,7 @@ def point_within_gca(pt, gca_cart, is_directed=False):
     is_directed : bool, optional, default = False
         If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected,
         and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False.
+        For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1
 
     Returns
     -------
@@ -79,7 +80,7 @@ def point_within_gca(pt, gca_cart, is_directed=False):
     ):
         return False
 
-    if GCRv0_lonlat[0] == GCRv1_lonlat[0]:
+    if np.isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0, atol=ERROR_TOLERANCE):
         # If the pt and the GCA are on the same longitude (the y coordinates are the same)
         if GCRv0_lonlat[0] == pt_lonlat[0]:
             # Now use the latitude to determine if the pt falls between the interval
@@ -88,6 +89,44 @@ def point_within_gca(pt, gca_cart, is_directed=False):
             # If the pt and the GCA are not on the same longitude when the GCA is a longnitude arc, then the pt is not on the GCA
             return False
 
+    # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point
+    if np.isclose(
+        abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0, atol=ERROR_TOLERANCE
+    ):
+        if not np.isclose(
+            GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE
+        ) and not np.isclose(
+            GCRv1_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE
+        ):
+            return False
+        else:
+            # Determine the pole latitude and latitude extension
+            if (GCRv0_lonlat[1] > 0 and GCRv1_lonlat[1] > 0) or (
+                GCRv0_lonlat[1] < 0 and GCRv1_lonlat[1] < 0
+            ):
+                pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0 else -np.pi / 2
+                lat_extend = (
+                    abs(np.pi / 2 - abs(GCRv0_lonlat[1]))
+                    + np.pi / 2
+                    + abs(GCRv1_lonlat[1])
+                )
+            else:
+                pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1])
+                lat_extend = (
+                    abs(np.pi / 2 - abs(GCRv0_lonlat[1]))
+                    + np.pi / 2
+                    + abs(GCRv1_lonlat[1])
+                )
+
+            if is_directed and lat_extend >= np.pi:
+                raise ValueError(
+                    "The input GCA spans more than 180 degrees. Please divide it into two."
+                )
+
+            return in_between(GCRv0_lonlat[1], pt_lonlat[1], pole_lat) or in_between(
+                pole_lat, pt_lonlat[1], GCRv1_lonlat[1]
+            )
+
     if is_directed:
         # The anti-meridian case Sufficient condition: absolute difference between the longitudes of the two
         # vertices is greater than 180 degrees (π radians): abs(GCRv1_lon - GCRv0_lon) > π
@@ -142,6 +181,45 @@ def in_between(p, q, r) -> bool:
     return p <= q <= r or r <= q <= p
 
 
+def _decide_pole_latitude(lat1, lat2):
+    """Determine the pole latitude based on the latitudes of two points on a
+    Great Circle Arc (GCA).
+
+    This function calculates the combined latitude span from each point to its nearest pole
+    and decides which pole (North or South) the smaller GCA will pass. This decision is crucial
+    for handling GCAs that span exactly or more than 180 degrees in longitude, indicating
+    the arc might pass through or close to one of the Earth's poles.
+
+    Parameters
+    ----------
+    lat1 : float
+        Latitude of the first point in radians. Positive for the Northern Hemisphere, negative for the Southern.
+    lat2 : float
+        Latitude of the second point in radians. Positive for the Northern Hemisphere, negative for the Southern.
+
+    Returns
+    -------
+    float
+        The latitude of the pole (np.pi/2 for the North Pole or -np.pi/2 for the South Pole) the GCA is closer to.
+
+    Notes
+    -----
+    The function assumes the input latitudes are valid (i.e., between -np.pi/2 and np.pi/2) and expressed in radians.
+    The determination of which pole a GCA is closer to is based on the sum of the latitudinal spans from each point
+    to its nearest pole, considering the shortest path on the sphere.
+    """
+    # Calculate the total latitudinal span to the nearest poles
+    lat_extend = abs(np.pi / 2 - abs(lat1)) + np.pi / 2 + abs(lat2)
+
+    # Determine the closest pole based on the latitudinal span
+    if lat_extend < np.pi:
+        closest_pole = np.pi / 2 if lat1 > 0 else -np.pi / 2
+    else:
+        closest_pole = -np.pi / 2 if lat1 > 0 else np.pi / 2
+
+    return closest_pole
+
+
 def _angle_of_2_vectors(u, v):
     """Calculate the angle between two 3D vectors u and v in radians. Can be
     used to calcualte the span of a GCR.