Compute Grid Centerpoint using Welzl's algorithm

benchmarks/mpas_ocean.py

@@ -137,6 +137,13 @@ class HoleEdgeIndices(DatasetBenchmark):
     def time_construct_hole_edge_indices(self, resolution):
         ux.grid.geometry._construct_hole_edge_indices(self.uxds.uxgrid.edge_face_connectivity)
 
+class ConstructFaceLatLon(GridBenchmark):
+    def time_welzl(self, resolution):
+        self.uxgrid.construct_face_centers(method='welzl')
+
+    def time_cartesian_averaging(self, resolution):
+        self.uxgrid.construct_face_centers(method='cartesian average')
+
 class CheckNorm:
     param_names = ['resolution']
     params = ['480km', '120km']

docs/internal_api/index.rst

@@ -148,17 +148,26 @@ Coordinates
 .. autosummary::
    :toctree: generated/
 
-   grid.coordinates._lonlat_rad_to_xyz
-   grid.coordinates._xyz_to_lonlat_rad
-   grid.coordinates._xyz_to_lonlat_deg
-   grid.coordinates._normalize_xyz
    grid.coordinates._populate_node_latlon
    grid.coordinates._populate_node_xyz
    grid.coordinates._populate_face_centroids
    grid.coordinates._populate_edge_centroids
    grid.coordinates._construct_face_centroids
    grid.coordinates._construct_edge_centroids
    grid.coordinates._set_desired_longitude_range
+   grid.coordinates._populate_face_centerpoints
+   grid.coordinates._circle_from_two_points
+   grid.coordinates._circle_from_three_points
+   grid.coordinates._is_inside_circle
+   grid.coordinates._welzl_recursive
+   grid.coordinates._smallest_enclosing_circle
+   grid.coordinates._construct_face_centerpoints
+   grid.coordinates._lonlat_rad_to_xyz
+   grid.coordinates._xyz_to_lonlat_rad
+   grid.coordinates._xyz_to_lonlat_rad_no_norm
+   grid.coordinates._xyz_to_lonlat_deg
+   grid.coordinates._normalize_xyz
+   grid.coordinates._normalize_xyz_scalar
 
 
 Arcs
@@ -177,12 +186,12 @@ Utils
 
    grid.utils._newton_raphson_solver_for_gca_constLat
    grid.utils._inv_jacobian
+   grid.utils._angle_of_two_vectors
    grid.utils._swap_first_fill_value_with_last
    grid.utils._get_cartesiain_face_edge_nodes
    grid.utils._get_lonlat_rad_face_edge_nodes
 
 
-
 Validation
 ----------
 .. autosummary::

docs/user_api/index.rst

@@ -238,6 +238,7 @@ Methods
    Grid.get_kd_tree
    Grid.copy
    Grid.isel
+   Grid.construct_face_centers
    Grid.chunk
 
 

test/test_centroids.py

@@ -4,7 +4,8 @@
 import numpy.testing as nt
 import uxarray as ux
 from pathlib import Path
-from uxarray.grid.coordinates import _populate_face_centroids, _populate_edge_centroids, _normalize_xyz
+from uxarray.grid.coordinates import _populate_face_centroids, _populate_edge_centroids, _populate_face_centerpoints, _is_inside_circle, _circle_from_three_points, _circle_from_two_points
+from uxarray.grid.coordinates import _normalize_xyz
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
@@ -63,7 +64,8 @@ def test_centroids_from_mean_verts_scrip(self):
         expected_face_x = uxgrid.face_lon.values
         expected_face_y = uxgrid.face_lat.values
 
-        _populate_face_centroids(uxgrid, repopulate=True)
+        # _populate_face_centroids(uxgrid, repopulate=True)
+        uxgrid.construct_face_centers(method="cartesian average")
 
         # computed_face_x = (uxgrid.face_lon.values + 180) % 360 - 180
         computed_face_x = uxgrid.face_lon.values
@@ -107,3 +109,60 @@ def test_edge_centroids_from_mpas(self):
 
         nt.assert_array_almost_equal(expected_edge_lon, computed_edge_lon)
         nt.assert_array_almost_equal(expected_edge_lat, computed_edge_lat)
