Skip to content

Commit

Permalink
fix tests and add test case for case at pole
Browse files Browse the repository at this point in the history
  • Loading branch information
philipc2 committed Dec 16, 2024
1 parent 60cf6d6 commit f783f61
Show file tree
Hide file tree
Showing 4 changed files with 65 additions and 43 deletions.
37 changes: 18 additions & 19 deletions test/test_arcs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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):
Expand All @@ -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],
Expand All @@ -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]))



Expand All @@ -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]))
24 changes: 24 additions & 0 deletions test/test_intersections.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
18 changes: 9 additions & 9 deletions uxarray/grid/arcs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
29 changes: 14 additions & 15 deletions uxarray/grid/intersections.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down

0 comments on commit f783f61

Please sign in to comment.