diff --git a/docs/api.rst b/docs/api.rst index 1b03b5d3b..f9146e69f 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -126,6 +126,8 @@ Descriptors Grid.descriptors Grid.face_areas Grid.bounds + Grid.face_bounds_lon + Grid.face_bounds_lat Grid.edge_node_distances Grid.edge_face_distances Grid.antimeridian_face_indices diff --git a/test/test_arcs.py b/test/test_arcs.py index fa6bc8499..c139e3907 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 +from uxarray.grid.arcs import point_within_gca, _point_within_gca_cartesian try: import constants @@ -38,7 +38,7 @@ def test_pt_within_gcr(self): pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) with self.assertRaises(ValueError): - point_within_gca(pt_same_lon_in, gcr_180degree_cart) + _point_within_gca_cartesian(pt_same_lon_in, gcr_180degree_cart) # Test when the point and the GCA all have the same longitude gcr_same_lon_cart = [ @@ -46,19 +46,19 @@ def test_pt_within_gcr(self): _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(pt_same_lon_in, gcr_same_lon_cart)) + self.assertTrue(_point_within_gca_cartesian(pt_same_lon_in, gcr_same_lon_cart)) pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001) - res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart) + res = _point_within_gca_cartesian(pt_same_lon_out, gcr_same_lon_cart) self.assertFalse(res) pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0) - res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart) + res = _point_within_gca_cartesian(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 = _lonlat_rad_to_xyz(0.0, 1.5000000000000001) - res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart) + res = _point_within_gca_cartesian(pt_same_lon_out_add_one_place, gcr_same_lon_cart) self.assertTrue(res) # Normal case @@ -69,7 +69,7 @@ def test_pt_within_gcr(self): -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)) + self.assertTrue(_point_within_gca_cartesian(pt_cart_within, gcr_cart_2, True)) # Test other more complicate cases : The anti-meridian case @@ -80,16 +80,16 @@ def test_pt_within_gcr_antimeridian(self): 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)) + self.assertTrue(_point_within_gca_cartesian(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) + _point_within_gca_cartesian(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)) + _point_within_gca_cartesian(pt_cart, gcr_cart_flip, is_directed=False)) # 2nd anti-meridian case # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828], @@ -100,9 +100,9 @@ def test_pt_within_gcr_antimeridian(self): 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)) + _point_within_gca_cartesian(pt_cart_within, gcr_cart_1, is_directed=True)) self.assertFalse( - point_within_gca(pt_cart_within, gcr_cart_1, is_directed=False)) + _point_within_gca_cartesian(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] @@ -112,10 +112,10 @@ def test_pt_within_gcr_antimeridian(self): gcr_cart = np.array([v1_cart, v2_cart]) pt_cart = _lonlat_rad_to_xyz(0.01, 0.0) with self.assertRaises(ValueError): - point_within_gca(pt_cart, gcr_cart, is_directed=True) + _point_within_gca_cartesian(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)) + _point_within_gca_cartesian(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]]) @@ -125,7 +125,7 @@ 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(pt_cart, gcr_cart, is_directed=False)) + self.assertTrue(_point_within_gca_cartesian(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( @@ -134,6 +134,6 @@ def test_pt_within_gcr_cross_pole(self): # 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)) + self.assertFalse(_point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=False)) with self.assertRaises(ValueError): - point_within_gca(pt_cart, gcr_cart, is_directed=True) + _point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=True) diff --git a/test/test_cross_sections.py b/test/test_cross_sections.py index 1847e69b3..6fbe8abcf 100644 --- a/test/test_cross_sections.py +++ b/test/test_cross_sections.py @@ -35,38 +35,34 @@ class TestQuadHex: All four faces intersect a constant latitude of 0.0 """ - @pytest.mark.parametrize("use_spherical_bounding_box", [True, False]) - def test_constant_lat_cross_section_grid(self, use_spherical_bounding_box): - - + def test_constant_lat_cross_section_grid(self): uxgrid = ux.open_grid(quad_hex_grid_path) - grid_top_two = uxgrid.cross_section.constant_latitude(lat=0.1, use_spherical_bounding_box=use_spherical_bounding_box) + grid_top_two = uxgrid.cross_section.constant_latitude(lat=0.1, ) assert grid_top_two.n_face == 2 - grid_bottom_two = uxgrid.cross_section.constant_latitude(lat=-0.1, use_spherical_bounding_box=use_spherical_bounding_box) + grid_bottom_two = uxgrid.cross_section.constant_latitude(lat=-0.1, ) assert grid_bottom_two.n_face == 2 - grid_all_four = uxgrid.cross_section.constant_latitude(lat=0.0, use_spherical_bounding_box=use_spherical_bounding_box) + grid_all_four = uxgrid.cross_section.constant_latitude(lat=0.0, ) assert grid_all_four.n_face == 4 with pytest.raises(ValueError): # no intersections found at this line - uxgrid.cross_section.constant_latitude(lat=10.0, use_spherical_bounding_box=use_spherical_bounding_box) + uxgrid.cross_section.constant_latitude(lat=10.0, ) - @pytest.mark.parametrize("use_spherical_bounding_box", [False]) - def test_constant_lon_cross_section_grid(self, use_spherical_bounding_box): + def test_constant_lon_cross_section_grid(self): uxgrid = ux.open_grid(quad_hex_grid_path) - grid_left_two = uxgrid.cross_section.constant_longitude(lon=-0.1, use_spherical_bounding_box=use_spherical_bounding_box) + grid_left_two = uxgrid.cross_section.constant_longitude(lon=-0.1, ) assert grid_left_two.n_face == 2 - grid_right_two = uxgrid.cross_section.constant_longitude(lon=0.2, use_spherical_bounding_box=use_spherical_bounding_box) + grid_right_two = uxgrid.cross_section.constant_longitude(lon=0.2, ) assert grid_right_two.n_face == 2 @@ -74,64 +70,63 @@ def test_constant_lon_cross_section_grid(self, use_spherical_bounding_box): # no intersections found at this line uxgrid.cross_section.constant_longitude(lon=10.0) - @pytest.mark.parametrize("use_spherical_bounding_box", [False]) - def test_constant_lat_cross_section_uxds(self, use_spherical_bounding_box): + def test_constant_lat_cross_section_uxds(self): uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path) uxds.uxgrid.normalize_cartesian_coordinates() - da_top_two = uxds['t2m'].cross_section.constant_latitude(lat=0.1, use_spherical_bounding_box=use_spherical_bounding_box) + da_top_two = uxds['t2m'].cross_section.constant_latitude(lat=0.1, ) nt.assert_array_equal(da_top_two.data, uxds['t2m'].isel(n_face=[1, 2]).data) - da_bottom_two = uxds['t2m'].cross_section.constant_latitude(lat=-0.1, use_spherical_bounding_box=use_spherical_bounding_box) + da_bottom_two = uxds['t2m'].cross_section.constant_latitude(lat=-0.1, ) nt.assert_array_equal(da_bottom_two.data, uxds['t2m'].isel(n_face=[0, 3]).data) - da_all_four = uxds['t2m'].cross_section.constant_latitude(lat=0.0, use_spherical_bounding_box=use_spherical_bounding_box) + da_all_four = uxds['t2m'].cross_section.constant_latitude(lat=0.0, ) nt.assert_array_equal(da_all_four.data , uxds['t2m'].data) with pytest.raises(ValueError): # no intersections found at this line - uxds['t2m'].cross_section.constant_latitude(lat=10.0, use_spherical_bounding_box=use_spherical_bounding_box) + uxds['t2m'].cross_section.constant_latitude(lat=10.0, ) + - @pytest.mark.parametrize("use_spherical_bounding_box", [False]) - def test_constant_lon_cross_section_uxds(self, use_spherical_bounding_box): + def test_constant_lon_cross_section_uxds(self): uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path) uxds.uxgrid.normalize_cartesian_coordinates() - da_left_two = uxds['t2m'].cross_section.constant_longitude(lon=-0.1, use_spherical_bounding_box=use_spherical_bounding_box) + da_left_two = uxds['t2m'].cross_section.constant_longitude(lon=-0.1, ) nt.assert_array_equal(da_left_two.data, uxds['t2m'].isel(n_face=[0, 2]).data) - da_right_two = uxds['t2m'].cross_section.constant_longitude(lon=0.2, use_spherical_bounding_box=use_spherical_bounding_box) + da_right_two = uxds['t2m'].cross_section.constant_longitude(lon=0.2, ) nt.assert_array_equal(da_right_two.data, uxds['t2m'].isel(n_face=[1, 3]).data) with pytest.raises(ValueError): # no intersections found at this line - uxds['t2m'].cross_section.constant_longitude(lon=10.0, use_spherical_bounding_box=use_spherical_bounding_box) + uxds['t2m'].cross_section.constant_longitude(lon=10.0, ) class TestCubeSphere: - @pytest.mark.parametrize("use_spherical_bounding_box", [True, False]) - def test_north_pole(self, use_spherical_bounding_box): + + def test_north_pole(self): uxgrid = ux.open_grid(cube_sphere_grid) lats = [89.85, 89.9, 89.95, 89.99] for lat in lats: - cross_grid = uxgrid.cross_section.constant_latitude(lat=lat, use_spherical_bounding_box=use_spherical_bounding_box) + cross_grid = uxgrid.cross_section.constant_latitude(lat=lat, ) # Cube sphere grid should have 4 faces centered around the pole assert cross_grid.n_face == 4 - @pytest.mark.parametrize("use_spherical_bounding_box", [True, False]) - def test_south_pole(self, use_spherical_bounding_box): + + def test_south_pole(self): uxgrid = ux.open_grid(cube_sphere_grid) lats = [-89.85, -89.9, -89.95, -89.99] for lat in lats: - cross_grid = uxgrid.cross_section.constant_latitude(lat=lat, use_spherical_bounding_box=use_spherical_bounding_box) + cross_grid = uxgrid.cross_section.constant_latitude(lat=lat, ) # Cube sphere grid should have 4 faces centered around the pole assert cross_grid.n_face == 4 @@ -152,8 +147,7 @@ def test_constant_lat(self): candidate_faces = constant_lat_intersections_face_bounds( lat=const_lat, - face_min_lat_rad=bounds_rad[:, 0, 0], - face_max_lat_rad=bounds_rad[:, 0, 1], + face_bounds_lat=bounds_rad[:, 0], ) # Expected output @@ -176,8 +170,7 @@ def test_constant_lat_out_of_bounds(self): candidate_faces = constant_lat_intersections_face_bounds( lat=const_lat, - face_min_lat_rad=bounds_rad[:, 0, 0], - face_max_lat_rad=bounds_rad[:, 0, 1], + face_bounds_lat=bounds_rad[:, 0], ) assert len(candidate_faces) == 0 diff --git a/test/test_geometry.py b/test/test_geometry.py index 67e9369b9..05e0930fc 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -10,12 +10,11 @@ from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE import uxarray.utils.computing as ac_utils from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad -from uxarray.grid.arcs import extreme_gca_latitude +from uxarray.grid.arcs import extreme_gca_latitude, _extreme_gca_latitude_cartesian from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes -from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, stereographic_projection, \ - inverse_stereographic_projection +from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, _pole_point_inside_polygon_cartesian, stereographic_projection, inverse_stereographic_projection + -from spatialpandas.geometry import MultiPolygon current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -84,12 +83,12 @@ def test_pole_point_inside_polygon_from_vertice_north(self): [vertices[3], vertices[0]]]) # Check if the North pole is inside the polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'North', face_edge_cart) self.assertTrue(result, "North pole should be inside the polygon") # Check if the South pole is inside the polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'South', face_edge_cart) self.assertFalse(result, "South pole should not be inside the polygon") @@ -111,12 +110,12 @@ def test_pole_point_inside_polygon_from_vertice_south(self): [vertices[2], vertices[0]]]) # Check if the North pole is inside the polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'North', face_edge_cart) self.assertFalse(result, "North pole should not be inside the polygon") # Check if the South pole is inside the polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'South', face_edge_cart) self.assertTrue(result, "South pole should be inside the polygon") @@ -138,12 +137,12 @@ def test_pole_point_inside_polygon_from_vertice_pole(self): [vertices[3], vertices[0]]]) # Check if the North pole is inside the polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'North', face_edge_cart) self.assertTrue(result, "North pole should be inside the polygon") # Check if the South pole is inside the polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'South', face_edge_cart) self.assertFalse(result, "South pole should not be inside the polygon") @@ -164,7 +163,7 @@ def test_pole_point_inside_polygon_from_vertice_cross(self): [vertices[3], vertices[0]]]) # Check if the North pole is inside the polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'North', face_edge_cart) self.assertTrue(result, "North pole should be inside the polygon") @@ -354,7 +353,7 @@ def test_extreme_gca_latitude_max(self): ]) # Calculate the maximum latitude - max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max') + max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max') # Check if the maximum latitude is correct expected_max_latitude = self._max_latitude_rad_iterative(gca_cart) @@ -366,7 +365,7 @@ def test_extreme_gca_latitude_max(self): gca_cart = np.array([[0.0, 0.0, 1.0], [1.0, 0.0, 0.0]]) # Calculate the maximum latitude - max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max') + max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max') # Check if the maximum latitude is correct expected_max_latitude = np.pi / 2 # 90 degrees in radians @@ -379,7 +378,7 @@ def test_extreme_gca_latitude_max_short(self): gca_cart = np.array([[0.65465367, -0.37796447, -0.65465367], [0.6652466, -0.33896007, -0.6652466]]) # Calculate the maximum latitude - max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max') + max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max') # Check if the maximum latitude is correct expected_max_latitude = self._max_latitude_rad_iterative(gca_cart) @@ -395,7 +394,7 @@ def test_extreme_gca_latitude_min(self): ]) # Calculate the minimum latitude - min_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'min') + min_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'min') # Check if the minimum latitude is correct expected_min_latitude = self._min_latitude_rad_iterative(gca_cart) @@ -407,7 +406,7 @@ def test_extreme_gca_latitude_min(self): gca_cart = np.array([[0.0, 0.0, -1.0], [1.0, 0.0, 0.0]]) # Calculate the minimum latitude - min_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'min') + min_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'min') # Check if the minimum latitude is correct expected_min_latitude = -np.pi / 2 # 90 degrees in radians @@ -432,7 +431,7 @@ def test_insert_pt_in_latlonbox_non_periodic(self): old_box = np.array([[0.1, 0.2], [0.3, 0.4]]) # Radians new_pt = np.array([0.15, 0.35]) expected = np.array([[0.1, 0.2], [0.3, 0.4]]) - result = ux.grid.geometry._insert_pt_in_latlonbox( + result = ux.grid.geometry.insert_pt_in_latlonbox( old_box, new_pt, False) np.testing.assert_array_equal(result, expected) @@ -440,22 +439,22 @@ def test_insert_pt_in_latlonbox_periodic(self): old_box = np.array([[0.1, 0.2], [6.0, 0.1]]) # Radians, periodic new_pt = np.array([0.15, 6.2]) expected = np.array([[0.1, 0.2], [6.0, 0.1]]) - result = ux.grid.geometry._insert_pt_in_latlonbox(old_box, new_pt, True) + result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt, True) np.testing.assert_array_equal(result, expected) def test_insert_pt_in_latlonbox_pole(self): old_box = np.array([[0.1, 0.2], [0.3, 0.4]]) - new_pt = np.array([np.pi / 2, INT_FILL_VALUE]) # Pole point + new_pt = np.array([np.pi / 2, np.nan]) # Pole point expected = np.array([[0.1, np.pi / 2], [0.3, 0.4]]) - result = ux.grid.geometry._insert_pt_in_latlonbox(old_box, new_pt) + result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt) np.testing.assert_array_equal(result, expected) def test_insert_pt_in_empty_state(self): - old_box = np.array([[INT_FILL_VALUE, INT_FILL_VALUE], - [INT_FILL_VALUE, INT_FILL_VALUE]]) # Empty state + old_box = np.array([[np.nan, np.nan], + [np.nan, np.nan]]) # Empty state new_pt = np.array([0.15, 0.35]) expected = np.array([[0.15, 0.15], [0.35, 0.35]]) - result = ux.grid.geometry._insert_pt_in_latlonbox(old_box, new_pt) + result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt) np.testing.assert_array_equal(result, expected) @@ -593,9 +592,9 @@ def test_populate_bounds_normal(self): vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T # vertices_cart = [_lonlat_rad_to_xyz(vertices_rad[0], vertices_rad[1])] lat_max = max(np.deg2rad(60.0), - extreme_gca_latitude(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")) + _extreme_gca_latitude_cartesian(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")) lat_min = min(np.deg2rad(10.0), - extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) + _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) lon_min = np.deg2rad(10.0) lon_max = np.deg2rad(50.0) grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) @@ -622,9 +621,9 @@ def test_populate_bounds_antimeridian(self): vertices_rad = np.radians(vertices_lonlat) vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = max(np.deg2rad(60.0), - extreme_gca_latitude(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")) + _extreme_gca_latitude_cartesian(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")) lat_min = min(np.deg2rad(10.0), - extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) + _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) lon_min = np.deg2rad(350.0) lon_max = np.deg2rad(50.0) grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) @@ -745,7 +744,7 @@ def test_populate_bounds_node_on_pole(self): vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = min(np.deg2rad(10.0), - extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) + _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) lon_min = np.deg2rad(10.0) lon_max = np.deg2rad(50.0) grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) @@ -773,7 +772,7 @@ def test_populate_bounds_edge_over_pole(self): vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = min(np.deg2rad(60.0), - extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) + _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) lon_min = np.deg2rad(210.0) lon_max = np.deg2rad(30.0) grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) @@ -801,7 +800,7 @@ def test_populate_bounds_pole_inside(self): vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = min(np.deg2rad(60.0), - extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) + _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) lon_min = 0 lon_max = 2 * np.pi grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) @@ -1396,7 +1395,10 @@ def test_populate_bounds_GCAList_mix(self): [[np.deg2rad(60.0), np.pi / 2], [0., 2 * np.pi]]] grid = ux.Grid.from_face_vertices(faces, latlon=True) - bounds_xarray = _populate_bounds(grid, is_face_GCA_list=[[True, False, True, False]] * 4, return_array=True) + bounds_xarray = _populate_bounds(grid, is_face_GCA_list=np.array([[True, False, True, False], + [True, False, True, False], + [True, False, True, False], + [True, False, True, False]]) , return_array=True) face_bounds = bounds_xarray.values for i in range(len(faces)): nt.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE) diff --git a/test/test_helpers.py b/test/test_helpers.py index c5b923a26..da9b8c3eb 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -15,7 +15,7 @@ from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad from uxarray.grid.arcs import point_within_gca, _angle_of_2_vectors, in_between from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes -from uxarray.grid.geometry import _pole_point_inside_polygon +from uxarray.grid.geometry import pole_point_inside_polygon, _pole_point_inside_polygon_cartesian try: import constants @@ -267,7 +267,7 @@ def test_get_cartesian_face_edge_nodes_pipeline(self): ) # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'North', face_edges_connectivity_cartesian[0] ) @@ -300,7 +300,7 @@ def test_get_cartesian_face_edge_nodes_filled_value(self): ) # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'North', face_edges_connectivity_cartesian[0] ) @@ -422,7 +422,7 @@ def test_get_lonlat_face_edge_nodes_pipeline(self): face_edges_connectivity_cartesian.append(edge_cart) # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'North', np.array(face_edges_connectivity_cartesian) ) @@ -461,7 +461,7 @@ def test_get_lonlat_face_edge_nodes_filled_value(self): face_edges_connectivity_cartesian.append(edge_cart) # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon - result = ux.grid.geometry._pole_point_inside_polygon( + result = _pole_point_inside_polygon_cartesian( 'North', np.array(face_edges_connectivity_cartesian) ) diff --git a/test/test_intersections.py b/test/test_intersections.py index 0f3c62607..c28c5e0bc 100644 --- a/test/test_intersections.py +++ b/test/test_intersections.py @@ -5,8 +5,9 @@ # from uxarray.grid.coordinates import node_lonlat_rad_to_xyz, node_xyz_to_lonlat_rad +from uxarray.grid.arcs import _extreme_gca_latitude_cartesian from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad -from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection +from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection, _gca_gca_intersection_cartesian class TestGCAGCAIntersection(TestCase): @@ -24,10 +25,10 @@ def test_get_GCA_GCA_intersections_antimeridian(self): _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0), _lonlat_rad_to_xyz(np.deg2rad(179.0), 0.0) ]) - res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart) # res_cart should be empty since these two GCRs are not intersecting - self.assertTrue(np.array_equal(res_cart, np.array([]))) + self.assertTrue(len(res_cart) == 0) GCR1_cart = np.array([ _lonlat_rad_to_xyz(np.deg2rad(170.0), @@ -40,7 +41,7 @@ def test_get_GCA_GCA_intersections_antimeridian(self): _lonlat_rad_to_xyz(np.deg2rad(175.0), 0.0) ]) - res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart) res_cart = res_cart[0] # Test if the result is normalized @@ -66,7 +67,7 @@ def test_get_GCA_GCA_intersections_parallel(self): _lonlat_rad_to_xyz(0.5 * np.pi, 0.0), _lonlat_rad_to_xyz(-0.5 * np.pi - 0.01, 0.0) ]) - res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart) res_cart = res_cart[0] expected_res = np.array(_lonlat_rad_to_xyz(0.5 * np.pi, 0.0)) # Test if two results are equal within the error tolerance @@ -84,7 +85,7 @@ def test_get_GCA_GCA_intersections_perpendicular(self): _lonlat_rad_to_xyz(*[0.5 * np.pi, 0.0]), _lonlat_rad_to_xyz(*[-0.5 * np.pi - 0.01, 0.0]) ]) - res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart) res_cart = res_cart[0] # Test if the result is normalized self.assertTrue( @@ -133,7 +134,7 @@ def test_GCA_constLat_intersections_two_pts(self): _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(10.0)) ]) - max_lat = ux.grid.arcs.extreme_gca_latitude(GCR1_cart, 'max') + max_lat = _extreme_gca_latitude_cartesian(GCR1_cart, 'max') query_lat = (np.deg2rad(10.0) + max_lat) / 2.0 diff --git a/uxarray/__init__.py b/uxarray/__init__.py index 857c65f6c..c9cdff9e9 100644 --- a/uxarray/__init__.py +++ b/uxarray/__init__.py @@ -1,4 +1,9 @@ import uxarray.constants +import sys +# +# # TODO: numba recursion limit ? + + from .core.api import open_grid, open_dataset, open_mfdataset from .core.dataset import UxDataset @@ -34,6 +39,7 @@ def disable_fma(): disable_fma() +sys.setrecursionlimit(10000) __all__ = ( diff --git a/uxarray/cross_sections/dataarray_accessor.py b/uxarray/cross_sections/dataarray_accessor.py index de35739ab..d8c357e55 100644 --- a/uxarray/cross_sections/dataarray_accessor.py +++ b/uxarray/cross_sections/dataarray_accessor.py @@ -17,63 +17,89 @@ def __repr__(self): prefix = "\n" methods_heading = "Supported Methods:\n" - methods_heading += " * constant_latitude(center_coord, k, element, **kwargs)\n" + methods_heading += " * constant_latitude(lat)\n" + methods_heading += " * constant_longitude(lon)\n" return prefix + methods_heading - def constant_latitude(self, lat: float, use_spherical_bounding_box=False): - """Extracts a cross-section of the data array at a specified constant - latitude. + def constant_latitude(self, lat: float): + """Extracts a cross-section of the data array by selecting all faces that + intersect with a specified line of constant latitude. Parameters ---------- lat : float The latitude at which to extract the cross-section, in degrees. - use_spherical_bounding_box : bool, optional - If True, uses a spherical bounding box for intersection calculations. + Must be between -90.0 and 90.0 + + Returns + ------- + uxarray.UxDataArray + A subset of the original data array containing only the faces that intersect + with the specified latitude. Raises ------ ValueError - If no intersections are found at the specified latitude, a ValueError is raised. + If no intersections are found at the specified longitude or the data variable is not face-centered. Examples -------- - >>> uxda.constant_latitude_cross_section(lat=-15.5) + >>> # Extract data at 15.5°S latitude + >>> cross_section = uxda.constant_latitude(lat=-15.5) + + Notes + ----- + The initial execution time may be significantly longer than subsequent runs + due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching. """ - faces = self.uxda.uxgrid.get_faces_at_constant_latitude( - lat, use_spherical_bounding_box - ) + if not self.uxda._face_centered(): + raise ValueError( + "Cross sections are only supported for face-centered data variables." + ) + + faces = self.uxda.uxgrid.get_faces_at_constant_latitude(lat) return self.uxda.isel(n_face=faces) - def constant_longitude(self, lon: float, use_spherical_bounding_box=False): - """Extracts a cross-section of the data array at a specified constant - longitude. + def constant_longitude(self, lon: float): + """Extracts a cross-section of the data array by selecting all faces that + intersect with a specified line of constant longitude. Parameters ---------- lon : float - The longitude at which to extract the cross-section, in degrees. - use_spherical_bounding_box : bool, optional - If True, uses a spherical bounding box for intersection calculations. + The latitude at which to extract the cross-section, in degrees. + Must be between -180.0 and 180.0 + + Returns + ------- + uxarray.UxDataArray + A subset of the original data array containing only the faces that intersect + with the specified longitude. Raises ------ ValueError - If no intersections are found at the specified longitude, a ValueError is raised. + If no intersections are found at the specified longitude or the data variable is not face-centered. Examples -------- - >>> uxda.constant_longitude_cross_section(lon=-15.5) + >>> # Extract data at 0° longitude + >>> cross_section = uxda.constant_latitude(lon=0.0) Notes ----- - The accuracy and performance of the function can be controlled using the `method` parameter. - For higher precision requreiments, consider using method='acurate'. + The initial execution time may be significantly longer than subsequent runs + due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching. """ + if not self.uxda._face_centered(): + raise ValueError( + "Cross sections are only supported for face-centered data variables." + ) + faces = self.uxda.uxgrid.get_faces_at_constant_longitude( - lon, use_spherical_bounding_box + lon, ) return self.uxda.isel(n_face=faces) diff --git a/uxarray/cross_sections/grid_accessor.py b/uxarray/cross_sections/grid_accessor.py index 5b9055166..146da10af 100644 --- a/uxarray/cross_sections/grid_accessor.py +++ b/uxarray/cross_sections/grid_accessor.py @@ -17,54 +17,50 @@ def __repr__(self): prefix = "\n" methods_heading = "Supported Methods:\n" - methods_heading += " * constant_latitude(lat, )\n" + methods_heading += " * constant_latitude(lat, return_face_indices)\n" + methods_heading += " * constant_longitude(lon, return_face_indices)\n" return prefix + methods_heading def constant_latitude( - self, lat: float, return_face_indices=False, use_spherical_bounding_box=False + self, + lat: float, + return_face_indices: bool = False, ): - """Extracts a cross-section of the grid at a specified constant - latitude. + """Extracts a cross-section of the grid by selecting all faces that + intersect with a specified line of constant latitude. Parameters ---------- lat : float The latitude at which to extract the cross-section, in degrees. + Must be between -90.0 and 90.0 return_face_indices : bool, optional - If True, returns both the grid at the specified latitude and the indices - of the intersecting faces. If False, only the grid is returned. - Default is False. - use_spherical_bounding_box : bool, optional - If True, uses a spherical bounding box for intersection calculations. - + If True, also returns the indices of the faces that intersect with the line of constant latitude. Returns ------- - grid_at_constant_lat : Grid - The grid with faces that interesected at a given lattitude - faces : array, optional - The indices of the faces that intersect with the specified latitude. This is only - returned if `return_face_indices` is set to True. + uxarray.Grid + A subset of the original grid containing only the faces that intersect + with the specified latitude. Raises ------ ValueError - If no intersections are found at the specified latitude, a ValueError is raised. + If no intersections are found at the specified latitude. Examples -------- - >>> grid, indices = grid.cross_section.constant_latitude( - ... lat=30.0, return_face_indices=True - ... ) - >>> grid = grid.cross_section.constant_latitude(lat=-15.5) + >>> # Extract data at 15.5°S latitude + >>> cross_section = grid.constant_latitude(lat=-15.5) Notes ----- - The accuracy and performance of the function can be controlled using the `method` parameter. - For higher precision requreiments, consider using method='acurate'. + The initial execution time may be significantly longer than subsequent runs + due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching. """ + faces = self.uxgrid.get_faces_at_constant_latitude( - lat, use_spherical_bounding_box + lat, ) if len(faces) == 0: @@ -78,49 +74,44 @@ def constant_latitude( return grid_at_constant_lat def constant_longitude( - self, lon: float, use_spherical_bounding_box=False, return_face_indices=False + self, + lon: float, + return_face_indices: bool = False, ): - """Extracts a cross-section of the grid at a specified constant - longitude. + """Extracts a cross-section of the grid by selecting all faces that + intersect with a specified line of constant longitude. Parameters ---------- lon : float The longitude at which to extract the cross-section, in degrees. - use_spherical_bounding_box : bool, optional - If True, uses a spherical bounding box for intersection calculations. + Must be between -90.0 and 90.0 return_face_indices : bool, optional - If True, returns both the grid at the specified longitude and the indices - of the intersecting faces. If False, only the grid is returned. - Default is False. + If True, also returns the indices of the faces that intersect with the line of constant longitude. Returns ------- - grid_at_constant_lon : Grid - The grid with faces that interesected at a given longitudes - faces : array, optional - The indices of the faces that intersect with the specified longitude. This is only - returned if `return_face_indices` is set to True. + uxarray.Grid + A subset of the original grid containing only the faces that intersect + with the specified longitude. Raises ------ ValueError - If no intersections are found at the specified longitude, a ValueError is raised. + If no intersections are found at the specified longitude. Examples -------- - >>> grid, indices = grid.cross_section.constant_longitude( - ... lat=0.0, return_face_indices=True - ... ) - >>> grid = grid.cross_section.constant_longitude(lat=20.0) + >>> # Extract data at 0° longitude + >>> cross_section = grid.constant_latitude(lon=0.0) Notes ----- - The accuracy and performance of the function can be controlled using the `method` parameter. - For higher precision requreiments, consider using method='acurate'. + The initial execution time may be significantly longer than subsequent runs + due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching. """ faces = self.uxgrid.get_faces_at_constant_longitude( - lon, use_spherical_bounding_box + lon, ) if len(faces) == 0: diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 6be158005..3ccb4501b 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -1,9 +1,10 @@ import numpy as np +import math -# from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place + +from uxarray.grid.coordinates import _xyz_to_lonlat_rad_scalar from uxarray.grid.coordinates import ( - _xyz_to_lonlat_rad_scalar, _normalize_xyz_scalar, ) @@ -11,7 +12,7 @@ from uxarray.constants import ERROR_TOLERANCE, MACHINE_EPSILON -from uxarray.utils.computing import isclose, cross, dot, allclose +from uxarray.utils.computing import isclose, dot from numba import njit @@ -27,111 +28,206 @@ def _to_list(obj): return obj +def _point_within_gca_cartesian(pt_xyz, gca_xyz, is_directed=False): + pt_xyz = np.asarray(pt_xyz) + gca_xyz = np.asarray(gca_xyz) + + pt_lonlat = np.array( + _xyz_to_lonlat_rad_scalar(pt_xyz[0], pt_xyz[1], pt_xyz[2], normalize=False) + ) + gca_a_xyz = gca_xyz[0] + + gca_a_lonlat = np.array( + _xyz_to_lonlat_rad_scalar( + gca_xyz[0][0], gca_xyz[0][1], gca_xyz[0][2], normalize=False + ) + ) + gca_b_xyz = gca_xyz[1] + + gca_b_lonlat = np.array( + _xyz_to_lonlat_rad_scalar( + gca_xyz[1][0], gca_xyz[1][1], gca_xyz[1][2], normalize=False + ) + ) + + return point_within_gca( + pt_xyz, + pt_lonlat, + gca_a_xyz, + gca_a_lonlat, + gca_b_xyz, + gca_b_lonlat, + is_directed=is_directed, + ) + + @njit(cache=True) -def _point_within_gca_body( - angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed +def point_within_gca( + pt_xyz, + pt_lonlat, + gca_a_xyz, + gca_a_lonlat, + gca_b_xyz, + gca_b_lonlat, + is_directed=False, ): - angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1]) - if allclose(angle, np.pi, rtol=0.0, atol=MACHINE_EPSILON): + """Check if a point lies on a given Great Circle Arc (GCA). The anti- + meridian case is also considered. + + Parameters + ---------- + pt : numpy.ndarray (float) + Cartesian coordinates of the point. + gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr) + Cartesian coordinates of the Great Circle Arc (GCR). + 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 + ------- + bool + True if the point lies between the two endpoints of the GCR, False otherwise. + + Raises + ------ + ValueError + If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes. + In such cases, consider breaking the GCR into two separate GCRs. + + ValueError + If the input GCR spans more than 180 degrees (π radians). + In such cases, consider breaking the GCR into two separate GCRs. + + Notes + ----- + The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and + accounting for the anti-meridian case. + + The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude). + In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation. + + The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations. + + Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. + """ + + # ================================================================================================================== + # 1. If the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes + angle = _angle_of_2_vectors(gca_a_xyz, gca_b_xyz) + + if np.allclose(angle, np.pi, rtol=0.0, atol=MACHINE_EPSILON): raise ValueError( "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. " "Consider breaking the Great Circle Arc" "into two Great Circle Arcs" ) - + # ================================================================================================================== # See if the point is on the plane of the GCA, because we are dealing with floating point numbers with np.dot now # just using the rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON, but consider using the more proper error tolerance # in the future - cross_product = cross(np.asarray(gca_cart[0]), np.asarray(gca_cart[1])) - if not allclose( - dot(np.asarray(cross_product), np.asarray(pt)), # Custom dot function + cross_product = np.cross(gca_a_xyz, gca_b_xyz) + + if not np.allclose( + np.dot(cross_product, pt_xyz), 0, rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON, ): return False - - if isclose( - GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON + # ================================================================================================================== + if np.isclose( + gca_a_lonlat[0], gca_b_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON ): # If the pt and the GCA are on the same longitude (the y coordinates are the same) - if isclose( - GCRv0_lonlat[0], pt_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON + if np.isclose( + gca_a_lonlat[0], pt_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON ): # Now use the latitude to determine if the pt falls between the interval - return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1]) + return in_between(gca_a_lonlat[1], pt_lonlat[1], gca_b_lonlat[1]) else: - # 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 + # If the pt and the GCA are not on the same longitude when the GCA is a longitude 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 the longitude span is exactly 180 degree, then the GCA goes through the pole point # Or if one of the endpoints is on the pole point, then the GCA goes through the pole point if ( - isclose( - abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), + np.isclose( + np.abs(gca_b_lonlat[0] - gca_a_lonlat[0]), np.pi, rtol=0.0, atol=MACHINE_EPSILON, ) - or isclose( - abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + or np.isclose( + np.abs(gca_a_lonlat[1]), + np.pi / 2, + rtol=ERROR_TOLERANCE, + atol=ERROR_TOLERANCE, ) - or isclose( - abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + or np.isclose( + np.abs(gca_b_lonlat[1]), + np.pi / 2, + rtol=ERROR_TOLERANCE, + atol=ERROR_TOLERANCE, ) ): + # ============================================================================================================== # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0] # Since the point is our calculated properly, we use the atol=ERROR_TOLERANCE and rtol=ERROR_TOLERANCE - if isclose( - abs(pt_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + if np.isclose( + np.abs(pt_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE ): - pt_lonlat[0] = GCRv0_lonlat[0] + pt_lonlat[0] = gca_a_lonlat[0] + # ============================================================================================================== # Special case, if one of the GCA endpoints is on the pole point, and another endpoint is not # then we need to check if the pt is on the GCA - if isclose( - abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 - ) or isclose(abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0): + if np.isclose( + abs(gca_a_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 + ) or np.isclose( + abs(gca_b_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 + ): # Identify the non-pole endpoint non_pole_endpoint = None - if not isclose( - abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 + if not np.isclose( + abs(gca_a_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 ): - non_pole_endpoint = GCRv0_lonlat - elif not isclose( - abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 + non_pole_endpoint = gca_a_lonlat + elif not np.isclose( + abs(gca_b_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 ): - non_pole_endpoint = GCRv1_lonlat + non_pole_endpoint = gca_b_lonlat - if non_pole_endpoint is not None and not isclose( + if non_pole_endpoint is not None and not np.isclose( non_pole_endpoint[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 ): return False - - if not isclose( - GCRv0_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 - ) and not isclose( - GCRv1_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 + # ============================================================================================================== + if not np.isclose( + gca_a_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 + ) and not np.isclose( + gca_b_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 ): return False else: # Determine the pole latitude and latitude extension - if (GCRv0_lonlat[1] > 0.0 and GCRv1_lonlat[1] > 0.0) or ( - GCRv0_lonlat[1] < 0.0 and GCRv1_lonlat[1] < 0.0 + if (gca_a_lonlat[1] > 0.0 and gca_b_lonlat[1] > 0.0) or ( + gca_a_lonlat[1] < 0.0 and gca_b_lonlat[1] < 0.0 ): - pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0.0 else -np.pi / 2 + pole_lat = np.pi / 2 if gca_a_lonlat[1] > 0.0 else -np.pi / 2 lat_extend = ( - abs(np.pi / 2 - abs(GCRv0_lonlat[1])) + abs(np.pi / 2 - abs(gca_a_lonlat[1])) + np.pi / 2 - + abs(GCRv1_lonlat[1]) + + abs(gca_b_lonlat[1]) ) else: - pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1]) + pole_lat = _decide_pole_latitude(gca_a_lonlat[1], gca_b_lonlat[1]) lat_extend = ( - abs(np.pi / 2 - abs(GCRv0_lonlat[1])) + abs(np.pi / 2 - abs(gca_a_lonlat[1])) + np.pi / 2 - + abs(GCRv1_lonlat[1]) + + abs(gca_b_lonlat[1]) ) if is_directed and lat_extend >= np.pi: @@ -139,109 +235,44 @@ def _point_within_gca_body( "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] + return in_between(gca_a_lonlat[1], pt_lonlat[1], pole_lat) or in_between( + pole_lat, pt_lonlat[1], gca_b_lonlat[1] ) - + # endIf + # ============================================================================================================== 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) > π - if abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]) > np.pi: + if abs(gca_b_lonlat[0] - gca_a_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) - if GCRv0_lonlat[0] <= np.pi <= GCRv1_lonlat[0]: + if gca_a_lonlat[0] <= np.pi <= gca_b_lonlat[0]: raise ValueError( "The input Great Circle Arc span is larger than 180 degree, please break it into two." ) # 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 lon - elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0.0: + elif 2 * np.pi > gca_a_lonlat[0] > np.pi > gca_b_lonlat[0] > 0.0: return in_between( - GCRv0_lonlat[0], pt_lonlat[0], 2 * np.pi - ) or in_between(0, pt_lonlat[0], GCRv1_lonlat[0]) + gca_a_lonlat[0], pt_lonlat[0], 2 * np.pi + ) or in_between(0, pt_lonlat[0], gca_b_lonlat[0]) # The non-anti-meridian case. else: - return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) + return in_between(gca_a_lonlat[0], pt_lonlat[0], gca_b_lonlat[0]) else: # The undirected case # sort the longitude - GCRv0_lonlat_min, GCRv1_lonlat_max = sorted([GCRv0_lonlat[0], GCRv1_lonlat[0]]) + GCRv0_lonlat_min, GCRv1_lonlat_max = sorted([gca_a_lonlat[0], gca_b_lonlat[0]]) if np.pi > GCRv1_lonlat_max - GCRv0_lonlat_min >= 0.0: - return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) + return in_between(gca_a_lonlat[0], pt_lonlat[0], gca_b_lonlat[0]) else: return in_between(GCRv1_lonlat_max, pt_lonlat[0], 2 * np.pi) or in_between( 0.0, pt_lonlat[0], GCRv0_lonlat_min ) - -def point_within_gca(pt, gca_cart, is_directed=False): - """Check if a point lies on a given Great Circle Arc (GCA). The anti- - meridian case is also considered. - - Parameters - ---------- - pt : numpy.ndarray (float) - Cartesian coordinates of the point. - gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr) - Cartesian coordinates of the Great Circle Arc (GCR). - 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 - ------- - bool - True if the point lies between the two endpoints of the GCR, False otherwise. - - Raises - ------ - ValueError - If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes. - In such cases, consider breaking the GCR into two separate GCRs. - - ValueError - If the input GCR spans more than 180 degrees (π radians). - In such cases, consider breaking the GCR into two separate GCRs. - - Notes - ----- - The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and - accounting for the anti-meridian case. - - The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude). - In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation. - - The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations. - - Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. - """ - # Convert the cartesian coordinates to lonlat coordinates - pt_lonlat = np.array( - _xyz_to_lonlat_rad_scalar(pt[0], pt[1], pt[2], normalize=False) - ) - GCRv0_lonlat = np.array( - _xyz_to_lonlat_rad_scalar( - gca_cart[0][0], gca_cart[0][1], gca_cart[0][2], normalize=False - ) - ) - GCRv1_lonlat = np.array( - _xyz_to_lonlat_rad_scalar( - gca_cart[1][0], gca_cart[1][1], gca_cart[1][2], normalize=False - ) - ) - gca_cart = np.asarray(gca_cart) - - # 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(gca_cart[0], gca_cart[1]) - - out = _point_within_gca_body( - angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed - ) - - return out + return None @njit(cache=True) @@ -306,8 +337,52 @@ def _decide_pole_latitude(lat1, lat2): return closest_pole -def extreme_gca_latitude(gca_cart, extreme_type): - """Calculate the maximum or minimum latitude of a great circle arc defined +@njit(cache=True) +def max3(a, b, c): + if a >= b and a >= c: + return a + elif b >= c: + return b + else: + return c + + +@njit(cache=True) +def min3(a, b, c): + if a <= b and a <= c: + return a + elif b <= c: + return b + else: + return c + + +@njit(cache=True) +def clip_scalar(a, a_min, a_max): + if a < a_min: + return a_min + elif a > a_max: + return a_max + else: + return a + + +def _extreme_gca_latitude_cartesian(gca_cart, extreme_type): + # should be shape [2, 2] + gca_lonlat = np.array( + [ + _xyz_to_lonlat_rad_scalar(gca_cart[0, 0], gca_cart[0, 1], gca_cart[0, 2]), + _xyz_to_lonlat_rad_scalar(gca_cart[1, 0], gca_cart[1, 1], gca_cart[1, 2]), + ] + ) + + return extreme_gca_latitude(gca_cart, gca_lonlat, extreme_type) + + +@njit(cache=True) +def extreme_gca_latitude(gca_cart, gca_lonlat, extreme_type): + """ + Calculate the maximum or minimum latitude of a great circle arc defined by two 3D points. Parameters @@ -315,6 +390,9 @@ def extreme_gca_latitude(gca_cart, extreme_type): gca_cart : numpy.ndarray An array containing two 3D vectors that define a great circle arc. + gca_lonlat : numpy.ndarray + An array containing the longitude and latitude of the two points. + extreme_type : str The type of extreme latitude to calculate. Must be either 'max' or 'min'. @@ -328,36 +406,57 @@ def extreme_gca_latitude(gca_cart, extreme_type): ValueError If `extreme_type` is not 'max' or 'min'. """ - extreme_type = extreme_type.lower() - - if extreme_type not in ("max", "min"): + # Validate extreme_type + if (extreme_type != "max") and (extreme_type != "min"): raise ValueError("extreme_type must be either 'max' or 'min'") - n1, n2 = gca_cart + # Extract the two points + n1 = gca_cart[0] + n2 = gca_cart[1] + + # Compute dot product + dot_n1_n2 = dot(n1, n2) - dot_n1_n2 = dot(np.asarray(n1), np.asarray(n2)) + # Compute denominator denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0) - d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom - d_a_max = ( - np.clip(d_a_max, 0, 1) - if isclose(d_a_max, 0, atol=ERROR_TOLERANCE) - or isclose(d_a_max, 1, atol=ERROR_TOLERANCE) - else d_a_max - ) - # Before we make sure the grid coordinates are normalized, do not try to skip the normalization steps! - _, lat_n1 = _xyz_to_lonlat_rad_scalar(n1[0], n1[1], n1[2], normalize=True) - _, lat_n2 = _xyz_to_lonlat_rad_scalar(n2[0], n2[1], n2[2], normalize=True) - - if 0 < d_a_max < 1: - node3 = (1 - d_a_max) * n1 + d_a_max * n2 - node3 = np.array(_normalize_xyz_scalar(node3[0], node3[1], node3[2])) - d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1)) - - return ( - max(d_lat_rad, lat_n1, lat_n2) - if extreme_type == "max" - else min(d_lat_rad, lat_n1, lat_n2) - ) + # Initialize latitudes + lon_n1, lat_n1 = gca_lonlat[0] + lon_n2, lat_n2 = gca_lonlat[1] + + # Check if denominator is zero to avoid division by zero + if denom != 0.0: + d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom + + # Handle cases where d_a_max is very close to 0 or 1 + if isclose(d_a_max, 0.0, atol=ERROR_TOLERANCE) or isclose( + d_a_max, 1.0, atol=ERROR_TOLERANCE + ): + d_a_max = clip_scalar(d_a_max, 0.0, 1.0) + + # Check if d_a_max is within the valid range + if (d_a_max > 0.0) and (d_a_max < 1.0): + # Compute the intermediate point on the GCA + node3 = (1.0 - d_a_max) * n1 + d_a_max * n2 + + # Normalize the intermediate point + x, y, z = _normalize_xyz_scalar(node3[0], node3[1], node3[2]) + node3_normalized = np.empty(3) + node3_normalized[0] = x + node3_normalized[1] = y + node3_normalized[2] = z + + # Compute latitude of the intermediate point + d_lat_rad = math.asin(clip_scalar(node3_normalized[2], -1.0, 1.0)) + + # Return the extreme latitude + if extreme_type == "max": + return max3(d_lat_rad, lat_n1, lat_n2) + else: + return min3(d_lat_rad, lat_n1, lat_n2) + + # If denom is zero or d_a_max is not in (0,1), return max or min of lat_n1 and lat_n2 + if extreme_type == "max": + return max(lat_n1, lat_n2) else: - return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2) + return min(lat_n1, lat_n2) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index cb3b2cadf..25e1d0101 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -1,6 +1,7 @@ import xarray as xr import numpy as np import warnings +import math from uxarray.conventions import ugrid @@ -67,36 +68,10 @@ def _xyz_to_lonlat_rad_no_norm( @njit(cache=True) -def _xyz_to_lonlat_rad_scalar( - x: Union[np.ndarray, float], - y: Union[np.ndarray, float], - z: Union[np.ndarray, float], - normalize: bool = True, -): - """Converts a Cartesian x,y,z coordinates into Spherical latitude and - longitude without normalization, decorated with Numba. - - Parameters - ---------- - x : float - Cartesian x coordinate - y: float - Cartesiain y coordinate - z: float - Cartesian z coordinate - - - Returns - ------- - lon : float - Longitude in radians - lat: float - Latitude in radians - """ - +def _xyz_to_lonlat_rad_scalar(x, y, z, normalize=True): if normalize: x, y, z = _normalize_xyz_scalar(x, y, z) - denom = np.abs(x * x + y * y + z * z) + denom = abs(x * x + y * y + z * z) x /= denom y /= denom z /= denom @@ -104,15 +79,19 @@ def _xyz_to_lonlat_rad_scalar( lon = np.atan2(y, x) lat = np.asin(z) - # set longitude range to [0, pi] - lon = np.mod(lon, 2 * np.pi) + # Set longitude range to [0, 2*pi] + lon = lon % (2 * math.pi) - z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE + z_abs = abs(z) + if z_abs > (1.0 - ERROR_TOLERANCE): + lat = math.copysign(math.pi / 2, z) + lon = 0.0 - lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) - lon = np.where(z_mask, 0.0, lon) + lonlat = np.empty(2) + lonlat[0] = lon + lonlat[1] = lat - return lon, lat + return lonlat def _xyz_to_lonlat_rad( diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index ed564aafb..8f05ee666 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -2,7 +2,7 @@ import cartopy.crs as ccrs import geopandas from matplotlib.collections import LineCollection, PolyCollection -from numba import njit +from numba import njit, prange import numpy as np import pandas as pd import shapely @@ -11,6 +11,7 @@ import spatialpandas from spatialpandas.geometry import MultiPolygonArray, PolygonArray import xarray as xr +import math from uxarray.constants import ( ERROR_TOLERANCE, @@ -18,6 +19,7 @@ INT_FILL_VALUE, ) from uxarray.grid.arcs import extreme_gca_latitude, point_within_gca +from uxarray.grid.coordinates import _xyz_to_lonlat_rad from uxarray.grid.intersections import gca_gca_intersection from uxarray.grid.utils import ( _get_cartesian_face_edge_nodes, @@ -25,62 +27,67 @@ ) from uxarray.utils.computing import allclose, isclose -POLE_POINTS = {"North": np.array([0.0, 0.0, 1.0]), "South": np.array([0.0, 0.0, -1.0])} +POLE_POINTS_XYZ = { + "North": np.array([0.0, 0.0, 1.0]), + "South": np.array([0.0, 0.0, -1.0]), +} +POLE_POINTS_LONLAT = { + "North": np.array([0.0, np.pi / 2]), + "South": np.array([0.0, -np.pi / 2]), +} # number of faces/polygons before raising a warning for performance GDF_POLYGON_THRESHOLD = 100000 -REFERENCE_POINT_EQUATOR = np.array([1.0, 0.0, 0.0]) +REFERENCE_POINT_EQUATOR_XYZ = np.array([1.0, 0.0, 0.0]) +REFERENCE_POINT_EQUATOR_LONLAT = np.array([0.0, 0.0]) # TODO ? +POLE_NAME_TO_INT = {"North": 1, "Equator": 0, "South": -1} + +@njit(cache=True) +def error_radius(p1, p2): + """Calculate the error radius between two points in 3D space.""" + numerator = math.sqrt( + (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2 + (p1[2] - p2[2]) ** 2 + ) + denominator = math.sqrt(p2[0] ** 2 + p2[1] ** 2 + p2[2] ** 2) + return numerator / denominator + + +@njit(cache=True) def _unique_points(points, tolerance=ERROR_TOLERANCE): """Identify unique intersection points from a list of points, considering floating point precision errors. Parameters ---------- - points : list of array-like - A list containing the intersection points, where each point is an array-like structure (e.g., list, tuple, or numpy array) containing x, y, and z coordinates. - tolerance : float, optional - The distance threshold within which two points are considered identical. Default is 1e-6. + points : np.ndarray + An array of shape (n_points, 3) containing the intersection points. + tolerance : float + The distance threshold within which two points are considered identical. Returns ------- - list of np.ndarray - A list of unique snapped points in Cartesian coordinates. - - Notes - ----- - Given the nature of the mathematical equations and the spherical surface, it is more reasonable to calculate the "error radius" of the results using the following formula. - In the equation below, \(\tilde{P}_x\) and \(\tilde{P}_y\) are the calculated results, and \(P_x\) and \(P_y\) are the actual intersection points for the \(x\) and \(y\) coordinates, respectively. - The \(z\) coordinate is always the input \(z_0\), the constant latitude, so it is not included in this error calculation. - - .. math:: - \begin{aligned} - &\frac{\sqrt{(\tilde{v}_x - v_x)^2 + (\tilde{v}_y - v_y)^2 + (\tilde{v}_z - v_z)^2}}{\sqrt{v_x^2 + v_y^2 + v_z^2}}\\ - &= \sqrt{(\tilde{v}_x - v_x)^2 + (\tilde{v}_y - v_y)^2 + (\tilde{v}_z - v_z)^2} - \end{aligned} - - This method ensures that small numerical inaccuracies do not lead to multiple close points being considered different. + np.ndarray + An array of unique points in Cartesian coordinates. """ - unique_points = [] - points = [np.array(point) for point in points] # Ensure all points are numpy arrays - - def error_radius(p1, p2): - """Calculate the error radius between two points in 3D space.""" - numerator = np.sqrt( - (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2 + (p1[2] - p2[2]) ** 2 - ) - denominator = np.sqrt(p2[0] ** 2 + p2[1] ** 2 + p2[2] ** 2) - return numerator / denominator + n_points = points.shape[0] + unique_points = np.empty((n_points, 3), dtype=points.dtype) + unique_count = 0 - for point in points: - if not any( - error_radius(point, unique_point) < tolerance - for unique_point in unique_points - ): - unique_points.append(point) + for i in range(n_points): + point = points[i] + is_unique = True + for j in range(unique_count): + unique_point = unique_points[j] + if error_radius(point, unique_point) < tolerance: + is_unique = False + break + if is_unique: + unique_points[unique_count] = point + unique_count += 1 - return unique_points + return unique_points[:unique_count] @njit(cache=True) @@ -670,145 +677,250 @@ def _convert_shells_to_polygons(shells): return polygons -def _pole_point_inside_polygon(pole, face_edge_cart): - """Determines if a pole point is inside a polygon. +def _pole_point_inside_polygon_cartesian(pole, face_edges_xyz): + if isinstance(pole, str): + pole = POLE_NAME_TO_INT[pole] - .. note:: - - If the pole point is on the edge of the polygon, it will be considered "inside the polygon". + x = face_edges_xyz[:, :, 0] + y = face_edges_xyz[:, :, 1] + z = face_edges_xyz[:, :, 2] - Parameters - ---------- - pole : str - Either 'North' or 'South'. - face_edge_cart : np.ndarray - A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3) + lon, lat = _xyz_to_lonlat_rad(x, y, z) - Returns - ------- - bool - True if pole point is inside polygon, False otherwise. + face_edges_lonlat = np.stack((lon, lat), axis=2) - Raises - ------ - ValueError - If the provided pole is neither 'North' nor 'South'. + return pole_point_inside_polygon(pole, face_edges_xyz, face_edges_lonlat) - Warning - ------- - UserWarning - Raised if the face contains both pole points. - """ - if pole not in POLE_POINTS: - raise ValueError('Pole point must be either "North" or "South"') + pass - # Classify the polygon's location - location = _classify_polygon_location(face_edge_cart) - pole_point = POLE_POINTS[pole] - if location == pole: - ref_edge = np.array([pole_point, REFERENCE_POINT_EQUATOR]) - return _check_intersection(ref_edge, face_edge_cart) % 2 != 0 - elif location == "Equator": - # smallest offset I can obtain when using the float64 type +@njit(cache=True) +def pole_point_inside_polygon(pole, face_edges_xyz, face_edges_lonlat): + """Determines if a pole point is inside a polygon.""" - ref_edge_north = np.array([pole_point, REFERENCE_POINT_EQUATOR]) - ref_edge_south = np.array([-pole_point, REFERENCE_POINT_EQUATOR]) + if pole != 1 and pole != -1: + raise ValueError("Pole must be 1 (North) or -1 (South)") - north_edges = face_edge_cart[np.any(face_edge_cart[:, :, 2] > 0, axis=1)] - south_edges = face_edge_cart[ - ~np.isin(face_edge_cart, north_edges).all(axis=(1, 2)) - ] - # The equator one is assigned to the south edges - return ( - _check_intersection(ref_edge_north, north_edges) - + _check_intersection(ref_edge_south, south_edges) - ) % 2 != 0 - elif ( - location == "North" - and pole == "South" - or location == "South" - and pole == "North" - ): + # Define constants within the function + pole_point_xyz = np.empty(3, dtype=np.float64) + pole_point_xyz[0] = 0.0 + pole_point_xyz[1] = 0.0 + pole_point_xyz[2] = 1.0 * pole + + pole_point_lonlat = np.empty(2, dtype=np.float64) + pole_point_lonlat[0] = 0.0 + pole_point_lonlat[1] = (math.pi / 2) * pole + + REFERENCE_POINT_EQUATOR_XYZ = np.empty(3, dtype=np.float64) + REFERENCE_POINT_EQUATOR_XYZ[0] = 1.0 + REFERENCE_POINT_EQUATOR_XYZ[1] = 0.0 + REFERENCE_POINT_EQUATOR_XYZ[2] = 0.0 + + REFERENCE_POINT_EQUATOR_LONLAT = np.empty(2, dtype=np.float64) + REFERENCE_POINT_EQUATOR_LONLAT[0] = 0.0 + REFERENCE_POINT_EQUATOR_LONLAT[1] = 0.0 + + # Classify the polygon's location + location = _classify_polygon_location(face_edges_xyz) + + if (location == 1 and pole == -1) or (location == -1 and pole == 1): return False + + elif location == -1 or location == 1: + # Initialize ref_edge_xyz + ref_edge_xyz = np.empty((2, 3), dtype=np.float64) + ref_edge_xyz[0, 0] = pole_point_xyz[0] + ref_edge_xyz[0, 1] = pole_point_xyz[1] + ref_edge_xyz[0, 2] = pole_point_xyz[2] + ref_edge_xyz[1, :] = REFERENCE_POINT_EQUATOR_XYZ + + # Initialize ref_edge_lonlat + ref_edge_lonlat = np.empty((2, 2), dtype=np.float64) + ref_edge_lonlat[0, 0] = pole_point_lonlat[0] + ref_edge_lonlat[0, 1] = pole_point_lonlat[1] + ref_edge_lonlat[1, :] = REFERENCE_POINT_EQUATOR_LONLAT + + intersection_count = _check_intersection( + ref_edge_xyz, ref_edge_lonlat, face_edges_xyz, face_edges_lonlat + ) + return (intersection_count % 2) != 0 + + elif location == 0: # Equator + # Initialize ref_edge_north_xyz and ref_edge_north_lonlat + ref_edge_north_xyz = np.empty((2, 3), dtype=np.float64) + ref_edge_north_xyz[0, 0] = 0.0 + ref_edge_north_xyz[0, 1] = 0.0 + ref_edge_north_xyz[0, 2] = 1.0 + ref_edge_north_xyz[1, :] = REFERENCE_POINT_EQUATOR_XYZ + + ref_edge_north_lonlat = np.empty((2, 2), dtype=np.float64) + ref_edge_north_lonlat[0, 0] = 0.0 + ref_edge_north_lonlat[0, 1] = math.pi / 2 + ref_edge_north_lonlat[1, :] = REFERENCE_POINT_EQUATOR_LONLAT + + # Initialize ref_edge_south_xyz and ref_edge_south_lonlat + ref_edge_south_xyz = np.empty((2, 3), dtype=np.float64) + ref_edge_south_xyz[0, 0] = 0.0 + ref_edge_south_xyz[0, 1] = 0.0 + ref_edge_south_xyz[0, 2] = -1.0 + ref_edge_south_xyz[1, :] = REFERENCE_POINT_EQUATOR_XYZ + + ref_edge_south_lonlat = np.empty((2, 2), dtype=np.float64) + ref_edge_south_lonlat[0, 0] = 0.0 + ref_edge_south_lonlat[0, 1] = -math.pi / 2 + ref_edge_south_lonlat[1, :] = REFERENCE_POINT_EQUATOR_LONLAT + + # Classify edges based on z-coordinate + n_edges = face_edges_xyz.shape[0] + north_edges_xyz = np.empty((n_edges, 2, 3), dtype=np.float64) + north_edges_lonlat = np.empty((n_edges, 2, 2), dtype=np.float64) + south_edges_xyz = np.empty((n_edges, 2, 3), dtype=np.float64) + south_edges_lonlat = np.empty((n_edges, 2, 2), dtype=np.float64) + north_count = 0 + south_count = 0 + + for i in range(n_edges): + edge_xyz = face_edges_xyz[i] + edge_lonlat = face_edges_lonlat[i] + if edge_xyz[0, 2] > 0 or edge_xyz[1, 2] > 0: + north_edges_xyz[north_count] = edge_xyz + north_edges_lonlat[north_count] = edge_lonlat + north_count += 1 + elif edge_xyz[0, 2] < 0 or edge_xyz[1, 2] < 0: + south_edges_xyz[south_count] = edge_xyz + south_edges_lonlat[south_count] = edge_lonlat + south_count += 1 + else: + # skip edges exactly on the equator + continue + + if north_count > 0: + north_edges_xyz = north_edges_xyz[:north_count] + north_edges_lonlat = north_edges_lonlat[:north_count] + else: + north_edges_xyz = np.empty((0, 2, 3), dtype=np.float64) + north_edges_lonlat = np.empty((0, 2, 2), dtype=np.float64) + + if south_count > 0: + south_edges_xyz = south_edges_xyz[:south_count] + south_edges_lonlat = south_edges_lonlat[:south_count] + else: + south_edges_xyz = np.empty((0, 2, 3), dtype=np.float64) + south_edges_lonlat = np.empty((0, 2, 2), dtype=np.float64) + + # Count south intersections + north_intersections = _check_intersection( + ref_edge_north_xyz, + ref_edge_north_lonlat, + north_edges_xyz, + north_edges_lonlat, + ) + + # Count south intersections + south_intersections = _check_intersection( + ref_edge_south_xyz, + ref_edge_south_lonlat, + south_edges_xyz, + south_edges_lonlat, + ) + + return ((north_intersections + south_intersections) % 2) != 0 + else: raise ValueError( - "Invalid pole point query. Current location: {}, query pole point: {}".format( - location, pole - ) + f"Invalid pole point query. Current location: {location}, query pole point: {pole}" ) -def _check_intersection(ref_edge, edges): - """Check the number of intersections of the reference edge with the given - edges. +@njit(cache=True) +def _classify_polygon_location(face_edge_cart): + """Classify the location of the polygon relative to the hemisphere.""" + z_coords = face_edge_cart[:, :, 2] + if np.all(z_coords > 0): # Use strict inequality + return 1 # North + elif np.all(z_coords < 0): # Use strict inequality + return -1 # South + else: + return 0 # Equator + + +@njit(cache=True) +def _check_intersection(ref_edge_xyz, ref_edge_lonlat, edges_xyz, edges_lonlat): + """Check the number of intersections of the reference edge with the given edges. Parameters ---------- - ref_edge : np.ndarray - Reference edge to check intersections against. - edges : np.ndarray + ref_edge_xyz : np.ndarray + Reference edge to check intersections against. Shape: (2, 3) + ref_edge_lonlat : np.ndarray + Reference edge longitude and latitude. Shape: (2, 2) + edges_xyz : np.ndarray Edges to check for intersections. Shape: (n_edges, 2, 3) + edges_lonlat : np.ndarray + Longitude and latitude of the edges. Shape: (n_edges, 2, 2) Returns ------- int Count of intersections. """ - pole_point, ref_point = ref_edge - intersection_points = [] + pole_point_xyz = ref_edge_xyz[0] + n_edges = edges_xyz.shape[0] - for edge in edges: - intersection_point = gca_gca_intersection(ref_edge, edge) + # Assuming at most 2 intersections per edge + max_intersections = n_edges * 2 + intersection_points = np.empty((max_intersections, 3), dtype=np.float64) + intersection_count = 0 - if intersection_point.size != 0: - # Handle both single point and multiple points case - if ( - intersection_point.ndim == 1 - ): # If there's only one point, make it a 2D array - intersection_point = [intersection_point] # Convert to list of points + for i in range(n_edges): + edge_xyz = edges_xyz[i] + edge_lonlat = edges_lonlat[i] - # for each intersection point, check if it is a pole point - for point in intersection_point: - if allclose(point, pole_point, atol=ERROR_TOLERANCE): - return True - intersection_points.append(point) + # compute intersection + intersection_point = gca_gca_intersection( + ref_edge_xyz, ref_edge_lonlat, edge_xyz, edge_lonlat + ) - # Only return the unique intersection points, the unique tolerance is set to ERROR_TOLERANCE - intersection_points = _unique_points(intersection_points, tolerance=ERROR_TOLERANCE) + if intersection_point.size != 0: + if intersection_point.ndim == 1: + # Only one point + point = intersection_point + if allclose(point, pole_point_xyz, atol=ERROR_TOLERANCE): + return True + intersection_points[intersection_count] = point + intersection_count += 1 + else: + # Multiple points + num_points = intersection_point.shape[0] + for j in range(num_points): + point = intersection_point[j] + if allclose(point, pole_point_xyz, atol=ERROR_TOLERANCE): + return True + intersection_points[intersection_count] = point + intersection_count += 1 + + if intersection_count == 0: + return 0 + + unique_intersection_points = _unique_points( + intersection_points[:intersection_count], tolerance=ERROR_TOLERANCE + ) + unique_count = unique_intersection_points.shape[0] - # If the unique intersection point is one and it is exactly one of the nodes of the face, return 0 - if len(intersection_points) == 1: - for edge in edges: + # If there's only one unique intersection point, check if it matches any edge nodes + if unique_count == 1: + intersection_point = unique_intersection_points[0] + for i in range(n_edges): + edge_xyz = edges_xyz[i] if allclose( - intersection_points[0], edge[0], atol=ERROR_TOLERANCE - ) or allclose(intersection_points[0], edge[1], atol=ERROR_TOLERANCE): + intersection_point, edge_xyz[0], atol=ERROR_TOLERANCE + ) or allclose(intersection_point, edge_xyz[1], atol=ERROR_TOLERANCE): return 0 - return len(intersection_points) - - -def _classify_polygon_location(face_edge_cart): - """Classify the location of the polygon relative to the hemisphere. - - Parameters - ---------- - face_edge_cart : np.ndarray - A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3) - - Returns - ------- - str - Returns either 'North', 'South' or 'Equator' based on the polygon's location. - """ - z_coords = face_edge_cart[:, :, 2] - if np.all(z_coords > 0): - return "North" - elif np.all(z_coords < 0): - return "South" - else: - return "Equator" + return unique_count +@njit(cache=True) def _get_latlonbox_width(latlonbox_rad): """Calculate the width of a latitude-longitude box in radians. The box should be represented by a 2x2 array in radians and lon0 represent the @@ -847,18 +959,16 @@ def _get_latlonbox_width(latlonbox_rad): If the input longitude range is flagged as periodic but in the form [lon0, lon1] where lon0 < lon1. The function will automatically use the is_lon_periodic=False instead. """ - + # TODO: ensure this function works properly lon0, lon1 = latlonbox_rad[1] - # Check longitude range validity - # Normalize the longitude so that it is within the range [0, 2π] - # Normalize the longitude if lon0 != INT_FILL_VALUE: lon0 = np.mod(lon0, 2 * np.pi) if lon1 != INT_FILL_VALUE: lon1 = np.mod(lon1, 2 * np.pi) if (lon0 < 0.0 or lon0 > 2.0 * np.pi) and lon0 != INT_FILL_VALUE: - raise Exception("lon0 out of range ({} not in [0, 2π])".format(lon0)) + # -1 used for exception + return -1 if lon0 <= lon1: return lon1 - lon0 @@ -867,18 +977,33 @@ def _get_latlonbox_width(latlonbox_rad): return 2 * np.pi - lon0 + lon1 -def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True): - """Update the latitude-longitude box to include a new point in radians. +@njit(cache=True) +def all_elements_nan(arr): + """Check if all elements in an array are np.nan.""" + for i in range(arr.shape[0]): + if not np.isnan(arr[i]): + return False + return True - This function compares the new point's latitude and longitude with the - existing latitude-longitude box and updates the box if necessary to include the new point. + +@njit(cache=True) +def any_close_lat(lat_pt, atol): + """Check if the latitude point is close to either the North or South Pole.""" + return isclose(lat_pt, 0.5 * math.pi, atol) or isclose(lat_pt, -0.5 * math.pi, atol) + + +@njit(cache=True) +def insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True): + """Update the latitude-longitude box to include a new point in radians. Parameters ---------- old_box : np.ndarray - The original latitude-longitude box in radian, a 2x2 array: [min_lat, max_lat],[left_lon, right_lon]]. + The original latitude-longitude box in radians, a 2x2 array: + [[min_lat, max_lat], + [left_lon, right_lon]]. new_pt : np.ndarray - The new latitude-longitude point in radian, an array: [lat, lon]. + The new latitude-longitude point in radians, an array: [lat, lon]. is_lon_periodic : bool, optional Flag indicating if the latitude-longitude box is periodic in longitude (default is True). @@ -891,301 +1016,340 @@ def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True): ------ Exception If logic errors occur in the calculation process. - - Examples - -------- - >>> _insert_pt_in_latlonbox( - ... np.array([[1.0, 2.0], [3.0, 4.0]]), np.array([1.5, 3.5]) - ... ) - array([[1.0, 2.0], [3.0, 4.0]]) """ - if np.all(new_pt == INT_FILL_VALUE): + # Check if the new point is a fill value + all_fill = all_elements_nan(new_pt) + if all_fill: return old_box - latlon_box = np.copy(old_box) # Create a copy of the old box - latlon_box = np.array( - latlon_box, dtype=np.float64 - ) # Cast to float64, otherwise the following update might fail + # Create a copy of the old box + latlon_box = np.copy(old_box) - lat_pt, lon_pt = new_pt + # Extract latitude and longitude from the new point + lat_pt = new_pt[0] + lon_pt = new_pt[1] - # Normalize the longitude - if lon_pt != INT_FILL_VALUE: - lon_pt = np.mod(lon_pt, 2 * np.pi) + # Normalize the longitude if it's not a fill value + if not np.isnan(lon_pt): + lon_pt = lon_pt % (2.0 * math.pi) # Check if the latitude range is uninitialized and update it - if old_box[0][0] == old_box[0][1] == INT_FILL_VALUE: - latlon_box[0] = np.array([lat_pt, lat_pt]) - - # Check if the longitude range is uninitialized and update it - if old_box[1][0] == old_box[1][1] == INT_FILL_VALUE: - latlon_box[1] = np.array([lon_pt, lon_pt]) - - if lon_pt != INT_FILL_VALUE and (lon_pt < 0.0 or lon_pt > 2.0 * np.pi): - raise Exception(f"lon_pt out of range ({lon_pt} not in [0, 2π])") - - # Check for pole points and update latitudes - is_pole_point = ( - lon_pt == INT_FILL_VALUE - and isclose( - new_pt[0], np.asarray([0.5 * np.pi, -0.5 * np.pi]), atol=ERROR_TOLERANCE - ).any() - ) - - if is_pole_point: - # Check if the new point is close to the North Pole - if isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE): - latlon_box[0][1] = 0.5 * np.pi - - # Check if the new point is close to the South Pole - elif isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE): - latlon_box[0][0] = -0.5 * np.pi - - return latlon_box + if np.isnan(old_box[0, 0]) and np.isnan(old_box[0, 1]): + latlon_box[0, 0] = lat_pt + latlon_box[0, 1] = lat_pt else: - latlon_box[0] = [min(latlon_box[0][0], lat_pt), max(latlon_box[0][1], lat_pt)] + # Update latitude range + if lat_pt < latlon_box[0, 0]: + latlon_box[0, 0] = lat_pt + if lat_pt > latlon_box[0, 1]: + latlon_box[0, 1] = lat_pt - # Update longitude range for non-periodic or periodic cases - if not is_lon_periodic: - latlon_box[1] = [min(latlon_box[1][0], lon_pt), max(latlon_box[1][1], lon_pt)] + # Check if the longitude range is uninitialized and update it + if np.isnan(old_box[1, 0]) and np.isnan(old_box[1, 1]): + latlon_box[1, 0] = lon_pt + latlon_box[1, 1] = lon_pt else: - if ( - latlon_box[1][0] > latlon_box[1][1] - and (lon_pt < latlon_box[1][0] and lon_pt > latlon_box[1][1]) - ) or ( - latlon_box[1][0] <= latlon_box[1][1] - and not (latlon_box[1][0] <= lon_pt <= latlon_box[1][1]) - ): - # Calculate and compare new box widths - box_a, box_b = np.copy(latlon_box), np.copy(latlon_box) - box_a[1][0], box_b[1][1] = lon_pt, lon_pt - d_width_a, d_width_b = ( - _get_latlonbox_width(box_a), - _get_latlonbox_width(box_b), - ) - - # The width should not be negative, if so, raise an exception - if d_width_a < 0 or d_width_b < 0: - raise Exception("logic error") - - # Return the arc with the smaller width - latlon_box = box_a if d_width_a < d_width_b else box_b + # Validate longitude point + if not np.isnan(lon_pt) and (lon_pt < 0.0 or lon_pt > 2.0 * math.pi): + raise Exception("Longitude point out of range") + + # Check for pole points + is_pole_point = False + if np.isnan(lon_pt): + if any_close_lat(lat_pt, ERROR_TOLERANCE): + is_pole_point = True + + if is_pole_point: + # Update latitude for pole points + if isclose(lat_pt, 0.5 * math.pi, ERROR_TOLERANCE): + latlon_box[0, 1] = 0.5 * math.pi # Update max_lat for North Pole + elif isclose(lat_pt, -0.5 * math.pi, ERROR_TOLERANCE): + latlon_box[0, 0] = -0.5 * math.pi # Update min_lat for South Pole + else: + # Update longitude range based on periodicity + if not is_lon_periodic: + # Non-periodic: straightforward min and max updates + if lon_pt < latlon_box[1, 0]: + latlon_box[1, 0] = lon_pt + if lon_pt > latlon_box[1, 1]: + latlon_box[1, 1] = lon_pt + else: + # Periodic longitude handling + # Determine if the new point extends the current longitude range + # considering the periodic boundary at 2*pi + left_lon = latlon_box[1, 0] + right_lon = latlon_box[1, 1] + + # Check if the current box wraps around + wraps_around = left_lon > right_lon + + if wraps_around: + # If the box wraps around, check if the new point is outside the current range + if not (left_lon <= lon_pt or lon_pt <= right_lon): + # Decide to extend either the left or the right + # Calculate the new width for both possibilities + # Option 1: Extend the left boundary to lon_pt + box_a = np.copy(latlon_box) + box_a[1, 0] = lon_pt + d_width_a = _get_latlonbox_width(box_a) + + # Option 2: Extend the right boundary to lon_pt + box_b = np.copy(latlon_box) + box_b[1, 1] = lon_pt + d_width_b = _get_latlonbox_width(box_b) + + # Ensure widths are non-negative + if (d_width_a < 0.0) or (d_width_b < 0.0): + raise Exception( + "Logic error in longitude box width calculation" + ) + + # Choose the box with the smaller width + if d_width_a < d_width_b: + latlon_box = box_a + else: + latlon_box = box_b + else: + # If the box does not wrap around, simply update min or max longitude + if lon_pt < left_lon or lon_pt > right_lon: + # Calculate the new width for both possibilities + # Option 1: Extend the left boundary to lon_pt + box_a = np.copy(latlon_box) + box_a[1, 0] = lon_pt + d_width_a = _get_latlonbox_width(box_a) + + # Option 2: Extend the right boundary to lon_pt + box_b = np.copy(latlon_box) + box_b[1, 1] = lon_pt + d_width_b = _get_latlonbox_width(box_b) + + # Ensure widths are non-negative + if (d_width_a < 0.0) or (d_width_b < 0.0): + raise Exception( + "Logic error in longitude box width calculation" + ) + + # Choose the box with the smaller width + if d_width_a < d_width_b: + latlon_box = box_a + else: + latlon_box = box_b return latlon_box +@njit(cache=True) def _populate_face_latlon_bound( - face_edges_cartesian, - face_edges_lonlat_rad, - is_latlonface: bool = False, + face_edges_xyz, + face_edges_lonlat, + is_latlonface=False, is_GCA_list=None, ): - """Populates the bounding box for each face in the grid by evaluating the - geographical bounds based on the Cartesian and latitudinal/longitudinal - edge connectivity. This function also considers the presence of pole points - within the face's bounds and adjusts the bounding box accordingly. - - Parameters - ---------- - face_edges_cartesian : np.ndarray, shape (n_edges, 2, 3) - An array holding the Cartesian coordinates for the edges of a face, where `n_edges` - is the number of edges for the specific face. Each edge is represented by two points - (start and end), and each point is a 3D vector (x, y, z) in Cartesian coordinates. - - face_edges_lonlat_connectivity_rad : np.ndarray, shape (n_edges, 2, 2) - An array holding the longitude and latitude in radians for the edges of a face, - formatted similarly to `face_edges_cartesian`. Each edge's start and - end points are represented by their longitude and latitude values in radians. - - is_latlonface : bool, optional - A flag indicating if the current face is a latitudinal/longitudinal (latlon) face, - meaning its edges align with lines of constant latitude or longitude. If `True`, - edges are treated as following constant latitudinal or longitudinal lines. If `False`, - edges are considered as great circle arcs (GCA). Default is `False`. - - is_GCA_list : list or np.ndarray, optional - A list or an array of boolean values corresponding to each edge within the face, - indicating whether the edge is a GCA. If `None`, the determination of whether an - edge is a GCA or a constant latitudinal/longitudinal line is based on `is_latlonface`. - Default is `None`. - - Returns - ------- - face_latlon_array : np.ndarray, shape (2, 2) - An array representing the bounding box for the face in latitude and longitude - coordinates (in radians). The first row contains the minimum and maximum latitude - values, while the second row contains the minimum and maximum longitude values. - - Notes - ----- - This function evaluates the presence of North or South pole points within the face's - bounds by inspecting the Cartesian coordinates of the face's edges. It then constructs - the face's bounding box by considering the extreme latitude and longitude values found - among the face's edges, adjusting for the presence of pole points as necessary. - - The bounding box is used to determine the face's geographical extent and is crucial - for spatial analyses involving the grid. - - Example - ------- - Assuming the existence of a grid face with edges defined in both Cartesian and - latitudinal/longitudinal coordinates: - - face_edges_cartesian = np.array([...]) # Cartesian coords - face_edges_connectivity_rad = np.array([...]) # Lon/Lat coords in radians - - Populate the bounding box for the face, treating it as a latlon face: - - face_latlon_bound = _populate_face_latlon_bound(face_edges_cartesian, - face_edges_lonlat_connectivity_rad, - is_latlonface=True) - """ - # Check if face_edges contains pole points - has_north_pole = _pole_point_inside_polygon("North", face_edges_cartesian) - has_south_pole = _pole_point_inside_polygon("South", face_edges_cartesian) + has_north_pole = pole_point_inside_polygon(1, face_edges_xyz, face_edges_lonlat) + has_south_pole = pole_point_inside_polygon(-1, face_edges_xyz, face_edges_lonlat) - face_latlon_array = np.full((2, 2), INT_FILL_VALUE, dtype=np.float64) + # Initialize face_latlon_array with INT_FILL_VALUE + face_latlon_array = np.full((2, 2), np.nan, dtype=np.float64) if has_north_pole or has_south_pole: # Initial assumption that the pole point is inside the face is_center_pole = True - # Define the pole point based on the hemisphere - pole_point = ( - np.array([0.0, 0.0, 1.0]) if has_north_pole else np.array([0.0, 0.0, -1.0]) - ) - # Pre-defined new point latitude based on the pole - new_pt_latlon = np.array( - [np.pi / 2 if has_north_pole else -np.pi / 2, INT_FILL_VALUE], - dtype=np.float64, - ) + pole_point_xyz = np.zeros(3, dtype=np.float64) + pole_point_lonlat = np.zeros(2, dtype=np.float64) + new_pt_latlon = np.zeros(2, dtype=np.float64) - for i in range(face_edges_cartesian.shape[0]): - edge_cart = face_edges_cartesian[i] - edge_lonlat = face_edges_lonlat_rad[i] + if has_north_pole: + pole_point_xyz[2] = 1.0 # [0.0, 0.0, 1.0] + pole_point_lonlat[1] = math.pi / 2 # [0.0, pi/2] + new_pt_latlon[0] = math.pi / 2 # [pi/2, INT_FILL_VALUE] + new_pt_latlon[1] = np.nan + else: + pole_point_xyz[2] = -1.0 # [0.0, 0.0, -1.0] + pole_point_lonlat[1] = -math.pi / 2 # [0.0, -pi/2] + new_pt_latlon[0] = -math.pi / 2 # [-pi/2, INT_FILL_VALUE] + new_pt_latlon[1] = np.nan - # Skip processing if the edge_cart is marked as a dummy with a fill value - if np.any(edge_cart == INT_FILL_VALUE): - continue + for i in range(face_edges_xyz.shape[0]): + edge_xyz = face_edges_xyz[i] + edge_lonlat = face_edges_lonlat[i] - # Extract cartesian coordinates of the edge_cart's endpoints - n1_cart, n2_cart = edge_cart - n1_lonlat, n2_lonlat = edge_lonlat + # Skip processing if the edge is marked as a dummy with a fill value + if np.any(edge_xyz == INT_FILL_VALUE): + continue - # Convert latitudes and longitudes of the nodes to radians - node1_lon_rad, node1_lat_rad = n1_lonlat + # Extract cartesian coordinates of the edge's endpoints + n1_cart = edge_xyz[0] + n2_cart = edge_xyz[1] + n1_lonlat = edge_lonlat[0] + n2_lonlat = edge_lonlat[1] - # Determine if the edge_cart's extreme latitudes need to be considered using the corrected logic - is_GCA = ( - is_GCA_list[i] - if is_GCA_list is not None - else not is_latlonface or n1_cart[2] != n2_cart[2] - ) + # Extract latitudes and longitudes of the nodes + node1_lon_rad = n1_lonlat[0] + node1_lat_rad = n1_lonlat[1] - # Check if the node matches the pole point or if the pole point is within the edge_cart - if allclose(n1_cart, pole_point, atol=ERROR_TOLERANCE) or point_within_gca( - pole_point, np.array([n1_cart, n2_cart]), is_directed=False + # Determine if the edge's extreme latitudes need to be considered + if is_GCA_list is not None: + is_GCA = is_GCA_list[i] + else: + is_GCA = not is_latlonface or n1_cart[2] != n2_cart[2] + + # Check if the node matches the pole point or if the pole point is within the edge + max_abs_diff = np.max(np.abs(n1_cart - pole_point_xyz)) + if max_abs_diff <= ERROR_TOLERANCE or point_within_gca( + pole_point_xyz, + pole_point_lonlat, + n1_cart, + n1_lonlat, + n2_cart, + n2_lonlat, + is_directed=False, ): is_center_pole = False - face_latlon_array = _insert_pt_in_latlonbox( + face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, new_pt_latlon ) # Insert the current node's lat/lon into the latlonbox - face_latlon_array = _insert_pt_in_latlonbox( + face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([node1_lat_rad, node1_lon_rad]) ) + # Create n1n2_cart and n1n2_lonlat arrays using preallocation + n1n2_cart = np.empty((2, 3), dtype=np.float64) + n1n2_cart[0, :] = n1_cart + n1n2_cart[1, :] = n2_cart + + n1n2_lonlat = np.empty((2, 2), dtype=np.float64) + n1n2_lonlat[0, :] = n1_lonlat + n1n2_lonlat[1, :] = n2_lonlat + # Determine extreme latitudes for GCA edges - lat_max, lat_min = ( - ( - extreme_gca_latitude(np.array([n1_cart, n2_cart]), "max"), - extreme_gca_latitude(np.array([n1_cart, n2_cart]), "min"), - ) - if is_GCA - else (node1_lat_rad, node1_lat_rad) - ) + if is_GCA: + lat_max = extreme_gca_latitude(n1n2_cart, n1n2_lonlat, "max") + lat_min = extreme_gca_latitude(n1n2_cart, n1n2_lonlat, "min") + else: + lat_max = node1_lat_rad + lat_min = node1_lat_rad # Insert latitudinal extremes based on pole presence if has_north_pole: - face_latlon_array = _insert_pt_in_latlonbox( + face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([lat_min, node1_lon_rad]) ) - face_latlon_array[0][1] = ( - np.pi / 2 - ) # Ensure north pole is the upper latitude bound + face_latlon_array[0, 1] = math.pi / 2 # Upper latitude bound else: - face_latlon_array = _insert_pt_in_latlonbox( + face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([lat_max, node1_lon_rad]) ) - face_latlon_array[0][0] = ( - -np.pi / 2 - ) # Ensure south pole is the lower latitude bound + face_latlon_array[0, 0] = -math.pi / 2 # Lower latitude bound # Adjust longitude bounds globally if the pole is centrally inside the polygon if is_center_pole: - face_latlon_array[1] = [0.0, 2 * np.pi] + face_latlon_array[1, 0] = 0.0 + face_latlon_array[1, 1] = 2 * math.pi else: # Normal Face - # Iterate through each edge_cart of a face to update the bounding box (latlonbox) with extreme latitudes and longitudes - for i in range(face_edges_cartesian.shape[0]): - edge_cart = face_edges_cartesian[i] - edge_lonlat = face_edges_lonlat_rad[i] + for i in range(face_edges_xyz.shape[0]): + edge_xyz = face_edges_xyz[i] + edge_lonlat = face_edges_lonlat[i] - # Skip processing if the edge_cart is marked as a dummy with a fill value - if np.any(edge_cart == INT_FILL_VALUE): + # Skip processing if the edge is marked as a dummy with a fill value + if np.any(edge_xyz == INT_FILL_VALUE): continue - # Extract cartesian coordinates of the edge_cart's endpoints - n1_cart, n2_cart = edge_cart - n1_lonlat, n2_lonlat = edge_lonlat + # Extract cartesian coordinates of the edge's endpoints + n1_cart = edge_xyz[0] + n2_cart = edge_xyz[1] + n1_lonlat = edge_lonlat[0] + n2_lonlat = edge_lonlat[1] + + # Extract latitudes and longitudes of the nodes + node1_lon_rad = n1_lonlat[0] + node1_lat_rad = n1_lonlat[1] + # node2_lon_rad = n2_lonlat[0] + node2_lat_rad = n2_lonlat[1] + + # Determine if the edge's extreme latitudes need to be considered + if is_GCA_list is not None: + is_GCA = is_GCA_list[i] + else: + is_GCA = not is_latlonface or n1_cart[2] != n2_cart[2] - # Convert latitudes and longitudes of the nodes to radians - node1_lon_rad, node1_lat_rad = n1_lonlat - node2_lon_rad, node2_lat_rad = n2_lonlat + # Create n1n2_cart and n1n2_lonlat arrays using preallocation + n1n2_cart = np.empty((2, 3), dtype=np.float64) + n1n2_cart[0, :] = n1_cart + n1n2_cart[1, :] = n2_cart - # Determine if the edge_cart's extreme latitudes need to be considered using the corrected logic - is_GCA = ( - is_GCA_list[i] - if is_GCA_list is not None - else not is_latlonface or n1_cart[2] != n2_cart[2] - ) + n1n2_lonlat = np.empty((2, 2), dtype=np.float64) + n1n2_lonlat[0, :] = n1_lonlat + n1n2_lonlat[1, :] = n2_lonlat - lat_max, lat_min = ( - ( - extreme_gca_latitude(np.array([n1_cart, n2_cart]), "max"), - extreme_gca_latitude(np.array([n1_cart, n2_cart]), "min"), - ) - if is_GCA - else (node1_lat_rad, node1_lat_rad) - ) + if is_GCA: + lat_max = extreme_gca_latitude(n1n2_cart, n1n2_lonlat, "max") + lat_min = extreme_gca_latitude(n1n2_cart, n1n2_lonlat, "min") + else: + lat_max = node1_lat_rad + lat_min = node1_lat_rad - # Insert extreme latitude points into the latlonbox if they differ from the node latitudes - if not isclose( - node1_lat_rad, lat_max, atol=ERROR_TOLERANCE - ) and not isclose(node2_lat_rad, lat_max, atol=ERROR_TOLERANCE): - # Insert the maximum latitude - face_latlon_array = _insert_pt_in_latlonbox( + # Insert extreme latitude points into the latlonbox + if ( + abs(node1_lat_rad - lat_max) > ERROR_TOLERANCE + and abs(node2_lat_rad - lat_max) > ERROR_TOLERANCE + ): + face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([lat_max, node1_lon_rad]) ) - elif not isclose( - node1_lat_rad, lat_min, atol=ERROR_TOLERANCE - ) and not isclose(node2_lat_rad, lat_min, atol=ERROR_TOLERANCE): - # Insert the minimum latitude - face_latlon_array = _insert_pt_in_latlonbox( + elif ( + abs(node1_lat_rad - lat_min) > ERROR_TOLERANCE + and abs(node2_lat_rad - lat_min) > ERROR_TOLERANCE + ): + face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([lat_min, node1_lon_rad]) ) else: - # Insert the node's latitude and longitude as it matches the extreme latitudes - face_latlon_array = _insert_pt_in_latlonbox( + face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([node1_lat_rad, node1_lon_rad]) ) return face_latlon_array +@njit(cache=True, parallel=True) +def compute_temp_latlon_array( + face_node_connectivity, + faces_edges_cartesian, + faces_edges_lonlat_rad, + n_nodes_per_face, + is_latlonface, + is_face_GCA_list, + INT_FILL_VALUE, +): + n_face = face_node_connectivity.shape[0] + temp_latlon_array = np.full((n_face, 2, 2), INT_FILL_VALUE, dtype=np.float64) + for face_idx in prange(n_face): + cur_face_edges_cartesian = faces_edges_cartesian[ + face_idx, 0 : n_nodes_per_face[face_idx] + ] + cur_faces_edges_lonlat_rad = faces_edges_lonlat_rad[ + face_idx, 0 : n_nodes_per_face[face_idx] + ] + if is_face_GCA_list is not None: + is_GCA_list = is_face_GCA_list[face_idx] + else: + is_GCA_list = None + + temp_latlon_array[face_idx] = _populate_face_latlon_bound( + cur_face_edges_cartesian, + cur_faces_edges_lonlat_rad, + is_latlonface=is_latlonface, + is_GCA_list=is_GCA_list, + ) + return temp_latlon_array + + def _populate_bounds( grid, is_latlonface: bool = False, is_face_GCA_list=None, return_array=False ): @@ -1245,13 +1409,11 @@ def _populate_bounds( This will calculate and store the bounds for each face within the grid, adjusting for any special conditions such as crossing the antimeridian, and return them as a DataArray. """ - temp_latlon_array = np.full((grid.n_face, 2, 2), INT_FILL_VALUE, dtype=np.float64) - # Because Pandas.IntervalIndex does not support naming for each interval, we need to create a mapping - # between the intervals and the face indices - intervals_tuple_list = [] - intervals_name_list = [] + # Ensure grid's cartesian coordinates are normalized + grid.normalize_cartesian_coordinates() + # Prepare data for Numba functions faces_edges_cartesian = _get_cartesian_face_edge_nodes( grid.face_node_connectivity.values, grid.n_face, @@ -1269,42 +1431,29 @@ def _populate_bounds( grid.node_lat.values, ) - for face_idx, face_nodes in enumerate(grid.face_node_connectivity): - face_edges_cartesian = faces_edges_cartesian[face_idx] - - # Remove the edge in the face that contains the fill value - face_edges_cartesian = face_edges_cartesian[ - np.all(face_edges_cartesian != INT_FILL_VALUE, axis=(1, 2)) - ] - - face_edges_lonlat_rad = faces_edges_lonlat_rad[face_idx] - - # Remove the edge in the face that contains the fill value - face_edges_lonlat_rad = face_edges_lonlat_rad[ - np.all(face_edges_lonlat_rad != INT_FILL_VALUE, axis=(1, 2)) - ] - - is_GCA_list = ( - is_face_GCA_list[face_idx] if is_face_GCA_list is not None else None - ) + n_nodes_per_face = grid.n_nodes_per_face.values - temp_latlon_array[face_idx] = _populate_face_latlon_bound( - face_edges_cartesian, - face_edges_lonlat_rad, - is_latlonface=is_latlonface, - is_GCA_list=is_GCA_list, - ) + # Call the Numba-compiled function + temp_latlon_array = compute_temp_latlon_array( + grid.face_node_connectivity.values, + faces_edges_cartesian, + faces_edges_lonlat_rad, + n_nodes_per_face, + is_latlonface, + is_face_GCA_list, + INT_FILL_VALUE, + ) + # Process results outside Numba + intervals_tuple_list = [] + intervals_name_list = [] + for face_idx in range(grid.n_face): assert temp_latlon_array[face_idx][0][0] != temp_latlon_array[face_idx][0][1] assert temp_latlon_array[face_idx][1][0] != temp_latlon_array[face_idx][1][1] lat_array = temp_latlon_array[face_idx][0] - - # Now store the latitude intervals in the tuples intervals_tuple_list.append((lat_array[0], lat_array[1])) intervals_name_list.append(face_idx) - # Because Pandas.IntervalIndex does not support naming for each interval, we need to create a mapping - # between the intervals and the face indices intervalsIndex = pd.IntervalIndex.from_tuples(intervals_tuple_list, closed="both") df_intervals_map = pd.DataFrame( index=intervalsIndex, data=intervals_name_list, columns=["face_id"] @@ -1312,7 +1461,7 @@ def _populate_bounds( bounds = xr.DataArray( temp_latlon_array, - dims=["n_face", "Two", "Two"], + dims=["n_face", "lon_lat", "min_max"], attrs={ "cf_role": "face_latlon_bounds", "_FillValue": INT_FILL_VALUE, diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 88a755b0c..8e0a95277 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -62,6 +62,7 @@ _grid_to_matplotlib_linecollection, _populate_bounds, _construct_boundary_edge_indices, + compute_temp_latlon_array, ) from uxarray.grid.neighbors import ( @@ -75,6 +76,7 @@ constant_lat_intersections_no_extreme, constant_lon_intersections_no_extreme, constant_lat_intersections_face_bounds, + constant_lon_intersections_face_bounds, ) from spatialpandas import GeoDataFrame @@ -1360,15 +1362,19 @@ def face_areas(self, value): @property def bounds(self): - """Latitude Longitude Bounds for each Face in degrees. + """Latitude Longitude Bounds for each Face in radians. Dimensions ``(n_face", two, two)`` """ if "bounds" not in self._ds: - warn( - "Computing 'Grid.bounds' for the first time. This may take some time...", - UserWarning, - ) + if hasattr(compute_temp_latlon_array, "inspect_llvm"): + if len(compute_temp_latlon_array.inspect_llvm()) == 0: + warn( + "Necessary functions for computing face bounds are not translated yet with Numba. This initial" + "translation may take some time.", + RuntimeWarning, + ) + _populate_bounds(self) return self._ds["bounds"] @@ -1379,6 +1385,38 @@ def bounds(self, value): assert isinstance(value, xr.DataArray) self._ds["bounds"] = value + @property + def face_bounds_lon(self): + """Longitude bounds for each face in degrees.""" + + if "face_bounds_lon" not in self._ds: + bounds = self.bounds.values + + bounds_deg = np.rad2deg(bounds[:, 1, :]) + bounds_normalized = (bounds_deg + 180.0) % 360.0 - 180.0 + bounds_lon = bounds_normalized + mask_zero = (bounds_lon[:, 0] == 0) & (bounds_lon[:, 1] == 0) + # for faces that span all longitudes (i.e. pole faces) + bounds_lon[mask_zero] = [-180.0, 180.0] + self._ds["face_bounds_lon"] = xr.DataArray( + data=bounds_lon, + dims=["n_face", "min_max"], + ) + + return self._ds["face_bounds_lon"] + + @property + def face_bounds_lat(self): + """Latitude bounds for each face in degrees.""" + if "face_bounds_lat" not in self._ds: + bounds = self.bounds.values + bounds_lat = np.sort(np.rad2deg(bounds[:, 0, :]), axis=-1) + self._ds["face_bounds_lat"] = xr.DataArray( + data=bounds_lat, + dims=["n_face", "min_max"], + ) + return self._ds["face_bounds_lat"] + @property def face_jacobian(self): """Declare face_jacobian as a property.""" @@ -2189,18 +2227,16 @@ def isel(self, **dim_kwargs): "Indexing must be along a grid dimension: ('n_node', 'n_edge', 'n_face')" ) - def get_edges_at_constant_latitude(self, lat, use_spherical_bounding_box=False): + def get_edges_at_constant_latitude(self, lat: float, use_face_bounds: bool = False): """Identifies the indices of edges that intersect with a line of constant latitude. Parameters ---------- - lon : float - The latitude at which to identify intersecting edges, in degrees. - use_spherical_bounding_box : bool, optional - If `True`, - computes the bounding box for each face using great circle arcs for edges - and considers extreme minimums or maximums to increase accuracy. - Defaults to `False`. + lat : float + The latitude at which to extract the cross-section, in degrees. + Must be between -90.0 and 90.0 + use_face_bounds : bool, optional + If True, uses the bounds of each face for computing intersections. Returns ------- @@ -2213,7 +2249,7 @@ def get_edges_at_constant_latitude(self, lat, use_spherical_bounding_box=False): f"Latitude must be between -90 and 90 degrees. Received {lat}" ) - if use_spherical_bounding_box: + if use_face_bounds: raise NotImplementedError( "Computing the intersection using the spherical bounding box" "is not yet supported." @@ -2225,23 +2261,18 @@ def get_edges_at_constant_latitude(self, lat, use_spherical_bounding_box=False): return edges.squeeze() - def get_faces_at_constant_latitude(self, lat, use_spherical_bounding_box=False): + def get_faces_at_constant_latitude( + self, + lat: float, + ): """ Identifies the indices of faces that intersect with a line of constant latitude. - When `use_spherical_bounding_box` is set to `True`, - the bounding box for each face is computed by representing each edge as a great circle arc. - This approach takes into account the extreme minimums or maximums along the arcs. - Parameters ---------- lat : float - The latitude at which to identify intersecting faces, in degrees. - use_spherical_bounding_box : bool, optional - If `True`, - computes the bounding box for each face using great circle arcs for edges - and considers extreme minimums or maximums to increase accuracy. - Defaults to `False`. + The latitude at which to extract the cross-section, in degrees. + Must be between -90.0 and 90.0 Returns ------- @@ -2254,32 +2285,25 @@ def get_faces_at_constant_latitude(self, lat, use_spherical_bounding_box=False): f"Latitude must be between -90 and 90 degrees. Received {lat}" ) - if use_spherical_bounding_box: - faces = constant_lat_intersections_face_bounds( - lat=lat, - face_min_lat_rad=self.bounds.values[:, 0, 0], - face_max_lat_rad=self.bounds.values[:, 0, 1], - ) - return faces - else: - edges = self.get_edges_at_constant_latitude(lat, use_spherical_bounding_box) - faces = np.unique(self.edge_face_connectivity[edges].data.ravel()) - - return faces[faces != INT_FILL_VALUE] + faces = constant_lat_intersections_face_bounds( + lat=lat, + face_bounds_lat=self.face_bounds_lat.values, + ) + return faces - def get_edges_at_constant_longitude(self, lon, use_spherical_bounding_box=False): + def get_edges_at_constant_longitude( + self, lon: float, use_face_bounds: bool = False + ): """ Identifies the indices of edges that intersect with a line of constant longitude. Parameters ---------- lon : float - The longitude at which to identify intersecting edges, in degrees. - use_spherical_bounding_box : bool, optional - If `True`, - computes the bounding box for each face using great circle arcs for edges - and considers extreme minimums or maximums to increase accuracy. - Defaults to `False`. + The longitude at which to extract the cross-section, in degrees. + Must be between -90.0 and 90.0 + use_face_bounds : bool, optional + If True, uses the bounds of each face for computing intersections. Returns ------- @@ -2292,7 +2316,7 @@ def get_edges_at_constant_longitude(self, lon, use_spherical_bounding_box=False) f"Longitude must be between -180 and 180 degrees. Received {lon}" ) - if use_spherical_bounding_box: + if use_face_bounds: raise NotImplementedError( "Computing the intersection using the spherical bounding box" "is not yet supported." @@ -2303,23 +2327,15 @@ def get_edges_at_constant_longitude(self, lon, use_spherical_bounding_box=False) ) return edges.squeeze() - def get_faces_at_constant_longitude(self, lon, use_spherical_bounding_box=False): + def get_faces_at_constant_longitude(self, lon: float): """ Identifies the indices of faces that intersect with a line of constant longitude. - When `use_spherical_bounding_box` is set to `True`, - the bounding box for each face is computed by representing each edge as a great circle arc. - This approach takes into account the extreme minimums or maximums along the arcs. - Parameters ---------- lon : float - The longitude at which to identify intersecting faces, in degrees. - use_spherical_bounding_box : bool, optional - If `True`, - computes the bounding box for each face using great circle arcs for edges - and considers extreme minimums or maximums to increase accuracy. - Defaults to `False`. + The longitude at which to extract the cross-section, in degrees. + Must be between -90.0 and 90.0 Returns ------- @@ -2327,15 +2343,10 @@ def get_faces_at_constant_longitude(self, lon, use_spherical_bounding_box=False) An array of face indices that intersect with the specified longitude. """ - if use_spherical_bounding_box: - raise NotImplementedError( - "Computing the intersection using the spherical bounding box is not" - "yet supported." - ) - else: - edges = self.get_edges_at_constant_longitude( - lon, use_spherical_bounding_box + if lon > 180.0 or lon < -180.0: + raise ValueError( + f"Longitude must be between -180 and 180 degrees. Received {lon}" ) - faces = np.unique(self.edge_face_connectivity[edges].data.ravel()) - return faces[faces != INT_FILL_VALUE] + faces = constant_lon_intersections_face_bounds(lon, self.face_bounds_lon.values) + return faces diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 147a4b156..528304d35 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -1,10 +1,16 @@ import numpy as np -from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE +from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE, INT_DTYPE from uxarray.grid.utils import _newton_raphson_solver_for_gca_constLat -from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between +from uxarray.grid.arcs import ( + point_within_gca, + in_between, + _extreme_gca_latitude_cartesian, + _point_within_gca_cartesian, +) +from uxarray.grid.coordinates import _xyz_to_lonlat_rad_scalar import platform import warnings -from uxarray.utils.computing import cross_fma, allclose, dot, cross, norm +from uxarray.utils.computing import cross_fma, allclose, cross, norm from numba import njit, prange @@ -122,8 +128,8 @@ def constant_lon_intersections_no_extreme(lon, edge_node_x, edge_node_y, n_edge) return np.unique(intersecting_edges) -@njit -def constant_lat_intersections_face_bounds(lat, face_min_lat_rad, face_max_lat_rad): +@njit(cache=True) +def constant_lat_intersections_face_bounds(lat, face_bounds_lat): """Identifies the candidate faces on a grid that intersect with a given constant latitude. @@ -136,24 +142,24 @@ def constant_lat_intersections_face_bounds(lat, face_min_lat_rad, face_max_lat_r ---------- lat : float The latitude in degrees for which to find intersecting faces. - face_min_lat_rad : numpy.ndarray - A 1D array containing the minimum latitude bounds (in radians) of each face. - face_max_lat_rad : numpy.ndarray - A 1D array containing the maximum latitude bounds (in radians) of each face. + TODO: Returns ------- candidate_faces : numpy.ndarray A 1D array containing the indices of the faces that intersect with the given latitude. """ - lat = np.deg2rad(lat) - within_bounds = (face_min_lat_rad <= lat) & (face_max_lat_rad >= lat) + + face_bounds_lat_min = face_bounds_lat[:, 0] + face_bounds_lat_max = face_bounds_lat[:, 1] + + within_bounds = (face_bounds_lat_min <= lat) & (face_bounds_lat_max >= lat) candidate_faces = np.where(within_bounds)[0] return candidate_faces @njit(cache=True) -def constant_lon_intersections_face_bounds(lon, face_min_lon_rad, face_max_lon_rad): +def constant_lon_intersections_face_bounds(lon, face_bounds_lon): """Identifies the candidate faces on a grid that intersect with a given constant longitude. @@ -166,10 +172,7 @@ def constant_lon_intersections_face_bounds(lon, face_min_lon_rad, face_max_lon_r ---------- lon : float The longitude in degrees for which to find intersecting faces. - face_min_lon_rad : numpy.ndarray - A 1D array containing the minimum longitude bounds (in radians) of each face. - face_max_lon_rad : numpy.ndarray - A 1D array containing the maximum longitude bounds (in radians) of each face. + TODO: Returns ------- @@ -177,121 +180,113 @@ def constant_lon_intersections_face_bounds(lon, face_min_lon_rad, face_max_lon_r A 1D array containing the indices of the faces that intersect with the given longitude. """ - # lon = np.deg2rad(lon) - # lon = (lon + 2 * np.pi) % (2 * np.pi) - # within_bounds = (face_min_lon_rad <= lon) & (face_max_lon_rad >= lon) - # candidate_faces = np.where(within_bounds)[0] - # return candidate_faces - - raise NotImplementedError - - -def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): - """Calculate the intersection point(s) of two Great Circle Arcs (GCAs) in a - Cartesian coordinate system. - - To reduce relative errors, the Fused Multiply-Add (FMA) operation is utilized. - A warning is raised if the given coordinates are not in the cartesian coordinates, or - they cannot be accurately handled using floating-point arithmetic. - - Parameters - ---------- - gca1_cart : [n, 3] np.ndarray where n is the number of intersection points - Cartesian coordinates of the first GCA. - gca2_cart : [n, 3] np.ndarray where n is the number of intersection points - Cartesian coordinates of the second GCA. - fma_disabled : bool, optional (default=True) - If True, the FMA operation is disabled. And a naive `np.cross` is used instead. - - Returns - ------- - np.ndarray - Cartesian coordinates of the intersection point(s). Returns an empty array if no intersections, - a 2D array with one row if one intersection, and a 2D array with two rows if two intersections. + face_bounds_lon_min = face_bounds_lon[:, 0] + face_bounds_lon_max = face_bounds_lon[:, 1] + n_face = face_bounds_lon.shape[0] + + candidate_faces = [] + for i in range(n_face): + cur_face_bounds_lon_min = face_bounds_lon_min[i] + cur_face_bounds_lon_max = face_bounds_lon_max[i] + + if cur_face_bounds_lon_min < cur_face_bounds_lon_max: + if (lon >= cur_face_bounds_lon_min) & (lon <= cur_face_bounds_lon_max): + candidate_faces.append(i) + else: + # antimeridian case + if (lon >= cur_face_bounds_lon_min) | (lon <= cur_face_bounds_lon_max): + candidate_faces.append(i) + + return np.array(candidate_faces, dtype=INT_DTYPE) + + +def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz): + gca_a_xyz = np.asarray(gca_a_xyz) + gca_b_xyz = np.asarray(gca_b_xyz) + + gca_a_lonlat = np.array( + [ + _xyz_to_lonlat_rad_scalar( + gca_a_xyz[0, 0], gca_a_xyz[0, 1], gca_a_xyz[0, 2] + ), + _xyz_to_lonlat_rad_scalar( + gca_a_xyz[1, 0], gca_a_xyz[1, 1], gca_a_xyz[1, 2] + ), + ] + ) + gca_b_lonlat = np.array( + [ + _xyz_to_lonlat_rad_scalar( + gca_b_xyz[0, 0], gca_b_xyz[0, 1], gca_b_xyz[0, 2] + ), + _xyz_to_lonlat_rad_scalar( + gca_b_xyz[1, 0], gca_b_xyz[1, 1], gca_b_xyz[1, 2] + ), + ] + ) + return gca_gca_intersection(gca_a_xyz, gca_a_lonlat, gca_b_xyz, gca_b_lonlat) - Raises - ------ - ValueError - If the input GCAs are not in the cartesian [x, y, z] format. - """ - # Support lists as an input - gca1_cart = np.asarray(gca1_cart) - gca2_cart = np.asarray(gca2_cart) - # Check if the two GCAs are in the cartesian format (size of three) - if gca1_cart.shape[1] != 3 or gca2_cart.shape[1] != 3: +@njit(cache=True) +def gca_gca_intersection(gca_a_xyz, gca_a_lonlat, gca_b_xyz, gca_b_lonlat): + if gca_a_xyz.shape[1] != 3 or gca_b_xyz.shape[1] != 3: raise ValueError("The two GCAs must be in the cartesian [x, y, z] format") - w0, w1 = gca1_cart - v0, v1 = gca2_cart + # Extract points + w0_xyz = gca_a_xyz[0] + w1_xyz = gca_a_xyz[1] + v0_xyz = gca_b_xyz[0] + v1_xyz = gca_b_xyz[1] - # Compute normals and orthogonal bases using FMA - if fma_disabled: - w0w1_norm = cross(w0, w1) - v0v1_norm = cross(v0, v1) - cross_norms = cross(w0w1_norm, v0v1_norm) - else: - w0w1_norm = cross_fma(w0, w1) - v0v1_norm = cross_fma(v0, v1) - cross_norms = cross_fma(w0w1_norm, v0v1_norm) + w0_lonlat = gca_a_lonlat[0] + w1_lonlat = gca_a_lonlat[1] + v0_lonlat = gca_b_lonlat[0] + v1_lonlat = gca_b_lonlat[1] - # Raise a warning for windows users - if platform.system() == "Windows": - warnings.warn( - "The C/C++ implementation of FMA in MS Windows is reportedly broken. Use with care. (bug report: " - "https://bugs.python.org/msg312480)" - "The single rounding cannot be guaranteed, hence the relative error bound of 3u cannot be guaranteed." - ) + # Compute normals and orthogonal bases + w0w1_norm = cross(w0_xyz, w1_xyz) + v0v1_norm = cross(v0_xyz, v1_xyz) + cross_norms = cross(w0w1_norm, v0v1_norm) - # Check perpendicularity conditions and floating-point arithmetic limitations - if not allclose(dot(w0w1_norm, w0), 0.0, atol=MACHINE_EPSILON) or not allclose( - dot(w0w1_norm, w1), 0.0, atol=MACHINE_EPSILON - ): - warnings.warn( - "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution." - ) + # Initialize result array and counter + res = np.empty((2, 3)) + count = 0 - if not allclose(dot(v0v1_norm, v0), 0.0, 0.0, atol=MACHINE_EPSILON) or not allclose( - dot(v0v1_norm, v1), 0.0, atol=MACHINE_EPSILON - ): - warnings.warn( - "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. " - ) - - if not allclose( - dot(cross_norms, v0v1_norm), 0.0, atol=MACHINE_EPSILON - ) or not allclose(dot(cross_norms, w0w1_norm), 0.0, atol=MACHINE_EPSILON): - warnings.warn( - "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. " - ) - - # If the cross_norms is zero, the two GCAs are parallel + # Check if the two GCAs are parallel if allclose(cross_norms, 0.0, atol=MACHINE_EPSILON): - res = [] - # Check if the two GCAs are overlapping - if point_within_gca(v0, [w0, w1]): - # The two GCAs are overlapping, return both end points - res.append(v0) + if point_within_gca(v0_xyz, v0_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat): + res[count, :] = v0_xyz + count += 1 + + if point_within_gca(v1_xyz, v1_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat): + res[count, :] = v1_xyz + count += 1 - if point_within_gca(v1, [w0, w1]): - res.append(v1) - return np.array(res) + return res[:count, :] # Normalize the cross_norms cross_norms = cross_norms / norm(cross_norms) - x1 = cross_norms - x2 = -x1 + x1_xyz = cross_norms + x2_xyz = -x1_xyz - res = [] + # Get lon/lat arrays + x1_lonlat = _xyz_to_lonlat_rad_scalar(x1_xyz[0], x1_xyz[1], x1_xyz[2]) + x2_lonlat = _xyz_to_lonlat_rad_scalar(x2_xyz[0], x2_xyz[1], x2_xyz[2]) - # Determine which intersection point is within the GCAs range - if point_within_gca(x1, [w0, w1]) and point_within_gca(x1, [v0, v1]): - res.append(x1) + # Check intersection points + if point_within_gca( + x1_xyz, x1_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat + ) and point_within_gca(x1_xyz, x1_lonlat, v0_xyz, v0_lonlat, v1_xyz, v1_lonlat): + res[count, :] = x1_xyz + count += 1 - if point_within_gca(x2, [w0, w1]) and point_within_gca(x2, [v0, v1]): - res.append(x2) + if point_within_gca( + x2_xyz, x2_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat + ) and point_within_gca(x2_xyz, x2_lonlat, v0_xyz, v0_lonlat, v1_xyz, v1_lonlat): + res[count, :] = x2_xyz + count += 1 - return np.array(res) + return res[:count, :] def gca_const_lat_intersection( @@ -351,8 +346,9 @@ def gca_const_lat_intersection( return res # If the constant latitude is not the same as the GCA endpoints, calculate the intersection point - lat_min = extreme_gca_latitude(gca_cart, extreme_type="min") - lat_max = extreme_gca_latitude(gca_cart, extreme_type="max") + # TODO: + lat_min = _extreme_gca_latitude_cartesian(gca_cart, extreme_type="min") + lat_max = _extreme_gca_latitude_cartesian(gca_cart, extreme_type="max") constLat_rad = np.arcsin(constZ) @@ -389,7 +385,7 @@ def gca_const_lat_intersection( res = None # Now test which intersection point is within the GCA range - if point_within_gca(p1, gca_cart, is_directed=is_directed): + if _point_within_gca_cartesian(p1, gca_cart, is_directed=is_directed): try: converged_pt = _newton_raphson_solver_for_gca_constLat( p1, gca_cart, verbose=verbose @@ -411,7 +407,7 @@ def gca_const_lat_intersection( except RuntimeError: raise RuntimeError(f"Error encountered with initial guess: {p1}") - if point_within_gca(p2, gca_cart, is_directed=is_directed): + if _point_within_gca_cartesian(p2, gca_cart, is_directed=is_directed): try: converged_pt = _newton_raphson_solver_for_gca_constLat( p2, gca_cart, verbose=verbose diff --git a/uxarray/plot/accessor.py b/uxarray/plot/accessor.py index 47cc5c9d7..9c0adc9b9 100644 --- a/uxarray/plot/accessor.py +++ b/uxarray/plot/accessor.py @@ -399,7 +399,6 @@ def polygons( if projection is not None: kwargs["projection"] = projection - print("Hello") kwargs["geo"] = True if "crs" not in kwargs: central_longitude = projection.proj4_params["lon_0"] diff --git a/uxarray/utils/computing.py b/uxarray/utils/computing.py index 588d921dd..b3db27569 100644 --- a/uxarray/utils/computing.py +++ b/uxarray/utils/computing.py @@ -4,6 +4,16 @@ from numba import njit +@njit(cache=True) +def clip(a, a_min, a_max): + return np.clip(a, a_min, a_max) + + +@njit(cache=True) +def arcsin(x): + return np.arcsin(x) + + @njit(cache=True) def all(a): """Numba decorated implementation of ``np.all()``