+
+class TestCenterPoints(TestCase):
+
+    def test_circle_from_two_points(self):
+        """Test creation of circle from 2 points."""
+        p1 = (0, 0)
+        p2 = (0, 90)
+        center, radius = _circle_from_two_points(p1, p2)
+
+        # The expected radius in radians should be half the angle between the two vectors
+        expected_center = (0.0, 45.0)
+        expected_radius = np.deg2rad(45.0)
+
+        assert np.allclose(center, expected_center), f"Expected center {expected_center}, but got {center}"
+        assert np.allclose(radius, expected_radius), f"Expected radius {expected_radius}, but got {radius}"
+
+    def test_circle_from_three_points(self):
+        """Test creation of circle from 3 points."""
+        p1 = (0, 0)
+        p2 = (0, 90)
+        p3 = (90, 0)
+        center, radius = _circle_from_three_points(p1, p2, p3)
+        expected_radius = np.deg2rad(45.0)
+        expected_center = (30.0, 30.0)
+
+        assert np.allclose(center, expected_center), f"Expected center {expected_center}, but got {center}"
+        assert np.allclose(radius, expected_radius), f"Expected radius {expected_radius}, but got {radius}"
+
+    def test_is_inside_circle(self):
+        """Test if a points is inside the circle."""
+        # Define the circle
+        circle = ((0.0, 0.0), 1)  # Center at lon/lat with a radius in radians (angular measure of the radius)
+
+        # Define test points
+        point_inside = (30.0, 30.0)  # Should be inside the circle
+        point_outside = (90.0, 0.0)  # Should be outside the circle
+
+        # Test _is_inside_circle function
+        assert _is_inside_circle(circle, point_inside), f"Point {point_inside} should be inside the circle."
+        assert not _is_inside_circle(circle, point_outside), f"Point {point_outside} should be outside the circle."
+
+    def test_face_centerpoint(self):
+        """Use points from an actual spherical face and get the centerpoint."""
+
+        points = np.array([(-35.26438968, -45.0), (-36.61769496, -42.0), (-33.78769181, -42.0), (-32.48416571, -45.0)])
+        uxgrid = ux.open_grid(points, latlon=True)
+
+        # Uses the @property from get face_lon/lat - default is average or centroid
+        ctr_lon = uxgrid.face_lon.values[0]
+        ctr_lat = uxgrid.face_lat.values[0]
+
+        # now explicitly get the centerpoints stored to face_lon/lat using welzl's centerpoint algorithm
+        uxgrid.construct_face_centers(method = "welzl")
+
+        # Test the values of the calculated centerpoint, giving high tolerance of two decimal place
+        nt.assert_array_almost_equal(ctr_lon, uxgrid.face_lon.values[0], decimal=2)
+        nt.assert_array_almost_equal(ctr_lat, uxgrid.face_lat.values[0], decimal=2)

uxarray/grid/arcs.py

@@ -6,6 +6,9 @@
     _xyz_to_lonlat_rad_scalar,
     _normalize_xyz_scalar,
 )
+
+from uxarray.grid.utils import _angle_of_2_vectors
+
 from uxarray.constants import ERROR_TOLERANCE, MACHINE_EPSILON
 
 from uxarray.utils.computing import isclose, cross, dot, allclose
@@ -303,31 +306,6 @@ def _decide_pole_latitude(lat1, lat2):
     return closest_pole
 
 
-@njit
-def _angle_of_2_vectors(u, v):
-    """Calculate the angle between two 3D vectors u and v in radians. Can be
-    used to calcualte the span of a GCR.
-
-    Parameters
-    ----------
-    u : numpy.ndarray (float)
-        The first 3D vector.
-    v : numpy.ndarray (float)
-        The second 3D vector.
-
-    Returns
-    -------
-    float
-        The angle between u and v in radians.
-    """
-    v_norm_times_u = np.linalg.norm(v) * u
-    u_norm_times_v = np.linalg.norm(u) * v
-    vec_minus = v_norm_times_u - u_norm_times_v
-    vec_sum = v_norm_times_u + u_norm_times_v
-    angle_u_v_rad = 2 * np.arctan2(np.linalg.norm(vec_minus), np.linalg.norm(vec_sum))
-    return angle_u_v_rad
-
-
 def extreme_gca_latitude(gca_cart, extreme_type):
     """Calculate the maximum or minimum latitude of a great circle arc defined
     by two 3D points.