Point on GCA helper function

docs/internal_api/index.rst

@@ -152,3 +152,4 @@ Helpers
 
    utils.helpers._is_ugrid
    utils.helpers._replace_fill_values
+   utils.helpers._angle_of_2_vectors

docs/user_api/index.rst

@@ -155,3 +155,5 @@ Helpers
    node_lonlat_rad_to_xyz
    normalize_in_place
    parse_grid_type
+   is_between
+   point_within_GCR

test/test_helpers.py

@@ -192,3 +192,100 @@ def test_replace_fill_values_invalid(self):
                                                    original_fill=-1,
                                                    new_fill=INT_FILL_VALUE,
                                                    new_dtype=np.int16)
+
+
+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.utils.helpers.node_lonlat_rad_to_xyz([0.0, 0.0]),
+            ux.utils.helpers.node_lonlat_rad_to_xyz([np.pi, 0.0])
+        ]
+        pt_same_lon_in = ux.utils.helpers.node_lonlat_rad_to_xyz([0.0, 0.0])
+        with self.assertRaises(ValueError):
+            ux.utils.helpers.point_within_GCR(pt_same_lon_in,
+                                              gcr_180degree_cart)
+
+        gcr_180degree_cart = [
+            ux.utils.helpers.node_lonlat_rad_to_xyz([0.0, np.pi / 2.0]),
+            ux.utils.helpers.node_lonlat_rad_to_xyz([0.0, -np.pi / 2.0])
+        ]
+
+        pt_same_lon_in = ux.utils.helpers.node_lonlat_rad_to_xyz([0.0, 0.0])
+        with self.assertRaises(ValueError):
+            ux.utils.helpers.point_within_GCR(pt_same_lon_in,
+                                              gcr_180degree_cart)
+
+        # Test when the point and the GCR all have the same longitude
+        gcr_same_lon_cart = [
+            ux.utils.helpers.node_lonlat_rad_to_xyz([0.0, 1.5]),
+            ux.utils.helpers.node_lonlat_rad_to_xyz([0.0, -1.5])
+        ]
+        pt_same_lon_in = ux.utils.helpers.node_lonlat_rad_to_xyz([0.0, 0.0])
+        self.assertTrue(
+            ux.utils.helpers.point_within_GCR(pt_same_lon_in,
+                                              gcr_same_lon_cart))
+
+        pt_same_lon_out = ux.utils.helpers.node_lonlat_rad_to_xyz(
+            [0.0, 1.500000000000001])
+        res = ux.utils.helpers.point_within_GCR(pt_same_lon_out,
+                                                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.utils.helpers.node_lonlat_rad_to_xyz(
+            [0.0, 1.5000000000000001])
+        res = ux.utils.helpers.point_within_GCR(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(
+            ux.utils.helpers.point_within_GCR(pt_cart_within, gcr_cart_2))
+
+        # 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(ux.utils.helpers.point_within_GCR(pt_cart, gcr_cart))
+        # If we swap the gcr, it should still be true
+        gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
+                                                          0.593]])
+        self.assertTrue(
+            ux.utils.helpers.point_within_GCR(pt_cart, gcr_cart_flip))
+
+        # 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(
+            ux.utils.helpers.point_within_GCR(pt_cart_within, gcr_cart_1))
+
+        # The following two case should work even swapping the GCR
+        v1_rad = [0.1, 0.0]
+        v2_rad = [2 * np.pi - 0.1, 0.0]
+        v1_cart = ux.utils.helpers.node_lonlat_rad_to_xyz(v1_rad)
+        v2_cart = ux.utils.helpers.node_lonlat_rad_to_xyz(v2_rad)
+        gcr_cart = np.array([v1_cart, v2_cart])
+        pt_cart = ux.utils.helpers.node_lonlat_rad_to_xyz([0.01, 0.0])
+        self.assertTrue(ux.utils.helpers.point_within_GCR(pt_cart, gcr_cart))
+        gcr_car_flipped = np.array([v2_cart, v1_cart])
+        self.assertTrue(
+            ux.utils.helpers.point_within_GCR(pt_cart, gcr_car_flipped))

uxarray/utils/constants.py

@@ -3,3 +3,4 @@
 # numpy indexing code is written for np.intp
 INT_DTYPE = np.intp
 INT_FILL_VALUE = np.iinfo(INT_DTYPE).min
+ERROR_TOLERANCE = 1.0e-15

uxarray/utils/helpers.py

@@ -4,8 +4,9 @@
 from .get_quadratureDG import get_gauss_quadratureDG, get_tri_quadratureDG
 from numba import njit, config
 import math
