From dff3971a67d7ff03f9af229ab1cd5feacf50537f Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Tue, 23 Apr 2024 19:16:42 -0700 Subject: [PATCH] Zonal face weights at constant lat (#555) * Initial Draft Algorithm * Revert "Initial Draft Algorithm" This reverts commit 99573b8edfc2dc2fc329e91d395da4d907b42e4a. * Initial Draft Algorithm * Initial Draft Algorithm * Updated documentation and comments, added test case, and updated index.rst * updated function name * Refactored function to get cartesian coords * Merged main, fixed conflicts updated doc string * Added test, addressed review comments * fixed pre-commit * Updated docstring and variable names * Initial commit, working in progress * Codes Refractor and Docstring updates * Update some layout * update * Fix the filled value case * changes * changes * changes * First draft of the function * working on the test case * commits * add directed option for the gca constlat intersection * working on it * Change the `is_directed` default as false * Change the `is_directed` default as false * Fix testcase * Finished `_get_zonal_face_weight_rad` * Allow users to specify if the input is GCA or constLat * Consider the equator case * Draft complete * structure refractor needed * Initial commit * Finished * precommit * Fix conflic * Add docstring * Abandon `pd.append` * delete duplicate * Fix testcases * try fixing testcases * try fixing testcases again * try fixing testcases again * try fixing testcases again again * Add error-tolerence in `point_within_gca` * pre-commit fix * Add the `is_latlonface` flag * Add the `is_latlonface` flag * Try to consider about the concave face case * WIP * Update uxarray/grid/integrate.py Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> * working for merge conflict * working for merge conflict * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * WIP * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Resolve TestFaceWeights cases * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Resolve All Testcases * clean up * Add doc * fix doc --------- Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Co-authored-by: aazedwick Co-authored-by: Aaron Zedwick <95507181+aaronzedwick@users.noreply.github.com> Co-authored-by: Orhan Eroglu <32553057+erogluorhan@users.noreply.github.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- docs/internal_api/index.rst | 13 ++ test/test_integrate.py | 414 ++++++++++++++++++++++++++++++++++ test/test_intersections.py | 10 +- uxarray/constants.py | 2 + uxarray/grid/arcs.py | 2 +- uxarray/grid/integrate.py | 352 +++++++++++++++++++++++++++++ uxarray/grid/intersections.py | 7 +- 7 files changed, 788 insertions(+), 12 deletions(-) create mode 100644 uxarray/grid/integrate.py diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst index 13b572854..fc90d681d 100644 --- a/docs/internal_api/index.rst +++ b/docs/internal_api/index.rst @@ -208,6 +208,19 @@ Accurate Computing Utils utils.computing._acc_sqrt utils.computing._split + +Integration +----------- +.. autosummary:: + :toctree: generated/ + + grid.integrate._get_zonal_faces_weight_at_constLat + grid.integrate._get_zonal_face_interval + grid.integrate._process_overlapped_intervals + grid.integrate._is_edge_gca + grid.integrate._get_faces_constLat_intersection_info + + Remapping ========= diff --git a/test/test_integrate.py b/test/test_integrate.py index 50864015c..af26941c6 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -3,10 +3,13 @@ from unittest import TestCase from pathlib import Path import numpy as np +import pandas as pd import numpy.testing as nt import uxarray as ux +from uxarray.grid.coordinates import node_lonlat_rad_to_xyz +from uxarray.grid.integrate import _get_zonal_face_interval, _process_overlapped_intervals, _get_zonal_faces_weight_at_constLat current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -54,3 +57,414 @@ def test_multi_dim(self): assert integral.ndim == len(dims) - 1 nt.assert_almost_equal(integral, np.ones((5, 5)) * 4 * np.pi) + + +class TestFaceWeights(TestCase): + + def test_get_zonal_face_interval(self): + """Test that the zonal face weights are correct.""" + vertices_lonlat = [[1.6 * np.pi, 0.25 * np.pi], + [1.6 * np.pi, -0.25 * np.pi], + [0.4 * np.pi, -0.25 * np.pi], + [0.4 * np.pi, 0.25 * np.pi]] + vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat] + + face_edge_nodes = np.array([[vertices[0], vertices[1]], + [vertices[1], vertices[2]], + [vertices[2], vertices[3]], + [vertices[3], vertices[0]]]) + + constZ = np.sin(0.20) + # The latlon bounds for the latitude is not necessarily correct below since we don't use the latitudes bound anyway + interval_df = _get_zonal_face_interval(face_edge_nodes, constZ, + np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi, + 0.4 * np.pi]]), + is_directed=False) + expected_interval_df = pd.DataFrame({ + 'start': [1.6 * np.pi, 0.0], + 'end': [2.0 * np.pi, 00.4 * np.pi] + }) + # Sort both DataFrames by 'start' column before comparison + expected_interval_df_sorted = expected_interval_df.sort_values(by='start').reset_index(drop=True) + + # Converting the sorted DataFrames to NumPy arrays + actual_values_sorted = interval_df[['start', 'end']].to_numpy() + expected_values_sorted = expected_interval_df_sorted[['start', 'end']].to_numpy() + + # Asserting almost equal arrays + nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) + + def test_get_zonal_face_interval_GCA_constLat(self): + """Test that the zonal face weights are correct.""" + vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi], + [-0.4 * np.pi, -0.25 * np.pi], + [0.4 * np.pi, -0.25 * np.pi], + [0.4 * np.pi, 0.25 * np.pi]] + + vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat] + + face_edge_nodes = np.array([[vertices[0], vertices[1]], + [vertices[1], vertices[2]], + [vertices[2], vertices[3]], + [vertices[3], vertices[0]]]) + + constZ = np.sin(0.20 * np.pi) + interval_df = _get_zonal_face_interval(face_edge_nodes, constZ, + np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi, + 0.4 * np.pi]]), + is_directed=False, is_GCA_list=np.array([True, False, True, False])) + expected_interval_df = pd.DataFrame({ + 'start': [1.6 * np.pi, 0.0], + 'end': [2.0 * np.pi, 00.4 * np.pi] + }) + # Sort both DataFrames by 'start' column before comparison + expected_interval_df_sorted = expected_interval_df.sort_values(by='start').reset_index(drop=True) + + # Converting the sorted DataFrames to NumPy arrays + actual_values_sorted = interval_df[['start', 'end']].to_numpy() + expected_values_sorted = expected_interval_df_sorted[['start', 'end']].to_numpy() + + # Asserting almost equal arrays + nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) + + def test_get_zonal_face_interval_equator(self): + """Test that the zonal face weights are correct.""" + vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi], [-0.4 * np.pi, 0.0], + [0.4 * np.pi, 0.0], [0.4 * np.pi, 0.25 * np.pi]] + + vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat] + + face_edge_nodes = np.array([[vertices[0], vertices[1]], + [vertices[1], vertices[2]], + [vertices[2], vertices[3]], + [vertices[3], vertices[0]]]) + + interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0, + np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi, + 0.4 * np.pi]]), + is_directed=False, is_GCA_list=np.array([True, True, True, True])) + expected_interval_df = pd.DataFrame({ + 'start': [1.6 * np.pi, 0.0], + 'end': [2.0 * np.pi, 00.4 * np.pi] + }) + # Sort both DataFrames by 'start' column before comparison + expected_interval_df_sorted = expected_interval_df.sort_values(by='start').reset_index(drop=True) + + # Converting the sorted DataFrames to NumPy arrays + actual_values_sorted = interval_df[['start', 'end']].to_numpy() + expected_values_sorted = expected_interval_df_sorted[['start', 'end']].to_numpy() + + # Asserting almost equal arrays + nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) + + # Even if we change the is_GCA_list to False, the result should be the same + interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0, + np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi, + 0.4 * np.pi]]), + is_directed=False, is_GCA_list=np.array([True, False, True, False])) + expected_interval_df = pd.DataFrame({ + 'start': [1.6 * np.pi, 0.0], + 'end': [2.0 * np.pi, 00.4 * np.pi] + }) + + # Sort both DataFrames by 'start' column before comparison + expected_interval_df_sorted = expected_interval_df.sort_values(by='start').reset_index(drop=True) + + # Converting the sorted DataFrames to NumPy arrays + actual_values_sorted = interval_df[['start', 'end']].to_numpy() + expected_values_sorted = expected_interval_df_sorted[['start', 'end']].to_numpy() + + # Asserting almost equal arrays + nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) + + def test_process_overlapped_intervals(self): + # Example data that has overlapping intervals and gap + intervals_data = [ + { + 'start': 0.0, + 'end': 100.0, + 'face_index': 0 + }, + { + 'start': 50.0, + 'end': 150.0, + 'face_index': 1 + }, + { + 'start': 140.0, + 'end': 150.0, + 'face_index': 2 + }, + { + 'start': 150.0, + 'end': 250.0, + 'face_index': 3 + }, + { + 'start': 260.0, + 'end': 350.0, + 'face_index': 4 + }, + ] + + df = pd.DataFrame(intervals_data) + df['interval'] = df.apply(lambda row: pd.Interval( + left=row['start'], right=row['end'], closed='both'), + axis=1) + df['interval'] = df['interval'].astype('interval[float64]') + + # Expected result + expected_overlap_contributions = np.array({ + 0: 75.0, + 1: 70.0, + 2: 5.0, + 3: 100.0, + 4: 90.0 + }) + overlap_contributions, total_length = _process_overlapped_intervals(df) + self.assertEqual(total_length, 340.0) + nt.assert_array_equal(overlap_contributions, + expected_overlap_contributions) + + def test_process_overlapped_intervals_antimerdian(self): + intervals_data = [ + { + 'start': 350.0, + 'end': 360.0, + 'face_index': 0 + }, + { + 'start': 0.0, + 'end': 100.0, + 'face_index': 0 + }, + { + 'start': 100.0, + 'end': 150.0, + 'face_index': 1 + }, + { + 'start': 100.0, + 'end': 300.0, + 'face_index': 2 + }, + { + 'start': 310.0, + 'end': 360.0, + 'face_index': 3 + }, + ] + + df = pd.DataFrame(intervals_data) + df['interval'] = df.apply(lambda row: pd.Interval( + left=row['start'], right=row['end'], closed='both'), + axis=1) + df['interval'] = df['interval'].astype('interval[float64]') + + # Expected result + expected_overlap_contributions = np.array({ + 0: 105.0, + 1: 25.0, + 2: 175.0, + 3: 45.0 + }) + overlap_contributions, total_length = _process_overlapped_intervals(df) + self.assertEqual(total_length, 350.0) + nt.assert_array_equal(overlap_contributions, + expected_overlap_contributions) + + def test_get_zonal_faces_weight_at_constLat_equator(self): + face_0 = [[1.7 * np.pi, 0.25 * np.pi], [1.7 * np.pi, 0.0], + [0.3 * np.pi, 0.0], [0.3 * np.pi, 0.25 * np.pi]] + face_1 = [[0.3 * np.pi, 0.0], [0.3 * np.pi, -0.25 * np.pi], + [0.6 * np.pi, -0.25 * np.pi], [0.6 * np.pi, 0.0]] + face_2 = [[0.3 * np.pi, 0.25 * np.pi], [0.3 * np.pi, 0.0], [np.pi, 0.0], + [np.pi, 0.25 * np.pi]] + face_3 = [[0.7 * np.pi, 0.0], [0.7 * np.pi, -0.25 * np.pi], + [np.pi, -0.25 * np.pi], [np.pi, 0.0]] + + # Convert the face vertices to xyz coordinates + face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0] + face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1] + face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2] + face_3 = [node_lonlat_rad_to_xyz(v) for v in face_3] + + face_0_edge_nodes = np.array([[face_0[0], face_0[1]], + [face_0[1], face_0[2]], + [face_0[2], face_0[3]], + [face_0[3], face_0[0]]]) + face_1_edge_nodes = np.array([[face_1[0], face_1[1]], + [face_1[1], face_1[2]], + [face_1[2], face_1[3]], + [face_1[3], face_1[0]]]) + face_2_edge_nodes = np.array([[face_2[0], face_2[1]], + [face_2[1], face_2[2]], + [face_2[2], face_2[3]], + [face_2[3], face_2[0]]]) + face_3_edge_nodes = np.array([[face_3[0], face_3[1]], + [face_3[1], face_3[2]], + [face_3[2], face_3[3]], + [face_3[3], face_3[0]]]) + + face_0_latlon_bound = np.array([[0.0, 0.25 * np.pi], + [1.7 * np.pi, 0.3 * np.pi]]) + face_1_latlon_bound = np.array([[-0.25 * np.pi, 0.0], + [0.3 * np.pi, 0.6 * np.pi]]) + face_2_latlon_bound = np.array([[0.0, 0.25 * np.pi], + [0.3 * np.pi, np.pi]]) + face_3_latlon_bound = np.array([[-0.25 * np.pi, 0.0], + [0.7 * np.pi, np.pi]]) + + latlon_bounds = np.array([ + face_0_latlon_bound, face_1_latlon_bound, face_2_latlon_bound, + face_3_latlon_bound + ]) + + expected_weight_df = pd.DataFrame({ + 'face_index': [0, 1, 2, 3], + 'weight': [0.46153, 0.11538, 0.30769, 0.11538] + }) + + # Assert the results is the same to the 3 decimal places + weight_df = _get_zonal_faces_weight_at_constLat(np.array([ + face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes, + face_3_edge_nodes + ]), + 0.0, + latlon_bounds, + is_directed=False) + + nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3) + + def test_get_zonal_faces_weight_at_constLat_regular(self): + face_0 = [[1.7 * np.pi, 0.25 * np.pi], [1.7 * np.pi, 0.0], + [0.3 * np.pi, 0.0], [0.3 * np.pi, 0.25 * np.pi]] + face_1 = [[0.4 * np.pi, 0.3 * np.pi], [0.4 * np.pi, 0.0], + [0.5 * np.pi, 0.0], [0.5 * np.pi, 0.3 * np.pi]] + face_2 = [[0.5 * np.pi, 0.25 * np.pi], [0.5 * np.pi, 0.0], [np.pi, 0.0], + [np.pi, 0.25 * np.pi]] + face_3 = [[1.2 * np.pi, 0.25 * np.pi], [1.2 * np.pi, 0.0], + [1.6 * np.pi, -0.01 * np.pi], [1.6 * np.pi, 0.25 * np.pi]] + + # Convert the face vertices to xyz coordinates + face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0] + face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1] + face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2] + face_3 = [node_lonlat_rad_to_xyz(v) for v in face_3] + + face_0_edge_nodes = np.array([[face_0[0], face_0[1]], + [face_0[1], face_0[2]], + [face_0[2], face_0[3]], + [face_0[3], face_0[0]]]) + face_1_edge_nodes = np.array([[face_1[0], face_1[1]], + [face_1[1], face_1[2]], + [face_1[2], face_1[3]], + [face_1[3], face_1[0]]]) + face_2_edge_nodes = np.array([[face_2[0], face_2[1]], + [face_2[1], face_2[2]], + [face_2[2], face_2[3]], + [face_2[3], face_2[0]]]) + face_3_edge_nodes = np.array([[face_3[0], face_3[1]], + [face_3[1], face_3[2]], + [face_3[2], face_3[3]], + [face_3[3], face_3[0]]]) + + face_0_latlon_bound = np.array([[0.0, 0.25 * np.pi], + [1.7 * np.pi, 0.3 * np.pi]]) + face_1_latlon_bound = np.array([[0, 0.3 * np.pi], + [0.4 * np.pi, 0.5 * np.pi]]) + face_2_latlon_bound = np.array([[0.0, 0.25 * np.pi], + [0.5 * np.pi, np.pi]]) + face_3_latlon_bound = np.array([[-0.01 * np.pi, 0.25 * np.pi], + [1.2 * np.pi, 1.6 * np.pi]]) + + latlon_bounds = np.array([ + face_0_latlon_bound, face_1_latlon_bound, face_2_latlon_bound, + face_3_latlon_bound + ]) + + expected_weight_df = pd.DataFrame({ + 'face_index': [0, 1, 2, 3], + 'weight': [0.375, 0.0625, 0.3125, 0.25] + }) + + + + # Assert the results is the same to the 3 decimal places + weight_df = _get_zonal_faces_weight_at_constLat(np.array([ + face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes, + face_3_edge_nodes + ]), + np.sin(0.1 * np.pi), + latlon_bounds, + is_directed=False) + + nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3) + + def test_get_zonal_faces_weight_at_constLat_latlonface(self): + face_0 = [[np.deg2rad(350), np.deg2rad(40)], [np.deg2rad(350), np.deg2rad(20)], + [np.deg2rad(10), np.deg2rad(20)], [np.deg2rad(10), np.deg2rad(40)]] + face_1 = [[np.deg2rad(5), np.deg2rad(20)], [np.deg2rad(5), np.deg2rad(10)], + [np.deg2rad(25), np.deg2rad(10)], [np.deg2rad(25), np.deg2rad(20)]] + face_2 = [[np.deg2rad(30), np.deg2rad(40)], [np.deg2rad(30), np.deg2rad(20)], + [np.deg2rad(40), np.deg2rad(20)], [np.deg2rad(40), np.deg2rad(40)]] + + # Convert the face vertices to xyz coordinates + face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0] + face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1] + face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2] + + + face_0_edge_nodes = np.array([[face_0[0], face_0[1]], + [face_0[1], face_0[2]], + [face_0[2], face_0[3]], + [face_0[3], face_0[0]]]) + face_1_edge_nodes = np.array([[face_1[0], face_1[1]], + [face_1[1], face_1[2]], + [face_1[2], face_1[3]], + [face_1[3], face_1[0]]]) + face_2_edge_nodes = np.array([[face_2[0], face_2[1]], + [face_2[1], face_2[2]], + [face_2[2], face_2[3]], + [face_2[3], face_2[0]]]) + + face_0_latlon_bound = np.array([[np.deg2rad(20), np.deg2rad(40)], + [np.deg2rad(350), np.deg2rad(10)]]) + face_1_latlon_bound = np.array([[np.deg2rad(10), np.deg2rad(20)], + [np.deg2rad(5), np.deg2rad(25)]]) + face_2_latlon_bound = np.array([[np.deg2rad(20), np.deg2rad(40)], + [np.deg2rad(30), np.deg2rad(40)]]) + + + latlon_bounds = np.array([ + face_0_latlon_bound, face_1_latlon_bound, face_2_latlon_bound + ]) + + sum = 17.5 + 17.5 + 10 + expected_weight_df = pd.DataFrame({ + 'face_index': [0, 1, 2], + 'weight': [17.5 / sum, 17.5/sum, 10/sum] + }) + + # Assert the results is the same to the 3 decimal places + weight_df = _get_zonal_faces_weight_at_constLat(np.array([ + face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes + ]), + np.sin(np.deg2rad(20)), + latlon_bounds, + is_directed=False, is_latlonface=True) + + + nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3) + + + + # A error will be raise if we don't set is_latlonface=True since the face_2 will be concave if + # It's edges are all GCA + with self.assertRaises(ValueError): + _get_zonal_faces_weight_at_constLat(np.array([ + face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes + ]), + np.deg2rad(20), + latlon_bounds, + is_directed=False) diff --git a/test/test_intersections.py b/test/test_intersections.py index 0faa028dd..f17b3c359 100644 --- a/test/test_intersections.py +++ b/test/test_intersections.py @@ -102,9 +102,7 @@ def test_GCA_constLat_intersections_antimeridian(self): np.deg2rad(10.0)]) ]) - res = gca_constLat_intersection(GCR1_cart, - np.deg2rad(60.0), - verbose=True) + res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(60.0)), verbose=True) res_lonlat_rad = node_xyz_to_lonlat_rad(res[0].tolist()) self.assertTrue( np.allclose(res_lonlat_rad, @@ -119,9 +117,7 @@ def test_GCA_constLat_intersections_empty(self): np.deg2rad(10.0)]) ]) - res = gca_constLat_intersection(GCR1_cart, - np.deg2rad(-10.0), - verbose=False) + res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(-10.0)), verbose=False) self.assertTrue(res.size == 0) def test_GCA_constLat_intersections_two_pts(self): @@ -135,5 +131,5 @@ def test_GCA_constLat_intersections_two_pts(self): query_lat = (np.deg2rad(10.0) + max_lat) / 2.0 - res = gca_constLat_intersection(GCR1_cart, query_lat, verbose=False) + res = gca_constLat_intersection(GCR1_cart, np.sin(query_lat), verbose=False) self.assertTrue(res.shape[0] == 2) diff --git a/uxarray/constants.py b/uxarray/constants.py index 0704db2f7..d42befac4 100644 --- a/uxarray/constants.py +++ b/uxarray/constants.py @@ -4,10 +4,12 @@ INT_DTYPE = np.intp INT_FILL_VALUE = np.iinfo(INT_DTYPE).min + # According to Ogita, Takeshi & Rump, Siegfried & Oishi, Shin’ichi. (2005). Accurate Sum and Dot Product. # SIAM J. Scientific Computing. 26. 1955-1988. 10.1137/030601818. # half of the working precision is already the most optimal value for the error tolerance, # more tailored values will be used in the future. + ERROR_TOLERANCE = 1.0e-8 ENABLE_JIT_CACHE = True diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 86c3e44fd..e2c4bddfb 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -82,7 +82,7 @@ def point_within_gca(pt, gca_cart, is_directed=False): if np.isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): # If the pt and the GCA are on the same longitude (the y coordinates are the same) - if GCRv0_lonlat[0] == pt_lonlat[0]: + if np.isclose(GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): # 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]) else: diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py new file mode 100644 index 000000000..7542fdbbb --- /dev/null +++ b/uxarray/grid/integrate.py @@ -0,0 +1,352 @@ +import numpy as np +from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE +from uxarray.grid.intersections import gca_constLat_intersection +from uxarray.grid.coordinates import node_xyz_to_lonlat_rad +import pandas as pd + +DUMMY_EDGE_VALUE = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE] + + +def _get_zonal_faces_weight_at_constLat( + faces_edges_cart, + latitude_cart, + face_latlon_bound, + is_directed=False, + is_latlonface=False, + is_face_GCA_list=None, +): + """Utilize the sweep line algorithm to calculate the weight of each face at + a constant latitude. + + Parameters + ---------- + face_edges_cart : np.ndarray + A list of face polygon represented by edges in Cartesian coordinates. Shape: (n_faces, n_edges, 2, 3) + + latitude_cart : float + The latitude in Cartesian coordinates (The normalized z coordinate) + + face_latlon_bound : np.ndarray + The list of latitude and longitude bounds of faces. Shape: (n_faces,2, 2), + [...,[lat_min, lat_max], [lon_min, lon_max],...] + + 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. + + is_latlonface : bool, optional, default=False + A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means + That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. + Default is False. + + is_face_GCA_list : np.ndarray, optional (default=None) + A list of boolean values that indicates if the edge in that face is a GCA. Shape: (n_faces,n_edges). + True means edge face is a GCA. + False mean this edge is a constant latitude. + If None, all edges are considered as GCA. This attribute will overwrite the is_latlonface attribute. + + Returns + ------- + weights_df : pandas.DataFrame + A DataFrame with the calculated weights of each face. The DataFrame has two columns: + - 'face_index': The index of the face (integer). + - 'weight': The calculated weight of the face in radian (float). + The DataFrame is indexed by the face indices, providing a mapping from each face to its corresponding weight. + """ + intervals_list = [] + + # Iterate through all faces and their edges + for face_index, face_edges in enumerate(faces_edges_cart): + if is_face_GCA_list is not None: + is_GCA_list = is_face_GCA_list[face_index] + else: + is_GCA_list = None + face_interval_df = _get_zonal_face_interval( + face_edges, + latitude_cart, + face_latlon_bound[face_index], + is_directed=is_directed, + is_latlonface=is_latlonface, + is_GCA_list=is_GCA_list, + ) + for _, row in face_interval_df.iterrows(): + intervals_list.append( + {"start": row["start"], "end": row["end"], "face_index": face_index} + ) + + intervals_df = pd.DataFrame(intervals_list) + overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) + + # Calculate weights for each face + weights = { + face_index: overlap_contributions.get(face_index, 0.0) / total_length + for face_index in range(len(faces_edges_cart)) + } + + # Convert weights to DataFrame + weights_df = pd.DataFrame(list(weights.items()), columns=["face_index", "weight"]) + return weights_df + + +def _is_edge_gca(is_GCA_list, is_latlonface, edges_z): + """Determine if each edge is a Great Circle Arc (GCA) or a constant + latitude line in a vectorized manner. + + Parameters: + ---------- + is_GCA_list : np.ndarray or None + An array indicating whether each edge is a GCA (True) or a constant latitude line (False). + Shape: (n_edges). If None, edge types are determined based on `is_latlonface` and the z-coordinates. + is_latlonface : bool + Flag indicating if all edges should be considered as lat-lon faces, which implies all edges + are either constant latitude or longitude lines. + edges_z : np.ndarray + Array containing the z-coordinates for each vertex of the edges. This is used to determine + whether edges are on the equator or if they are aligned in latitude when `is_GCA_list` is None. + Shape should be (n_edges, 2). + + Returns: + ------- + np.ndarray + A boolean array where each element indicates whether the corresponding edge is considered a GCA. + True for GCA, False for constant latitude line. + """ + if is_GCA_list is not None: + return is_GCA_list + if is_latlonface: + return ~np.isclose(edges_z[:, 0], edges_z[:, 1], atol=ERROR_TOLERANCE) + return ~( + np.isclose(edges_z[:, 0], 0, atol=ERROR_TOLERANCE) + & np.isclose(edges_z[:, 1], 0, atol=ERROR_TOLERANCE) + ) + + +def _get_faces_constLat_intersection_info( + face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed +): + """Processes each edge of a face polygon in a vectorized manner to + determine overlaps and calculate the intersections for a given latitude and + the faces. + + Parameters: + ---------- + face_edges_cart : np.ndarray + A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3). + latitude_cart : float + The latitude in Cartesian coordinates to which intersections or overlaps are calculated. + is_GCA_list : np.ndarray or None + An array indicating whether each edge is a GCA (True) or a constant latitude line (False). + Shape: (n_edges). If None, the function will determine edge types based on `is_latlonface`. + is_latlonface : bool + Flag indicating if all faces are considered as lat-lon faces, meaning all edges are either + constant latitude or longitude lines. This parameter overwrites the `is_GCA_list` if set to True. + is_directed : bool + Flag indicating if the GCA should be considered as directed (from v0 to v1). If False, + the smaller circle (less than 180 degrees) side of the GCA is used. + + Returns: + ------- + tuple + A tuple containing: + - intersections_pts_list_cart (list): A list of intersection points where each point is where an edge intersects with the latitude. + - overlap_flag (bool): A boolean indicating if any overlap with the latitude was detected. + - overlap_edge (np.ndarray or None): The edge that overlaps with the latitude, if any; otherwise, None. + """ + valid_edges_mask = ~(np.any(face_edges_cart == DUMMY_EDGE_VALUE, axis=(1, 2))) + + # Apply mask to filter out dummy edges + valid_edges = face_edges_cart[valid_edges_mask] + + # Extract Z coordinates for edge determination + edges_z = valid_edges[:, :, 2] + + # Determine if each edge is GCA or constant latitude + is_GCA = _is_edge_gca(is_GCA_list, is_latlonface, edges_z) + + # Check overlap with latitude + overlaps_with_latitude = np.all( + np.isclose(edges_z, latitude_cart, atol=ERROR_TOLERANCE), axis=1 + ) + overlap_flag = np.any(overlaps_with_latitude & ~is_GCA) + + # Identify overlap edges if needed + intersections_pts_list_cart = [] + if overlap_flag: + overlap_index = np.where(overlaps_with_latitude & ~is_GCA)[0][0] + intersections_pts_list_cart.extend(valid_edges[overlap_index]) + else: + # Calculate intersections (assuming a batch-capable intersection function) + for idx, edge in enumerate(valid_edges): + if is_GCA[idx]: + intersections = gca_constLat_intersection( + edge, latitude_cart, is_directed=is_directed + ) + if intersections.size == 0: + continue + elif intersections.shape[0] == 2: + intersections_pts_list_cart.extend(intersections) + else: + intersections_pts_list_cart.append(intersections[0]) + + # Handle unique intersections and check for convex or concave cases + unique_intersections = np.unique(intersections_pts_list_cart, axis=0) + if len(unique_intersections) == 2: + unique_intersection_lonlat = np.array( + [node_xyz_to_lonlat_rad(pt.tolist()) for pt in unique_intersections] + ) + sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0) + pt_lon_min, pt_lon_max = sorted_lonlat[:, 0] + return unique_intersections, pt_lon_min, pt_lon_max + elif len(unique_intersections) != 0: + raise ValueError( + "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly" + ) + elif len(unique_intersections) == 0: + raise ValueError( + "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" + ) + + +def _get_zonal_face_interval( + face_edges_cart, + latitude_cart, + face_latlon_bound, + is_directed=False, + is_latlonface=False, + is_GCA_list=None, +): + """Processes a face polygon represented by edges in Cartesian coordinates + to find intervals where the face intersects with a given latitude. This + function handles directed and undirected Great Circle Arcs (GCAs) and edges + at constant latitude. + + Requires the face edges to be sorted in counter-clockwise order, and the span of the + face in longitude should be less than pi. Also, all arcs/edges length should be within pi. + + Users can specify which edges are GCAs and which are constant latitude using `is_GCA_list`. + However, edges on the equator are always treated as constant latitude edges regardless of + `is_GCA_list`. + + Parameters + ---------- + face_edges_cart : np.ndarray + A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3) + latitude_cart : float + The latitude in cartesian, the normalized Z coordinates. + face_latlon_bound : np.ndarray + The latitude and longitude bounds of the face. Shape: (2, 2), [[lat_min, lat_max], [lon_min, lon_max]] + is_directed : bool, optional + If True, the GCA is considered to be directed (from v0 to v1). If False, the GCA is undirected, + and the smaller circle (less than 180 degrees) side of the GCA is used. Default is False. + is_latlonface : bool, optional, default=False + A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means + That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. + Default is False. This attribute will overwrite the is_latlonface attribute. + is_GCA_list : np.ndarray, optional, default=False + An array indicating if each edge is a GCA (True) or a constant latitude (False). Shape: (n_edges,). + If None, all edges are considered as GCAs. Default is None. + + Returns + ------- + Intervals_df : pd.DataFrame + A DataFrame containing the intervals stored as pandas Intervals for the face. + The columns of the DataFrame are: ['start', 'end'] + """ + face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1] + + # Call the vectorized function to process all edges + ( + unique_intersections, + pt_lon_min, + pt_lon_max, + ) = _get_faces_constLat_intersection_info( + face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed + ) + + # Convert intersection points to longitude-latitude + longitudes = np.array( + [node_xyz_to_lonlat_rad(pt.tolist())[0] for pt in unique_intersections] + ) + + # Handle special wrap-around cases by checking the face bounds + if face_lon_bound_left >= face_lon_bound_right: + if not ( + (pt_lon_max >= np.pi and pt_lon_min >= np.pi) + or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi) + ): + # They're at different sides of the 0-lon, adding wrap-around points + longitudes = np.append(longitudes, [0.0, 2 * np.pi]) + + # Ensure longitudes are sorted + longitudes.sort() + + # Split the sorted longitudes into start and end points of intervals + starts = longitudes[::2] # Start points + ends = longitudes[1::2] # End points + + # Create the intervals DataFrame + Intervals_df = pd.DataFrame({"start": starts, "end": ends}) + # For consistency, sort the intervals by start + interval_df_sorted = Intervals_df.sort_values(by="start").reset_index(drop=True) + return interval_df_sorted + + +def _process_overlapped_intervals(intervals_df): + """Process the overlapped intervals using the sweep line algorithm, + considering multiple intervals per face. + + Parameters + ---------- + intervals_df : pd.DataFrame + The DataFrame that contains the intervals and the corresponding face index. + The columns of the DataFrame are: ['start', 'end', 'face_index'] + + Returns + ------- + dict + A dictionary where keys are face indices and values are their respective contributions to the total length. + float + The total length of all intervals considering their overlaps. + + Example + ------- + >>> intervals_data = [ + ... {'start': 0.0, 'end': 100.0, 'face_index': 0}, + ... {'start': 50.0, 'end': 150.0, 'face_index': 1}, + ... {'start': 140.0, 'end': 150.0, 'face_index': 2} + ... ] + >>> intervals_df = pd.DataFrame(intervals_data) + >>> overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) + >>> print(overlap_contributions) + >>> print(total_length) + """ + events = [] + for idx, row in intervals_df.iterrows(): + events.append((row["start"], "start", row["face_index"])) + events.append((row["end"], "end", row["face_index"])) + + events.sort() # Sort by position and then by start/end + + active_faces = set() + last_position = None + total_length = 0.0 + overlap_contributions = {} + + for position, event_type, face_idx in events: + if last_position is not None and active_faces: + segment_length = position - last_position + segment_weight = segment_length / len(active_faces) if active_faces else 0 + for active_face in active_faces: + overlap_contributions[active_face] = ( + overlap_contributions.get(active_face, 0) + segment_weight + ) + total_length += segment_length + + if event_type == "start": + active_faces.add(face_idx) + elif event_type == "end": + active_faces.remove(face_idx) + + last_position = position + + return overlap_contributions, total_length diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index b0d52cdb7..3d3b03f42 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -116,7 +116,7 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=False): def gca_constLat_intersection( - gca_cart, constLat, fma_disabled=False, verbose=False, is_directed=False + gca_cart, constZ, fma_disabled=False, verbose=False, is_directed=False ): """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a constant latitude line in a Cartesian coordinate system. @@ -128,8 +128,8 @@ def gca_constLat_intersection( Parameters ---------- gca_cart : [2, 3] np.ndarray Cartesian coordinates of the two end points GCA. - constLat : float - The constant latitude of the latitude line. + constZ : float + The constant latitude represented in cartesian of the latitude line. fma_disabled : bool, optional (default=False) If True, the FMA operation is disabled. And a naive `np.cross` is used instead. verbose : bool, optional (default=False) @@ -153,7 +153,6 @@ def gca_constLat_intersection( If running on the Windows system with fma_disabled=False since the C/C++ implementation of FMA in MS Windows is fundamentally broken. (bug report: https://bugs.python.org/msg312480) """ - constZ = np.sin(constLat) x1, x2 = gca_cart if fma_disabled: