DRAFT: Bilinear Remapping

docs/internal_api/index.rst

@@ -142,6 +142,9 @@ Geometry
    grid.geometry._populate_face_latlon_bound
    grid.geometry._populate_bounds
    grid.geometry._construct_hole_edge_indices
+   grid.geometry._point_to_plane
+   grid.geometry._point_to_sphere
+   grid.geometry._ray_casting_plane
 
 Coordinates
 -----------

test/test_geometry.py

@@ -12,8 +12,8 @@
 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, _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, _pole_point_inside_polygon_cartesian, 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
 
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
@@ -30,6 +30,7 @@
 grid_quad_hex = current_path / "meshfiles" / "ugrid" / "quad-hexagon" / "grid.nc"
 grid_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
 grid_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
+grid_mpas_2 = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
 
 # List of grid files to test
 grid_files_latlonBound = [grid_quad_hex, grid_geoflow, gridfile_CSne8, grid_mpas]
@@ -81,7 +82,7 @@ def test_pole_point_inside_polygon_from_vertice_north(self):
                                    [vertices[1], vertices[2]],
                                    [vertices[2], vertices[3]],
                                    [vertices[3], vertices[0]]])
-
+        print(face_edge_cart)
         # Check if the North pole is inside the polygon
         result = _pole_point_inside_polygon_cartesian(
             'North', face_edge_cart)

test/test_remap.py

@@ -22,6 +22,7 @@
 dsfile_v2_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "v2.nc"
 dsfile_v3_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "v3.nc"
 mpasfile_QU = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
+mpasfile_QU_2 = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
 
 
 class TestNearestNeighborRemap(TestCase):
@@ -224,7 +225,6 @@ def test_preserve_coordinates(self):
         assert "time" in res.coords
 
 
-
 class TestInverseDistanceWeightedRemapping(TestCase):
     """Testing for inverse distance weighted remapping."""
 
@@ -429,3 +429,13 @@ def test_value_errors(self):
                 destination_grid=destination_grid,
                 remap_to="nodes", power=6
             )
+
+
+class TestBilinearRemapping(TestCase):
+
+    def test_uxda_remap(self):
+        source_uxds = ux.open_dataset(mpasfile_QU, mpasfile_QU)
+        destination = ux.open_dataset(mpasfile_QU_2, mpasfile_QU_2)
+
+        bilinear_remap = source_uxds['latVertex'].remap.bilinear(destination.uxgrid, remap_to='face centers',
+                                                                                     coord_type='spherical')

uxarray/grid/geometry.py

@@ -20,6 +20,7 @@
 )
 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,
@@ -1236,11 +1237,13 @@ def _populate_face_latlon_bound(
                 face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_min, node1_lon_rad])
                 )
+
                 face_latlon_array[0, 1] = math.pi / 2  # Upper latitude bound
             else:
                 face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_max, node1_lon_rad])
                 )
+
                 face_latlon_array[0, 0] = -math.pi / 2  # Lower latitude bound
 
         # Adjust longitude bounds globally if the pole is centrally inside the polygon