diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst index 84380a678..0e0c6e287 100644 --- 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 new file mode 100644 index 000000000..c2e24f5aa --- /dev/null +++ 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 index 549ae95c5..da0fec36a 100644 --- 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 index 5f918d94a..c594a6753 100644 --- 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.