+from typing import List, Union
 
-from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.utils.constants import INT_DTYPE, INT_FILL_VALUE, ERROR_TOLERANCE
 
 config.DISABLE_JIT = False
 
@@ -153,17 +154,13 @@ def calculate_face_area(x,
         node3 = np.array([x[j + 2], y[j + 2], z[j + 2]], dtype=np.float64)
 
         if (coords_type == "spherical"):
-            node1 = np.array(
-                node_lonlat_rad_to_xyz([np.deg2rad(x[0]),
-                                        np.deg2rad(y[0])]))
-            node2 = np.array(
-                node_lonlat_rad_to_xyz(
-                    [np.deg2rad(x[j + 1]),
-                     np.deg2rad(y[j + 1])]))
-            node3 = np.array(
-                node_lonlat_rad_to_xyz(
-                    [np.deg2rad(x[j + 2]),
-                     np.deg2rad(y[j + 2])]))
+            node1 = node_lonlat_rad_to_xyz([np.deg2rad(x[0]), np.deg2rad(y[0])])
+            node2 = node_lonlat_rad_to_xyz(
+                [np.deg2rad(x[j + 1]),
+                 np.deg2rad(y[j + 1])])
+            node3 = node_lonlat_rad_to_xyz(
+                [np.deg2rad(x[j + 2]),
+                 np.deg2rad(y[j + 2])])
 
         for p in range(len(dW)):
             if quadrature_rule == "gaussian":
@@ -473,25 +470,29 @@ def node_lonlat_rad_to_xyz(node_coord):
 
     Parameters
     ----------
-    node: float list
+    node: float list/np.array
         2D coordinates[longitude, latitude] in radiance
 
     Returns
     ----------
-    float list
+    float np.array
         the result array of the unit 3D coordinates [x, y, z] vector where :math:`x^2 + y^2 + z^2 = 1`
 
     Raises
     ----------
     RuntimeError
         The input array doesn't have the size of 3.
     """
-    if len(node_coord) != 2:
+    node_coord = np.asarray(node_coord)
+    if node_coord.shape[0] != 2:
         raise RuntimeError(
             "Input array should have a length of 2: [longitude, latitude]")
     lon = node_coord[0]
     lat = node_coord[1]
-    return [np.cos(lon) * np.cos(lat), np.sin(lon) * np.cos(lat), np.sin(lat)]
+    return np.array(
+        [np.cos(lon) * np.cos(lat),
+         np.sin(lon) * np.cos(lat),
+         np.sin(lat)])
 
 
 @njit
@@ -501,20 +502,21 @@ def node_xyz_to_lonlat_rad(node_coord):
 
     Parameters
     ----------
-    node_coord: float list
+    node_coord: float list/np.array
         3D Cartesian Coordinates [x, y, z] of the node
 
     Returns
     ----------
-    float list
+    float np.array
         the result array of longitude and latitude in radian [longitude_rad, latitude_rad]
 
     Raises
     ----------
     RuntimeError
         The input array doesn't have the size of 3.
     """
-    if len(node_coord) != 3:
+    node_coord = np.asarray(node_coord)
+    if node_coord.shape[0] != 3:
         raise RuntimeError("Input array should have a length of 3: [x, y, z]")
     reference_tolerance = 1.0e-12
     [dx, dy, dz] = normalize_in_place(node_coord)
@@ -539,30 +541,33 @@ def node_xyz_to_lonlat_rad(node_coord):
 
 
 @njit
-def normalize_in_place(node):
+def normalize_in_place(
+        node: Union[np.ndarray, List[Union[float]]]) -> np.ndarray:
     """Helper function to project an arbitrary node in 3D coordinates [x, y, z]
     on the unit sphere. It uses the `np.linalg.norm` internally to calculate
     the magnitude.
 
     Parameters
     ----------
-    node: float list
+    node: np.array or python list of `float` or gmpy2.mpfr
         3D Cartesian Coordinates [x, y, z]
 
     Returns
     ----------
-    float list
+    normalized_node: np.array of `float` or gmpy2.mpfr
         the result unit vector [x, y, z] where :math:`x^2 + y^2 + z^2 = 1`
 
     Raises
     ----------
     RuntimeError
         The input array doesn't have the size of 3.
     """
-    if len(node) != 3:
-        raise RuntimeError("Input array should have a length of 3: [x, y, z]")
+    node = np.asarray(node)
 
-    return np.array(node) / np.linalg.norm(np.array(node), ord=2)
+    if node.shape[0] == 3:
+        return node / np.linalg.norm(node, ord=2)
+    else:
+        raise RuntimeError("Input node should be a 3D array.")
 
 
 def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
@@ -673,3 +678,114 @@ def close_face_nodes(Mesh2_face_nodes, nMesh2_face, nMaxMesh2_face_nodes):
     np.put(closed.ravel(), first_fv_idx_1d, first_node_value)
 
     return closed
+
+
+def point_within_GCR(pt, gcr_cart):
+    """Check if a point is on a given Great Circle Arc. The anti-meridian case
+    is already considered.
+
+    Parameters
+    ----------
+    pt : np.ndarray of float
+        Cartesian coordinates of the point
+    gcr_cart : np.ndarray of shape (2, 3), (np.float or gmpy2.mpfr)
+        Cartesian coordinates of the GCR
+
+    Returns
+    -------
+    bool
+        True if the point is on the GCR (Lies between the two endpoints), False otherwise
+    """
+    # Convert the cartesian coordinates to lonlat coordinates, the node_xyz_to_lonlat_rad is already overloaded
+    # with gmpy2 data type
+    pt_lonlat = node_xyz_to_lonlat_rad(pt)
+    GCRv0_lonlat = node_xyz_to_lonlat_rad(gcr_cart[0])
+    GCRv1_lonlat = node_xyz_to_lonlat_rad(gcr_cart[1])
+
+    # First if the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes
+    angle = _angle_of_2_vectors(gcr_cart[0], gcr_cart[1])
+    if np.allclose(angle, np.pi, rtol=0, atol=ERROR_TOLERANCE):
+        raise ValueError(
+            "The input GCR is exactly 180 degree, this GCR can have multiple planes. Consider breaking the GCR "
+            "into two GCRs")
+
+    if not np.allclose(np.dot(np.cross(gcr_cart[0], gcr_cart[1]), pt),
+                       0,
+                       rtol=0,
+                       atol=ERROR_TOLERANCE):
+        return False
+
+    # Special case: If the gcr goes across the anti-meridian
+    # If the pt and the GCR are on the same longitude (the y coordinates are the same)
+    if GCRv0_lonlat[0] == GCRv1_lonlat[0] == pt_lonlat[0]:
+        # Now use the latitude to determine if the pt falls between the interval
+        return is_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1])
+
+    # 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) > π
+    if abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]) > np.pi:
+
+        # The necessary condition: the pt longitude is on the opposite side of the anti-meridian
+        # Case 1: where 0 --> x0--> 180 -->x1 -->0 case is lager than the 180degrees (pi radians)
+        # Need to redirect such that 180 --> x0 --> 0 --> x1 --> 180
+        if GCRv0_lonlat[0] <= np.pi <= GCRv1_lonlat[0]:
+            return is_between(0, pt_lonlat[0], GCRv0_lonlat[0]) or is_between(
+                GCRv1_lonlat[0], pt_lonlat[0], 2 * np.pi)
+
+        # The necessary condition: the pt longitude is on the opposite side of the anti-meridian
+        # Case 2: The anti-meridian case where 180 -->x0 --> 0 lon --> x1 --> 180 that doesn't require redirection.
+        elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0:
+            return is_between(GCRv0_lonlat[0],
+                              pt_lonlat[0], 2 * np.pi) or is_between(
+                                  0, pt_lonlat[0], GCRv1_lonlat[0])
+
+    # The non-anti-meridian case.
+    else:
+        is_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0])
+        return is_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0])
+
+
+def is_between(p, q, r) -> bool:
+    """Determines whether the number q is between p and r.
+
+    Parameters
+    ----------
+    p : float
+        The lower bound.
+    q : float
+        The number to check.
+    r : float
+        The upper bound.
+
+    Returns
+    -------
+    bool
+        True if q is between p and r, False otherwise.
+    """
+
+    return p <= q <= r or r <= q <= p
+
+
+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.
+
+    Parameters
+    ----------
+    u : numpy.array
+        The first 3D vector.
+    v : numpy.array
+        The second 3D vector.
+
+    Returns
+    -------
+    float
+        The angle between u and v in radians.
+    """
+    v_norm_times_u = np.linalg.norm(v) * u
+    u_norm_times_v = np.linalg.norm(u) * v
+    vec_minus = v_norm_times_u - u_norm_times_v
+    vec_sum = v_norm_times_u + u_norm_times_v
+    angle_u_v_rad = 2 * np.arctan2(np.linalg.norm(vec_minus),
+                                   np.linalg.norm(vec_sum))
+    return angle_u_v_rad