From f783f612c0dc2d6a42f93d4621789e7376327fc5 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 16 Dec 2024 14:48:10 -0600 Subject: [PATCH] fix tests and add test case for case at pole --- test/test_arcs.py | 37 +++++++++++++++++------------------ test/test_intersections.py | 24 +++++++++++++++++++++++ uxarray/grid/arcs.py | 18 ++++++++--------- uxarray/grid/intersections.py | 29 +++++++++++++-------------- 4 files changed, 65 insertions(+), 43 deletions(-) diff --git a/test/test_arcs.py b/test/test_arcs.py index 4461d2f01..4e9262982 100644 --- a/test/test_arcs.py +++ b/test/test_arcs.py @@ -10,7 +10,7 @@ import uxarray as ux from uxarray.grid.coordinates import _lonlat_rad_to_xyz -from uxarray.grid.arcs import point_within_gca, _point_within_gca_cartesian +from uxarray.grid.arcs import point_within_gca try: import constants @@ -31,29 +31,29 @@ class TestIntersectionPoint(TestCase): def test_pt_within_gcr(self): # The GCR that's eexactly 180 degrees will have Value Error raised - gcr_180degree_cart = [ + gcr_180degree_cart = np.asarray([ _lonlat_rad_to_xyz(0.0, np.pi / 2.0), _lonlat_rad_to_xyz(0.0, -np.pi / 2.0) - ] + ]) + pt_same_lon_in = np.asarray(_lonlat_rad_to_xyz(0.0, 0.0)) - pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) with self.assertRaises(ValueError): - _point_within_gca_cartesian(pt_same_lon_in, gcr_180degree_cart) - + point_within_gca(pt_same_lon_in, gcr_180degree_cart[0],gcr_180degree_cart[1] ) + # # Test when the point and the GCA all have the same longitude - gcr_same_lon_cart = [ + gcr_same_lon_cart = np.asarray([ _lonlat_rad_to_xyz(0.0, 1.5), _lonlat_rad_to_xyz(0.0, -1.5) - ] - pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) - self.assertTrue(_point_within_gca_cartesian(pt_same_lon_in, gcr_same_lon_cart)) + ]) + pt_same_lon_in = np.asarray(_lonlat_rad_to_xyz(0.0, 0.0)) + self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart[0], gcr_same_lon_cart[1])) - pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.5000001) - res = _point_within_gca_cartesian(pt_same_lon_out, gcr_same_lon_cart) + pt_same_lon_out = np.asarray(_lonlat_rad_to_xyz(0.0, 1.5000001)) + res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart[0], gcr_same_lon_cart[1]) self.assertFalse(res) - pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0) - res = _point_within_gca_cartesian(pt_same_lon_out_2, gcr_same_lon_cart) + pt_same_lon_out_2 = np.asarray(_lonlat_rad_to_xyz(0.1, 1.0)) + res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart[0], gcr_same_lon_cart[1]) self.assertFalse(res) def test_pt_within_gcr_antimeridian(self): @@ -64,13 +64,13 @@ def test_pt_within_gcr_antimeridian(self): pt_cart = np.array( [0.9438777657502077, 0.1193199333436068, 0.922714737029319]) self.assertTrue( - _point_within_gca_cartesian(pt_cart, gcr_cart)) + point_within_gca(pt_cart, gcr_cart[0], gcr_cart[1])) gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724, 0.593]]) # If we flip the gcr in the undirected mode, it should still work self.assertTrue( - _point_within_gca_cartesian(pt_cart, gcr_cart_flip)) + point_within_gca(pt_cart, gcr_cart_flip[0], gcr_cart_flip[1])) # 2nd anti-meridian case # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828], @@ -81,7 +81,7 @@ def test_pt_within_gcr_antimeridian(self): pt_cart_within = np.array( [0.6136726305712109, 0.28442243941920053, -0.365605190899831]) self.assertFalse( - _point_within_gca_cartesian(pt_cart_within, gcr_cart_1)) + point_within_gca(pt_cart_within, gcr_cart_1[0], gcr_cart_1[1])) @@ -93,5 +93,4 @@ def test_pt_within_gcr_cross_pole(self): # 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_cartesian(pt_cart, gcr_cart)) - + self.assertTrue(point_within_gca(pt_cart, gcr_cart[0], gcr_cart[1])) diff --git a/test/test_intersections.py b/test/test_intersections.py index 8983859d0..5bb542f9f 100644 --- a/test/test_intersections.py +++ b/test/test_intersections.py @@ -92,6 +92,30 @@ def test_get_GCA_GCA_intersections_perpendicular(self): self.assertTrue(len(res_cart) == 0) + def test_GCA_GCA_pole(self): + face_lonlat = np.deg2rad(np.array([-175, 26.5])) + + # this fails when the pole is set to exactly -90.0 + ref_point_lonlat = np.deg2rad(np.array([0.0, -89.9])) + face_xyz = np.array(_lonlat_rad_to_xyz(*face_lonlat)) + ref_point_xyz = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat)) + + edge_a_lonlat = np.deg2rad(np.array((-175, -24.5))) + edge_b_lonlat = np.deg2rad(np.array((-173, 25.7))) + + edge_a_xyz = np.array(_lonlat_rad_to_xyz(*edge_a_lonlat)) + edge_b_xyz = np.array(_lonlat_rad_to_xyz(*edge_b_lonlat)) + + gca_a_xyz = np.array([face_xyz, ref_point_xyz]) + + gca_b_xyz = np.array([edge_a_xyz, edge_b_xyz]) + + # The edge should intersect + self.assertTrue(len(gca_gca_intersection(gca_a_xyz, gca_b_xyz))) + + + + class TestGCAconstLatIntersection(TestCase): def test_GCA_constLat_intersections_antimeridian(self): diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 3d54ba0be..90df5f0d9 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -28,15 +28,15 @@ def _to_list(obj): return obj -def _point_within_gca_cartesian(pt_xyz, gca_xyz): - pt_xyz = np.asarray(pt_xyz) - gca_xyz = np.asarray(gca_xyz) - - gca_a_xyz = gca_xyz[0] - - gca_b_xyz = gca_xyz[1] - - return point_within_gca(pt_xyz, gca_a_xyz, gca_b_xyz) +# def _point_within_gca_cartesian(pt_xyz, gca_xyz): +# pt_xyz = np.asarray(pt_xyz) +# gca_xyz = np.asarray(gca_xyz) +# +# gca_a_xyz = gca_xyz[0] +# +# gca_b_xyz = gca_xyz[1] +# +# return point_within_gca(pt_xyz, gca_a_xyz, gca_b_xyz) @njit(cache=True) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index de865d4db..8189acc86 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -1,13 +1,14 @@ import numpy as np from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE, INT_DTYPE -from uxarray.grid.utils import _newton_raphson_solver_for_gca_constLat, _angle_of_2_vectors +from uxarray.grid.utils import ( + _newton_raphson_solver_for_gca_constLat, + _angle_of_2_vectors, +) from uxarray.grid.arcs import ( - point_within_gca, in_between, _extreme_gca_latitude_cartesian, - _point_within_gca_cartesian, + point_within_gca, ) -from uxarray.grid.coordinates import _xyz_to_lonlat_rad_scalar import platform import warnings from uxarray.utils.computing import cross_fma, allclose, cross, norm @@ -252,24 +253,22 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz): x1_xyz = cross_norms x2_xyz = -x1_xyz # Check intersection points - if point_within_gca( - x1_xyz, w0_xyz, w1_xyz - ) and point_within_gca(x1_xyz, v0_xyz, v1_xyz): + if point_within_gca(x1_xyz, w0_xyz, w1_xyz) and point_within_gca( + x1_xyz, v0_xyz, v1_xyz + ): res[count, :] = x1_xyz count += 1 - if point_within_gca( - x2_xyz, w0_xyz, w1_xyz - ) and point_within_gca(x2_xyz, v0_xyz, v1_xyz): + if point_within_gca(x2_xyz, w0_xyz, w1_xyz) and point_within_gca( + x2_xyz, v0_xyz, v1_xyz + ): res[count, :] = x2_xyz count += 1 return res[:count, :] -def gca_const_lat_intersection( - gca_cart, constZ, fma_disabled=True, verbose=False -): +def gca_const_lat_intersection(gca_cart, constZ, fma_disabled=True, verbose=False): """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a constant latitude line in a Cartesian coordinate system. @@ -360,7 +359,7 @@ def gca_const_lat_intersection( res = None # Now test which intersection point is within the GCA range - if _point_within_gca_cartesian(p1, gca_cart): + if point_within_gca(p1, gca_cart[0], gca_cart[1]): try: converged_pt = _newton_raphson_solver_for_gca_constLat( p1, gca_cart, verbose=verbose @@ -382,7 +381,7 @@ def gca_const_lat_intersection( except RuntimeError: raise RuntimeError(f"Error encountered with initial guess: {p1}") - if _point_within_gca_cartesian(p2, gca_cart): + if point_within_gca(p2, gca_cart[0], gca_cart[1]): try: converged_pt = _newton_raphson_solver_for_gca_constLat( p2, gca_cart, verbose=verbose