Compute Grid Centerpoint using Welzl's algorithm

benchmarks/mpas_ocean.py

@@ -178,3 +178,34 @@ def teardown(self, resolution):
 
     def time_construct_hole_edge_indices(self, resolution):
         ux.grid.geometry._construct_hole_edge_indices(self.uxds.uxgrid.edge_face_connectivity)
+
+
+class ConstructFaceLatLon:
+
+    param_names = ['resolution']
+    params = ['480km', '120km']
+
+    def setup(self, resolution):
+        self.uxgrid = ux.open_grid(file_path_dict[resolution][0])
+
+    def teardown(self, resolution):
+        del self.uxgrid
+
+    def time_welzl(self, resolution):
+        from uxarray.grid.coordinates import _construct_face_centerpoints
+
+
+        _construct_face_centerpoints(self.uxgrid.node_lon.values,
+                                     self.uxgrid.node_lat.values,
+                                     self.uxgrid.face_node_connectivity.values,
+                                     self.uxgrid.n_nodes_per_face.values)
+
+
+    def time_cartesian_averaging(self, resolution):
+        from uxarray.grid.coordinates import _construct_face_centroids
+
+        _construct_face_centroids(self.uxgrid.node_x.values,
+                                  self.uxgrid.node_y.values,
+                                  self.uxgrid.node_z.values,
+                                  self.uxgrid.face_node_connectivity.values,
+                                  self.uxgrid.n_nodes_per_face.values)

benchmarks/quad_hexagon.py

@@ -34,3 +34,39 @@ def time_open_dataset(self):
     def peakmem_open_dataset(self):
         """Peak memory usage of a `UxDataset`"""
         uxds = ux.open_dataset(grid_path, data_path)
+
+
+# from uxarray.grid.coordinates import _construct_face_centerpoints, _construct_face_centroids
+
+class ConstructFaceLatLon:
+
+
+    def time_welzl(self):
+        uxgrid = ux.open_grid("/Users/mbook/cp_uxarray/benchmarks/oQU480.grid.nc")
+
+        ux.grid.coordinates._construct_face_centerpoints(uxgrid.node_lon.values,
+                                     uxgrid.node_lat.values,
+                                     uxgrid.face_node_connectivity.values,
+                                     uxgrid.n_nodes_per_face.values)
+
+
+
+
+
+class QuadHexagon2:
+    def time_open_grid(self):
+        """Time to open a `Grid`"""
+        ux.open_grid(grid_path)
+
+    # def mem_open_grid(self):
+    #     """Memory Occupied by a `Grid`"""
+    #     return ux.open_grid(grid_path)
+
+    def peakmem_open_grid2(self):
+        """Peak memory usage of a `Grid`"""
+        uxgrid = ux.open_grid(grid_path)
+
+
+    def time_open_dataset(self):
+        """Time to open a `UxDataset`"""
+        ux.open_dataset(grid_path, data_path)

docs/internal_api/index.rst

@@ -148,17 +148,20 @@ 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
 
 
 Arcs
@@ -180,7 +183,12 @@ Utils
    grid.utils._swap_first_fill_value_with_last
    grid.utils._get_cartesiain_face_edge_nodes
    grid.utils._get_lonlat_rad_face_edge_nodes
-
+   grid.utils._lonlat_rad_to_xyz
+   grid.utils._xyz_to_lonlat_rad
+   grid.utils._xyz_to_lonlat_rad_no_norm
+   grid.utils._xyz_to_lonlat_deg
+   grid.utils._normalize_xyz
+   grid.utils._normalize_xyz_scalar
 
 
 Validation

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.utils 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)

test/test_helpers.py

@@ -12,7 +12,7 @@
 from uxarray.grid.connectivity import _replace_fill_values
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
-from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad
+from uxarray.grid.utils 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

uxarray/grid/arcs.py

@@ -2,7 +2,7 @@
 
 # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place
 
-from uxarray.grid.coordinates import _xyz_to_lonlat_rad_no_norm, _normalize_xyz_scalar
+from uxarray.grid.utils import _xyz_to_lonlat_rad_no_norm, _normalize_xyz_scalar
 from uxarray.constants import ERROR_TOLERANCE
 
 from uxarray.utils.computing import isclose, cross, dot

uxarray/grid/area.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from uxarray.grid.coordinates import _lonlat_rad_to_xyz
+from uxarray.grid.utils import _lonlat_rad_to_xyz
 
 from numba import njit, config
 from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT