From 60a2cbbcb2aac5f47b31fc28e6c183d675e5b36b Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Tue, 11 Jun 2024 06:08:18 -0500 Subject: [PATCH 01/29] o #803 initial implementation of Welzl's algorithm to calculate grid centerpoint. Need to use great circle distance and add/fix tests and data types in the algo --- test/test_centroids.py | 10 ++- uxarray/grid/coordinates.py | 145 ++++++++++++++++++++++++++++++++++++ uxarray/grid/grid.py | 54 ++++++++++++++ 3 files changed, 208 insertions(+), 1 deletion(-) diff --git a/test/test_centroids.py b/test/test_centroids.py index 6eb69542f..2e2563d62 100644 --- a/test/test_centroids.py +++ b/test/test_centroids.py @@ -4,7 +4,7 @@ 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, _normalize_xyz, _populate_face_centerpoints current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -107,3 +107,11 @@ 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) + + 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) + _populate_face_centerpoints(uxgrid) + print(uxgrid.face_lon_ctrpt) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index ba0b27e57..3f59d65e3 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -3,6 +3,10 @@ import warnings +import random +from typing import List, Tuple +from shapely.geometry import Point, LineString, Polygon + from uxarray.constants import ERROR_TOLERANCE from uxarray.conventions import ugrid @@ -10,6 +14,10 @@ from numba import njit +# Type aliases for better readability +PointT = Tuple[float, float] +Circle = Tuple[Point, float] + @njit(cache=True) def _lonlat_rad_to_xyz( @@ -227,6 +235,56 @@ def _populate_face_centroids(grid, repopulate=False): ) +def _populate_face_centerpoints(grid, repopulate=False): + """Finds the centerpoints of faces using Welzl's algorithm based. + + Parameters + ---------- + repopulate : bool, optional + Bool used to turn on/off repopulating the face coordinates of the centerpoints + """ + warnings.warn("This cannot be guaranteed to work correctly on concave polygons") + + node_lon = grid.node_lon.values + node_lat = grid.node_lat.values + + centerpoint_lat = [] + centerpoint_lon = [] + + face_nodes = grid.face_node_connectivity.values + n_nodes_per_face = grid.n_nodes_per_face.values + + if "face_lon_ctrpt" not in grid._ds or repopulate: + # Construct the centerpoints if there are none stored + if "face_x_ctrpt" not in grid._ds: + centerpoint_lon, centerpoint_lat = _construct_face_centerpoints( + node_lon, node_lat, face_nodes, n_nodes_per_face + ) + ctrpt_x, ctrpt_y, ctrpt_z = _lonlat_rad_to_xyz(centerpoint_lon, centerpoint_lat) + + # Populate the centerpoints + if "face_lon_ctrpt" not in grid._ds or repopulate: + grid._ds["face_lon_ctrpt"] = xr.DataArray( + centerpoint_lon, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LON_ATTRS + ) + grid._ds["face_lat_ctrpt"] = xr.DataArray( + centerpoint_lat, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LAT_ATTRS + ) + + if "face_x_ctrpt" not in grid._ds or repopulate: + grid._ds["face_x_ctrpt"] = xr.DataArray( + ctrpt_x, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_X_ATTRS + ) + + grid._ds["face_y_ctrpt"] = xr.DataArray( + ctrpt_y, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Y_ATTRS + ) + + grid._ds["face_z_ctrpt"] = xr.DataArray( + ctrpt_z, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Z_ATTRS + ) + + def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face): """Constructs the xyz centroid coordinate for each face using Cartesian Averaging.""" @@ -244,6 +302,93 @@ def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_fa return _normalize_xyz(centroid_x, centroid_y, centroid_z) +def distance(a: PointT, b: PointT) -> float: + return Point(a).distance(Point(b)) + + +def is_in_circle(c: Circle, p: PointT) -> bool: + center, radius = c + return Point(center).distance(Point(p)) <= radius + + +def circle_from_two_points(p1: PointT, p2: PointT) -> Circle: + print("p1: ", p1, "p2: ", p2) + line = LineString([Point(p1), Point(p2)]) + center = line.centroid + radius = center.distance(Point(p1)) + return (center.x, center.y), radius + + +def circle_from_three_points(p1: PointT, p2: PointT, p3: PointT) -> Circle: + triangle = Polygon([p1, p2, p3]) + center = triangle.centroid + radius = ( + center.distance(Point(p1)) + + center.distance(Point(p2)) + + center.distance(Point(p3)) + ) / 3 + return (center.x, center.y), radius + + +def welzl(points: List[PointT], boundary: List[PointT] = []) -> Circle: + if not points or len(boundary) == 3: + if len(boundary) == 0: + return ((0, 0), 0) + elif len(boundary) == 1: + return boundary[0], 0 + elif len(boundary) == 2: + return circle_from_two_points(boundary[0], boundary[1]) + elif len(boundary) == 3: + return circle_from_three_points(boundary[0], boundary[1], boundary[2]) + + p = points.pop() + circle = welzl(points, boundary) + + if is_in_circle(circle, p): + points.append(p) + return circle + else: + boundary.append(p) + result = welzl(points, boundary) + boundary.pop() + points.append(p) + return result + + +def smallest_enclosing_circle(points: List[PointT]) -> Circle: + shuffled_points = points[:] + random.shuffle(shuffled_points) + return welzl(shuffled_points) + + +def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_face): + """Constructs the face centerpoint using Welzl's algorithm3.""" + + points = np.column_stack((node_lon, node_lat)).tolist() + + ctrpt_lon = np.zeros((face_nodes.shape[0]), dtype=np.float64) + ctrpt_lat = np.zeros((face_nodes.shape[0]), dtype=np.float64) + + for face_idx, n_max_nodes in enumerate(n_nodes_per_face): + points_array = np.column_stack( + ( + node_lon[face_nodes[face_idx, 0:n_max_nodes]], + node_lat[face_nodes[face_idx, 0:n_max_nodes]], + ) + ) + points = [tuple(point) for point in points_array.tolist()] + circle = smallest_enclosing_circle(points) + ctrpt_lon[face_idx] = circle[0][0] + ctrpt_lat[face_idx] = circle[0][1] + + return ctrpt_lon, ctrpt_lat + # centerpoint_x, centerpoint_y, centerpoint_z = _lonlat_rad_to_xyz(float(circle[0][0]), float(circle[0][1])) + # print(points) + # print("centerpoint_x: ", centerpoint_x, "centerpoint_y: ", centerpoint_y, "centerpoint_z: ", centerpoint_z) + + # return _normalize_xyz(centerpoint_x, centerpoint_y, centerpoint_z) + + def _populate_edge_centroids(grid, repopulate=False): """Finds the centroids using cartesian averaging of the edges based off the vertices. The centroid is defined as the average of the x, y, z diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 716a1983b..105846586 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -707,6 +707,60 @@ def face_z(self) -> xr.DataArray: _populate_face_centroids(self) return self._ds["face_z"] + @property + def face_lon_ctrpt(self) -> xr.DataArray: + """Longitude of the center of each face in degrees. + + Dimensions: ``(n_face, )`` + """ + if "face_lon_ctrpt" not in self._ds: + _populate_face_centroids(self) + _set_desired_longitude_range(self._ds) + return self._ds["face_lon_ctrpt"] + + @property + def face_lat_ctrpt(self) -> xr.DataArray: + """Latitude of the center of each face in degrees. + + Dimensions: ``(n_face, )`` + """ + if "face_lat_ctrpt" not in self._ds: + _populate_face_centroids(self) + _set_desired_longitude_range(self._ds) + + return self._ds["face_lat_ctrpt"] + + @property + def face_x_ctrpt(self) -> xr.DataArray: + """Cartesian x location of the center of each face in meters. + + Dimensions: ``(n_face, )`` + """ + if "face_x_ctrpt" not in self._ds: + _populate_face_centroids(self) + + return self._ds["face_x_ctrpt"] + + @property + def face_y_ctrpt(self) -> xr.DataArray: + """Cartesian y location of the center of each face in meters. + + Dimensions: ``(n_face, )`` + """ + if "face_y_ctrpt" not in self._ds: + _populate_face_centroids(self) + return self._ds["face_y_ctrpt"] + + @property + def face_z_ctrpt(self) -> xr.DataArray: + """Cartesian z location of the center of each face in meters. + + Dimensions: ``(n_face, )`` + """ + if "face_z_ctrpt" not in self._ds: + _populate_face_centroids(self) + return self._ds["face_z_ctrpt"] + @property def face_node_connectivity(self) -> xr.DataArray: """Indices of the nodes that make up each face. From 3c41af31405edeb3145f2a524121635044edf9c0 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Tue, 18 Jun 2024 14:14:58 -0500 Subject: [PATCH 02/29] o Update grid class and add asserts for test --- test/test_centroids.py | 10 +++++++++- uxarray/grid/coordinates.py | 7 ------- uxarray/grid/grid.py | 11 ++++++----- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/test/test_centroids.py b/test/test_centroids.py index 2e2563d62..0bb7cc23d 100644 --- a/test/test_centroids.py +++ b/test/test_centroids.py @@ -114,4 +114,12 @@ def test_face_centerpoint(self): 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) _populate_face_centerpoints(uxgrid) - print(uxgrid.face_lon_ctrpt) + + # the expected centerpoint + ctr_lon = -34.55093034 + ctr_lat = -43.5 + + # Test the values of the calculated centerpoint + nt.assert_array_almost_equal(ctr_lon, uxgrid.face_lon_ctrpt.values[0], decimal=3) + nt.assert_array_almost_equal(ctr_lat, uxgrid.face_lat_ctrpt.values[0], decimal=3) + \ No newline at end of file diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 3f59d65e3..034dcc55f 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -312,7 +312,6 @@ def is_in_circle(c: Circle, p: PointT) -> bool: def circle_from_two_points(p1: PointT, p2: PointT) -> Circle: - print("p1: ", p1, "p2: ", p2) line = LineString([Point(p1), Point(p2)]) center = line.centroid radius = center.distance(Point(p1)) @@ -382,12 +381,6 @@ def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_fac ctrpt_lat[face_idx] = circle[0][1] return ctrpt_lon, ctrpt_lat - # centerpoint_x, centerpoint_y, centerpoint_z = _lonlat_rad_to_xyz(float(circle[0][0]), float(circle[0][1])) - # print(points) - # print("centerpoint_x: ", centerpoint_x, "centerpoint_y: ", centerpoint_y, "centerpoint_z: ", centerpoint_z) - - # return _normalize_xyz(centerpoint_x, centerpoint_y, centerpoint_z) - def _populate_edge_centroids(grid, repopulate=False): """Finds the centroids using cartesian averaging of the edges based off the diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 105846586..c7ccb55e2 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -27,6 +27,7 @@ from uxarray.grid.coordinates import ( _populate_face_centroids, _populate_edge_centroids, + _populate_face_centerpoints, _set_desired_longitude_range, _populate_node_latlon, _populate_node_xyz, @@ -714,7 +715,7 @@ def face_lon_ctrpt(self) -> xr.DataArray: Dimensions: ``(n_face, )`` """ if "face_lon_ctrpt" not in self._ds: - _populate_face_centroids(self) + _populate_face_centerpoints(self) _set_desired_longitude_range(self._ds) return self._ds["face_lon_ctrpt"] @@ -725,7 +726,7 @@ def face_lat_ctrpt(self) -> xr.DataArray: Dimensions: ``(n_face, )`` """ if "face_lat_ctrpt" not in self._ds: - _populate_face_centroids(self) + _populate_face_centerpoints(self) _set_desired_longitude_range(self._ds) return self._ds["face_lat_ctrpt"] @@ -737,7 +738,7 @@ def face_x_ctrpt(self) -> xr.DataArray: Dimensions: ``(n_face, )`` """ if "face_x_ctrpt" not in self._ds: - _populate_face_centroids(self) + _populate_face_centerpoints(self) return self._ds["face_x_ctrpt"] @@ -748,7 +749,7 @@ def face_y_ctrpt(self) -> xr.DataArray: Dimensions: ``(n_face, )`` """ if "face_y_ctrpt" not in self._ds: - _populate_face_centroids(self) + _populate_face_centerpoints(self) return self._ds["face_y_ctrpt"] @property @@ -758,7 +759,7 @@ def face_z_ctrpt(self) -> xr.DataArray: Dimensions: ``(n_face, )`` """ if "face_z_ctrpt" not in self._ds: - _populate_face_centroids(self) + _populate_face_centerpoints(self) return self._ds["face_z_ctrpt"] @property From 33c445fc8a85cc762223a6d60c6ee6d1491c330b Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Wed, 19 Jun 2024 12:46:55 -0500 Subject: [PATCH 03/29] o typo fix --- test/test_centroids.py | 1 - uxarray/grid/coordinates.py | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_centroids.py b/test/test_centroids.py index 0bb7cc23d..628276b70 100644 --- a/test/test_centroids.py +++ b/test/test_centroids.py @@ -122,4 +122,3 @@ def test_face_centerpoint(self): # Test the values of the calculated centerpoint nt.assert_array_almost_equal(ctr_lon, uxgrid.face_lon_ctrpt.values[0], decimal=3) nt.assert_array_almost_equal(ctr_lat, uxgrid.face_lat_ctrpt.values[0], decimal=3) - \ No newline at end of file diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 034dcc55f..278f32d58 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -382,6 +382,7 @@ def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_fac return ctrpt_lon, ctrpt_lat + def _populate_edge_centroids(grid, repopulate=False): """Finds the centroids using cartesian averaging of the edges based off the vertices. The centroid is defined as the average of the x, y, z From 9b3f06643fb328dcb3a9de8445e9e1964699cb91 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Wed, 26 Jun 2024 15:07:56 -0500 Subject: [PATCH 04/29] o overhaul to not use tuples and go with numpy array, document and fix test case asserts --- test/test_centroids.py | 12 +-- uxarray/grid/coordinates.py | 186 +++++++++++++++++++++++++----------- 2 files changed, 138 insertions(+), 60 deletions(-) diff --git a/test/test_centroids.py b/test/test_centroids.py index 628276b70..f609822fd 100644 --- a/test/test_centroids.py +++ b/test/test_centroids.py @@ -115,10 +115,10 @@ def test_face_centerpoint(self): uxgrid = ux.open_grid(points, latlon=True) _populate_face_centerpoints(uxgrid) - # the expected centerpoint - ctr_lon = -34.55093034 - ctr_lat = -43.5 + # the expected centerpoint should be close to centroid for this case + ctr_lon = uxgrid.face_lon.values[0] + ctr_lat = uxgrid.face_lat.values[0] - # Test the values of the calculated centerpoint - nt.assert_array_almost_equal(ctr_lon, uxgrid.face_lon_ctrpt.values[0], decimal=3) - nt.assert_array_almost_equal(ctr_lat, uxgrid.face_lat_ctrpt.values[0], decimal=3) + # Test the values of the calculated centerpoint, giving high tolerance of one decimal place + nt.assert_array_almost_equal(ctr_lon, uxgrid.face_lon_ctrpt.values[0], decimal=1) + nt.assert_array_almost_equal(ctr_lat, uxgrid.face_lat_ctrpt.values[0], decimal=1) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 278f32d58..252451365 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -3,9 +3,6 @@ import warnings -import random -from typing import List, Tuple -from shapely.geometry import Point, LineString, Polygon from uxarray.constants import ERROR_TOLERANCE from uxarray.conventions import ugrid @@ -14,10 +11,6 @@ from numba import njit -# Type aliases for better readability -PointT = Tuple[float, float] -Circle = Tuple[Point, float] - @njit(cache=True) def _lonlat_rad_to_xyz( @@ -236,14 +229,18 @@ def _populate_face_centroids(grid, repopulate=False): def _populate_face_centerpoints(grid, repopulate=False): - """Finds the centerpoints of faces using Welzl's algorithm based. + """Calculates the face centerpoints using Welzl's algorithm. It is a + randomized algorithm for finding the center and radius of the smallest + circle that encloses a set of points. It is here adapted to work on a unit + sphere. Also, this algorithm cannot be guaranteed to work on concave + polygons. Parameters ---------- repopulate : bool, optional - Bool used to turn on/off repopulating the face coordinates of the centerpoints + Bool used to turn on/off repopulating the face coordinates of the centerpoints, default is False """ - warnings.warn("This cannot be guaranteed to work correctly on concave polygons") + # warnings.warn("This cannot be guaranteed to work correctly on concave polygons") node_lon = grid.node_lon.values node_lat = grid.node_lat.values @@ -254,15 +251,17 @@ def _populate_face_centerpoints(grid, repopulate=False): face_nodes = grid.face_node_connectivity.values n_nodes_per_face = grid.n_nodes_per_face.values + # Check if the centerpoints are already populated if "face_lon_ctrpt" not in grid._ds or repopulate: # Construct the centerpoints if there are none stored if "face_x_ctrpt" not in grid._ds: centerpoint_lon, centerpoint_lat = _construct_face_centerpoints( node_lon, node_lat, face_nodes, n_nodes_per_face ) + # get the cartesian coordinates of the centerpoints ctrpt_x, ctrpt_y, ctrpt_z = _lonlat_rad_to_xyz(centerpoint_lon, centerpoint_lat) - # Populate the centerpoints + # set the grid variables for centerpoints if "face_lon_ctrpt" not in grid._ds or repopulate: grid._ds["face_lon_ctrpt"] = xr.DataArray( centerpoint_lon, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LON_ATTRS @@ -302,81 +301,160 @@ def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_fa return _normalize_xyz(centroid_x, centroid_y, centroid_z) -def distance(a: PointT, b: PointT) -> float: - return Point(a).distance(Point(b)) +def haversine_distance(point1, point2): + """Calculate the great-circle distance between two points on a unit sphere + using the Haversine formula. + + Parameters: + - point1: A tuple containing the latitude and longitude of the first point. + - point2: A tuple containing the latitude and longitude of the second point. + + Returns: + - The distance between the two points on the unit sphere. + """ + R = 1.0 # Radius of the Earth assumed to be 1 (unit sphere) + lat1, lon1 = np.radians(point1) + lat2, lon2 = np.radians(point2) + + dlat = lat2 - lat1 + dlon = lon2 - lon1 + + a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2 + c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)) + + distance = R * c + return distance -def is_in_circle(c: Circle, p: PointT) -> bool: - center, radius = c - return Point(center).distance(Point(p)) <= radius +def circle_from_two_points(p1, p2): + """Calculate the smallest circle that encloses two points on a unit sphere. + Parameters: + - p1: The first point as a tuple of (latitude, longitude). + - p2: The second point as a tuple of (latitude, longitude). + + Returns: + - A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + """ + center_lat = (p1[0] + p2[0]) / 2 + center_lon = (p1[1] + p2[1]) / 2 + center = (center_lat, center_lon) + radius = haversine_distance(p1, p2) / 2 + return center, radius -def circle_from_two_points(p1: PointT, p2: PointT) -> Circle: - line = LineString([Point(p1), Point(p2)]) - center = line.centroid - radius = center.distance(Point(p1)) - return (center.x, center.y), radius +def circle_from_three_points(p1, p2, p3): + """Calculate the smallest circle that encloses three points on a unit + sphere. This is a placeholder implementation. -def circle_from_three_points(p1: PointT, p2: PointT, p3: PointT) -> Circle: - triangle = Polygon([p1, p2, p3]) - center = triangle.centroid + Parameters: + - p1, p2, p3: Three points as tuples of (latitude, longitude). + + Returns: + - A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + """ + center = p1 # Placeholder center radius = ( - center.distance(Point(p1)) - + center.distance(Point(p2)) - + center.distance(Point(p3)) - ) / 3 - return (center.x, center.y), radius + max( + haversine_distance(p1, p2), + haversine_distance(p1, p3), + haversine_distance(p2, p3), + ) + / 2 + ) + return center, radius + + +def is_inside_circle(circle, point): + """Check if a point is inside a given circle on a unit sphere. + + Parameters: + - circle: A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + - point: The point to check, as a tuple of (latitude, longitude). + + Returns: + - True if the point is inside the circle, False otherwise. + """ + center, radius = circle + return haversine_distance(center, point) <= radius + +def welzl_recursive(points, boundary, R): + """Recursive helper function for Welzl's algorithm to find the smallest + enclosing circle. -def welzl(points: List[PointT], boundary: List[PointT] = []) -> Circle: - if not points or len(boundary) == 3: + Parameters: + - points: The set of points to consider. + - boundary: The current boundary points of the minimal enclosing circle. + - R: The current minimal enclosing circle. + + Returns: + - The smallest enclosing circle as a tuple of center and radius. + """ + # Base case: no points or boundary has 3 points + if len(points) == 0 or len(boundary) == 3: + # Construct the minimal circle based on the number of boundary points if len(boundary) == 0: - return ((0, 0), 0) + return R elif len(boundary) == 1: - return boundary[0], 0 + return (boundary[0], 0) elif len(boundary) == 2: return circle_from_two_points(boundary[0], boundary[1]) elif len(boundary) == 3: return circle_from_three_points(boundary[0], boundary[1], boundary[2]) - p = points.pop() - circle = welzl(points, boundary) + # Choose a point from the set and remove it + p = points[-1] + temp_points = np.delete(points, -1, axis=0) + circle = welzl_recursive(temp_points, boundary, R) - if is_in_circle(circle, p): - points.append(p) + # Check if the chosen point is inside the current circle + if circle and is_inside_circle(circle, p): return circle else: - boundary.append(p) - result = welzl(points, boundary) - boundary.pop() - points.append(p) - return result + # If not, the point must be on the boundary of the minimal enclosing circle + return welzl_recursive(temp_points, np.append(boundary, [p], axis=0), R) + +def smallest_enclosing_circle(points): + """Find the smallest circle that encloses all given points on a unit sphere + using Welzl's algorithm. -def smallest_enclosing_circle(points: List[PointT]) -> Circle: - shuffled_points = points[:] - random.shuffle(shuffled_points) - return welzl(shuffled_points) + Parameters: + - points: An array of points as tuples of (latitude, longitude). + + Returns: + - The smallest enclosing circle as a tuple of center and radius. + """ + np.random.shuffle( + points + ) # Randomize the input to increase the chance of an optimal solution + return welzl_recursive(points, np.empty((0, 2)), None) def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_face): - """Constructs the face centerpoint using Welzl's algorithm3.""" + """Constructs the face centerpoint using Welzl's algorithm. - points = np.column_stack((node_lon, node_lat)).tolist() + Parameters: + - node_lon: Longitudes of the nodes. + - node_lat: Latitudes of the nodes. + - face_nodes: Indices of nodes per face. + - n_nodes_per_face: Number of nodes per face. - ctrpt_lon = np.zeros((face_nodes.shape[0]), dtype=np.float64) - ctrpt_lat = np.zeros((face_nodes.shape[0]), dtype=np.float64) + Returns: + - Two arrays containing the longitudes and latitudes of the centerpoints. + """ + ctrpt_lon = np.zeros(face_nodes.shape[0], dtype=np.float64) + ctrpt_lat = np.zeros(face_nodes.shape[0], dtype=np.float64) for face_idx, n_max_nodes in enumerate(n_nodes_per_face): points_array = np.column_stack( ( - node_lon[face_nodes[face_idx, 0:n_max_nodes]], - node_lat[face_nodes[face_idx, 0:n_max_nodes]], + node_lon[face_nodes[face_idx, :n_max_nodes]], + node_lat[face_nodes[face_idx, :n_max_nodes]], ) ) - points = [tuple(point) for point in points_array.tolist()] - circle = smallest_enclosing_circle(points) + circle = smallest_enclosing_circle(points_array) ctrpt_lon[face_idx] = circle[0][0] ctrpt_lat[face_idx] = circle[0][1] From 504bc08617c97d7a58cac79f93e9b2489fb8fb63 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Wed, 26 Jun 2024 18:09:03 -0500 Subject: [PATCH 05/29] o Conform to formatting standards --- uxarray/grid/coordinates.py | 141 ++++++++++++++++++++++++++---------- 1 file changed, 103 insertions(+), 38 deletions(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 252451365..5cae660e3 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -237,8 +237,10 @@ def _populate_face_centerpoints(grid, repopulate=False): Parameters ---------- + grid : Grid + The grid containing the nodes and faces. repopulate : bool, optional - Bool used to turn on/off repopulating the face coordinates of the centerpoints, default is False + Bool used to turn on/off repopulating the face coordinates of the centerpoints, default is False. """ # warnings.warn("This cannot be guaranteed to work correctly on concave polygons") @@ -286,7 +288,26 @@ def _populate_face_centerpoints(grid, repopulate=False): def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face): """Constructs the xyz centroid coordinate for each face using Cartesian - Averaging.""" + Averaging. + + Parameters + ---------- + node_x : numpy.ndarray + X coordinates of the nodes. + node_y : numpy.ndarray + Y coordinates of the nodes. + node_z : numpy.ndarray + Z coordinates of the nodes. + face_nodes : numpy.ndarray + Indices of nodes per face. + n_nodes_per_face : numpy.ndarray + Number of nodes per face. + + Returns + ------- + tuple + The x, y, and z coordinates of the centroids. + """ centroid_x = np.zeros((face_nodes.shape[0]), dtype=np.float64) centroid_y = np.zeros((face_nodes.shape[0]), dtype=np.float64) @@ -305,12 +326,17 @@ def haversine_distance(point1, point2): """Calculate the great-circle distance between two points on a unit sphere using the Haversine formula. - Parameters: - - point1: A tuple containing the latitude and longitude of the first point. - - point2: A tuple containing the latitude and longitude of the second point. + Parameters + ---------- + point1 : tuple + A tuple containing the latitude and longitude of the first point. + point2 : tuple + A tuple containing the latitude and longitude of the second point. - Returns: - - The distance between the two points on the unit sphere. + Returns + ------- + float + The distance between the two points on the unit sphere. """ R = 1.0 # Radius of the Earth assumed to be 1 (unit sphere) lat1, lon1 = np.radians(point1) @@ -329,12 +355,17 @@ def haversine_distance(point1, point2): def circle_from_two_points(p1, p2): """Calculate the smallest circle that encloses two points on a unit sphere. - Parameters: - - p1: The first point as a tuple of (latitude, longitude). - - p2: The second point as a tuple of (latitude, longitude). + Parameters + ---------- + p1 : tuple + The first point as a tuple of (latitude, longitude). + p2 : tuple + The second point as a tuple of (latitude, longitude). - Returns: - - A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + Returns + ------- + tuple + A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. """ center_lat = (p1[0] + p2[0]) / 2 center_lon = (p1[1] + p2[1]) / 2 @@ -347,11 +378,19 @@ def circle_from_three_points(p1, p2, p3): """Calculate the smallest circle that encloses three points on a unit sphere. This is a placeholder implementation. - Parameters: - - p1, p2, p3: Three points as tuples of (latitude, longitude). + Parameters + ---------- + p1 : tuple + The first point. + p2 : tuple + The second point. + p3 : tuple + The third point. - Returns: - - A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + Returns + ------- + tuple + A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. """ center = p1 # Placeholder center radius = ( @@ -368,12 +407,17 @@ def circle_from_three_points(p1, p2, p3): def is_inside_circle(circle, point): """Check if a point is inside a given circle on a unit sphere. - Parameters: - - circle: A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. - - point: The point to check, as a tuple of (latitude, longitude). + Parameters + ---------- + circle : tuple + A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + point : tuple + The point to check, as a tuple of (latitude, longitude). - Returns: - - True if the point is inside the circle, False otherwise. + Returns + ------- + bool + True if the point is inside the circle, False otherwise. """ center, radius = circle return haversine_distance(center, point) <= radius @@ -383,13 +427,19 @@ def welzl_recursive(points, boundary, R): """Recursive helper function for Welzl's algorithm to find the smallest enclosing circle. - Parameters: - - points: The set of points to consider. - - boundary: The current boundary points of the minimal enclosing circle. - - R: The current minimal enclosing circle. + Parameters + ---------- + points : numpy.ndarray + The set of points to consider. + boundary : numpy.ndarray + The current boundary points of the minimal enclosing circle. + R : tuple + The current minimal enclosing circle. - Returns: - - The smallest enclosing circle as a tuple of center and radius. + Returns + ------- + tuple + The smallest enclosing circle as a tuple of center and radius. """ # Base case: no points or boundary has 3 points if len(points) == 0 or len(boundary) == 3: @@ -420,11 +470,15 @@ def smallest_enclosing_circle(points): """Find the smallest circle that encloses all given points on a unit sphere using Welzl's algorithm. - Parameters: - - points: An array of points as tuples of (latitude, longitude). + Parameters + ---------- + points : numpy.ndarray + An array of points as tuples of (latitude, longitude). - Returns: - - The smallest enclosing circle as a tuple of center and radius. + Returns + ------- + tuple + The smallest enclosing circle as a tuple of center and radius. """ np.random.shuffle( points @@ -435,14 +489,25 @@ def smallest_enclosing_circle(points): def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_face): """Constructs the face centerpoint using Welzl's algorithm. - Parameters: - - node_lon: Longitudes of the nodes. - - node_lat: Latitudes of the nodes. - - face_nodes: Indices of nodes per face. - - n_nodes_per_face: Number of nodes per face. + Parameters + ---------- + node_lon : array_like + Longitudes of the nodes. + node_lat : array_like + Latitudes of the nodes. + face_nodes : array_like + Indices of nodes per face. + n_nodes_per_face : int + Number of nodes per face. + + Returns + ------- + tuple of numpy.ndarray + Two arrays containing the longitudes and latitudes of the centerpoints. - Returns: - - Two arrays containing the longitudes and latitudes of the centerpoints. + Notes + ----- + This function calculates the centerpoints of faces defined by nodes on a sphere, using Welzl's algorithm to find the smallest enclosing circle for each face. """ ctrpt_lon = np.zeros(face_nodes.shape[0], dtype=np.float64) ctrpt_lat = np.zeros(face_nodes.shape[0], dtype=np.float64) From c62d2b7de0375798ce952a3f4bb62e22e57b8c29 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Fri, 28 Jun 2024 22:39:47 -0500 Subject: [PATCH 06/29] o fix bugs that reversed the coordinate ordering --- uxarray/grid/coordinates.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 5cae660e3..485d57593 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -17,8 +17,8 @@ def _lonlat_rad_to_xyz( lon: Union[np.ndarray, float], lat: Union[np.ndarray, float], ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Converts Spherical latitude and longitude coordinates into Cartesian x, - y, z coordinates.""" + """Converts Spherical lon and lat coordinates into Cartesian x, y, z + coordinates.""" x = np.cos(lon) * np.cos(lat) y = np.sin(lon) * np.cos(lat) z = np.sin(lat) @@ -128,7 +128,7 @@ def _normalize_xyz( def _populate_node_latlon(grid) -> None: - """Populates the latitude and longitude coordinates of a Grid (`node_lon`, + """Populates the lon and lat coordinates of a Grid (`node_lon`, `node_lat`)""" lon_rad, lat_rad = _xyz_to_lonlat_rad( grid.node_x.values, grid.node_y.values, grid.node_z.values @@ -329,9 +329,9 @@ def haversine_distance(point1, point2): Parameters ---------- point1 : tuple - A tuple containing the latitude and longitude of the first point. + A tuple containing the lon and lat of the first point. point2 : tuple - A tuple containing the latitude and longitude of the second point. + A tuple containing the lon and lat of the second point. Returns ------- @@ -339,8 +339,8 @@ def haversine_distance(point1, point2): The distance between the two points on the unit sphere. """ R = 1.0 # Radius of the Earth assumed to be 1 (unit sphere) - lat1, lon1 = np.radians(point1) - lat2, lon2 = np.radians(point2) + lon1, lat1 = np.radians(point1) + lon2, lat2 = np.radians(point2) dlat = lat2 - lat1 dlon = lon2 - lon1 @@ -365,11 +365,11 @@ def circle_from_two_points(p1, p2): Returns ------- tuple - A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + A tuple containing the center (as a tuple of lon and lat) and the radius of the circle. """ - center_lat = (p1[0] + p2[0]) / 2 - center_lon = (p1[1] + p2[1]) / 2 - center = (center_lat, center_lon) + center_lon = (p1[0] + p2[0]) / 2 + center_lat = (p1[1] + p2[1]) / 2 + center = (center_lon, center_lat) radius = haversine_distance(p1, p2) / 2 return center, radius @@ -390,7 +390,7 @@ def circle_from_three_points(p1, p2, p3): Returns ------- tuple - A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + A tuple containing the center (as a tuple of lon and lat) and the radius of the circle. """ center = p1 # Placeholder center radius = ( @@ -410,9 +410,9 @@ def is_inside_circle(circle, point): Parameters ---------- circle : tuple - A tuple containing the center (as a tuple of latitude and longitude) and the radius of the circle. + A tuple containing the center (as a tuple of lon and lat) and the radius of the circle. point : tuple - The point to check, as a tuple of (latitude, longitude). + The point to check, as a tuple of (lon, lat). Returns ------- @@ -473,7 +473,7 @@ def smallest_enclosing_circle(points): Parameters ---------- points : numpy.ndarray - An array of points as tuples of (latitude, longitude). + An array of points as tuples of (lon, lat). Returns ------- @@ -507,7 +507,8 @@ def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_fac Notes ----- - This function calculates the centerpoints of faces defined by nodes on a sphere, using Welzl's algorithm to find the smallest enclosing circle for each face. + This function calculates the centerpoints of faces defined by nodes on a sphere, using Welzl's algorithm to + find the smallest enclosing circle for each face. """ ctrpt_lon = np.zeros(face_nodes.shape[0], dtype=np.float64) ctrpt_lat = np.zeros(face_nodes.shape[0], dtype=np.float64) From d09499c977289592fef0345b6fe311bf467374c2 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Mon, 8 Jul 2024 19:03:28 -0500 Subject: [PATCH 07/29] o Move some fns from coordinates.py to utils as they caused circular dependency (use with arcs and arcs use coordinates). o Remove new routine in favor of using the existing angle b/w vectors to calculate distance. --- test/test_centroids.py | 3 +- test/test_helpers.py | 2 +- uxarray/grid/arcs.py | 2 +- uxarray/grid/area.py | 2 +- uxarray/grid/coordinates.py | 172 ++++-------------------------------- uxarray/grid/utils.py | 119 +++++++++++++++++++++++++ uxarray/io/_exodus.py | 2 +- 7 files changed, 142 insertions(+), 160 deletions(-) diff --git a/test/test_centroids.py b/test/test_centroids.py index f609822fd..17556b62f 100644 --- a/test/test_centroids.py +++ b/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, _populate_face_centerpoints +from uxarray.grid.coordinates import _populate_face_centroids, _populate_edge_centroids, _populate_face_centerpoints +from uxarray.grid.utils import _normalize_xyz current_path = Path(os.path.dirname(os.path.realpath(__file__))) diff --git a/test/test_helpers.py b/test/test_helpers.py index 202c51bfe..c5beff4e7 100644 --- a/test/test_helpers.py +++ b/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 diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 381d02e41..d082df709 100644 --- a/uxarray/grid/arcs.py +++ b/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, _normalize_xyz +from uxarray.grid.utils import _xyz_to_lonlat_rad, _normalize_xyz from uxarray.constants import ERROR_TOLERANCE diff --git a/uxarray/grid/area.py b/uxarray/grid/area.py index 1dc6de52e..1503071f7 100644 --- a/uxarray/grid/area.py +++ b/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 diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 485d57593..248df8e7f 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -4,127 +4,14 @@ import warnings -from uxarray.constants import ERROR_TOLERANCE from uxarray.conventions import ugrid - -from typing import Union - -from numba import njit - - -@njit(cache=True) -def _lonlat_rad_to_xyz( - lon: Union[np.ndarray, float], - lat: Union[np.ndarray, float], -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Converts Spherical lon and lat coordinates into Cartesian x, y, z - coordinates.""" - x = np.cos(lon) * np.cos(lat) - y = np.sin(lon) * np.cos(lat) - z = np.sin(lat) - - return x, y, z - - -def _xyz_to_lonlat_rad( - x: Union[np.ndarray, float], - y: Union[np.ndarray, float], - z: Union[np.ndarray, float], - normalize: bool = True, -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Converts Cartesian x, y, z coordinates in Spherical latitude and - longitude coordinates in degrees. - - Parameters - ---------- - x : Union[np.ndarray, float] - Cartesian x coordinates - y: Union[np.ndarray, float] - Cartesiain y coordinates - z: Union[np.ndarray, float] - Cartesian z coordinates - normalize: bool - Flag to select whether to normalize the coordinates - - Returns - ------- - lon : Union[np.ndarray, float] - Longitude in radians - lat: Union[np.ndarray, float] - Latitude in radians - """ - - if normalize: - x, y, z = _normalize_xyz(x, y, z) - denom = np.abs(x * x + y * y + z * z) - x /= denom - y /= denom - z /= denom - - lon = np.arctan2(y, x, dtype=np.float64) - lat = np.arcsin(z, dtype=np.float64) - - # set longitude range to [0, pi] - lon = np.mod(lon, 2 * np.pi) - - z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE - - lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) - lon = np.where(z_mask, 0.0, lon) - - return lon, lat - - -def _xyz_to_lonlat_deg( - x: Union[np.ndarray, float], - y: Union[np.ndarray, float], - z: Union[np.ndarray, float], - normalize: bool = True, -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Converts Cartesian x, y, z coordinates in Spherical latitude and - longitude coordinates in degrees. - - Parameters - ---------- - x : Union[np.ndarray, float] - Cartesian x coordinates - y: Union[np.ndarray, float] - Cartesiain y coordinates - z: Union[np.ndarray, float] - Cartesian z coordinates - normalize: bool - Flag to select whether to normalize the coordinates - - Returns - ------- - lon : Union[np.ndarray, float] - Longitude in degrees - lat: Union[np.ndarray, float] - Latitude in degrees - """ - lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z, normalize=normalize) - - lon = np.rad2deg(lon_rad) - lat = np.rad2deg(lat_rad) - - lon = (lon + 180) % 360 - 180 - return lon, lat - - -def _normalize_xyz( - x: Union[np.ndarray, float], - y: Union[np.ndarray, float], - z: Union[np.ndarray, float], -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Normalizes a set of Cartesiain coordinates.""" - denom = np.linalg.norm( - np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2, axis=0 - ) - - x_norm = x / denom - y_norm = y / denom - z_norm = z / denom - return x_norm, y_norm, z_norm +from uxarray.grid.arcs import _angle_of_2_vectors +from uxarray.grid.utils import ( + _xyz_to_lonlat_rad, + _lonlat_rad_to_xyz, + _xyz_to_lonlat_deg, + _normalize_xyz, +) def _populate_node_latlon(grid) -> None: @@ -322,36 +209,6 @@ def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_fa return _normalize_xyz(centroid_x, centroid_y, centroid_z) -def haversine_distance(point1, point2): - """Calculate the great-circle distance between two points on a unit sphere - using the Haversine formula. - - Parameters - ---------- - point1 : tuple - A tuple containing the lon and lat of the first point. - point2 : tuple - A tuple containing the lon and lat of the second point. - - Returns - ------- - float - The distance between the two points on the unit sphere. - """ - R = 1.0 # Radius of the Earth assumed to be 1 (unit sphere) - lon1, lat1 = np.radians(point1) - lon2, lat2 = np.radians(point2) - - dlat = lat2 - lat1 - dlon = lon2 - lon1 - - a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2 - c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)) - - distance = R * c - return distance - - def circle_from_two_points(p1, p2): """Calculate the smallest circle that encloses two points on a unit sphere. @@ -370,7 +227,9 @@ def circle_from_two_points(p1, p2): center_lon = (p1[0] + p2[0]) / 2 center_lat = (p1[1] + p2[1]) / 2 center = (center_lon, center_lat) - radius = haversine_distance(p1, p2) / 2 + v1, v2 = (np.array(_lonlat_rad_to_xyz(*np.radians(p))) for p in (p1, p2)) + distance = _angle_of_2_vectors(v1, v2) + radius = distance / 2 return center, radius @@ -393,11 +252,12 @@ def circle_from_three_points(p1, p2, p3): A tuple containing the center (as a tuple of lon and lat) and the radius of the circle. """ center = p1 # Placeholder center + v1, v2, v3 = (np.array(_lonlat_rad_to_xyz(*np.radians(p))) for p in (p1, p2, p3)) radius = ( max( - haversine_distance(p1, p2), - haversine_distance(p1, p3), - haversine_distance(p2, p3), + _angle_of_2_vectors(v1, v2), + _angle_of_2_vectors(v1, v3), + _angle_of_2_vectors(v2, v3), ) / 2 ) @@ -420,7 +280,9 @@ def is_inside_circle(circle, point): True if the point is inside the circle, False otherwise. """ center, radius = circle - return haversine_distance(center, point) <= radius + v1, v2 = (np.array(_lonlat_rad_to_xyz(*np.radians(p))) for p in (center, point)) + distance = _angle_of_2_vectors(v1, v2) + return distance <= radius def welzl_recursive(points, boundary, R): diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py index 020aae1d9..e9142ca9a 100644 --- a/uxarray/grid/utils.py +++ b/uxarray/grid/utils.py @@ -3,6 +3,10 @@ import warnings import uxarray.utils.computing as ac_utils +from typing import Union + +from numba import njit + def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None): """Replaces all instances of the current fill value (``original_fill``) in @@ -423,3 +427,118 @@ def _get_lonlat_rad_face_edge_nodes( face_edges_lonlat_rad[valid_mask, 1] = node_lat_rad[valid_edges] return face_edges_lonlat_rad.reshape(n_face, n_max_face_edges, 2, 2) + + +def _xyz_to_lonlat_rad( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], + normalize: bool = True, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Converts Cartesian x, y, z coordinates in Spherical latitude and + longitude coordinates in degrees. + + Parameters + ---------- + x : Union[np.ndarray, float] + Cartesian x coordinates + y: Union[np.ndarray, float] + Cartesiain y coordinates + z: Union[np.ndarray, float] + Cartesian z coordinates + normalize: bool + Flag to select whether to normalize the coordinates + + Returns + ------- + lon : Union[np.ndarray, float] + Longitude in radians + lat: Union[np.ndarray, float] + Latitude in radians + """ + + if normalize: + x, y, z = _normalize_xyz(x, y, z) + denom = np.abs(x * x + y * y + z * z) + x /= denom + y /= denom + z /= denom + + lon = np.arctan2(y, x, dtype=np.float64) + lat = np.arcsin(z, dtype=np.float64) + + # set longitude range to [0, pi] + lon = np.mod(lon, 2 * np.pi) + + z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE + + lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) + lon = np.where(z_mask, 0.0, lon) + + return lon, lat + + +def _normalize_xyz( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Normalizes a set of Cartesiain coordinates.""" + denom = np.linalg.norm( + np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2, axis=0 + ) + + x_norm = x / denom + y_norm = y / denom + z_norm = z / denom + return x_norm, y_norm, z_norm + + +@njit(cache=True) +def _lonlat_rad_to_xyz( + lon: Union[np.ndarray, float], + lat: Union[np.ndarray, float], +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Converts Spherical lon and lat coordinates into Cartesian x, y, z + coordinates.""" + x = np.cos(lon) * np.cos(lat) + y = np.sin(lon) * np.cos(lat) + z = np.sin(lat) + + return x, y, z + + +def _xyz_to_lonlat_deg( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], + normalize: bool = True, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Converts Cartesian x, y, z coordinates in Spherical latitude and + longitude coordinates in degrees. + + Parameters + ---------- + x : Union[np.ndarray, float] + Cartesian x coordinates + y: Union[np.ndarray, float] + Cartesiain y coordinates + z: Union[np.ndarray, float] + Cartesian z coordinates + normalize: bool + Flag to select whether to normalize the coordinates + + Returns + ------- + lon : Union[np.ndarray, float] + Longitude in degrees + lat: Union[np.ndarray, float] + Latitude in degrees + """ + lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z, normalize=normalize) + + lon = np.rad2deg(lon_rad) + lat = np.rad2deg(lat_rad) + + lon = (lon + 180) % 360 - 180 + return lon, lat diff --git a/uxarray/io/_exodus.py b/uxarray/io/_exodus.py index 8c91913fa..c0c5d2778 100644 --- a/uxarray/io/_exodus.py +++ b/uxarray/io/_exodus.py @@ -6,7 +6,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, _xyz_to_lonlat_deg +from uxarray.grid.utils import _lonlat_rad_to_xyz, _xyz_to_lonlat_deg from uxarray.conventions import ugrid From 478f99f54892ea320410a9e2eb5ee4f68c4269da Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Fri, 2 Aug 2024 12:45:07 -0500 Subject: [PATCH 08/29] add benchmark for faceLatLon construction --- benchmarks/mpas_ocean.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index d2baaf86b..7c9a1fd1a 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -164,3 +164,29 @@ def time_nearest_neighbor_remapping(self): def time_inverse_distance_weighted_remapping(self): self.uxds_480["bottomDepth"].remap.inverse_distance_weighted(self.uxds_120.uxgrid) + +from uxarray.grid.coordinates import _construct_face_centerpoints, _construct_face_centroids + +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): + _construct_face_centerpoints(self.uxgrid.node_lon, + self.uxgrid.node_lat, + self.uxgrid.face_nodes, + self.uxgrid.n_nodes_per_face) + + + def time_cartesian_averaging(self, resolution): + _construct_face_centroids(self.uxgrid.node_x, + self.uxgrid.node_y, + self.uxgrid.node_z, + self.uxgrid.face_nodes, + self.uxgrid.n_nodes_per_face) From 277b1ce47308f40bbbf9ce963cfde475c15f80e0 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Fri, 2 Aug 2024 13:19:21 -0500 Subject: [PATCH 09/29] fix face nodes typo --- benchmarks/mpas_ocean.py | 4 ++-- uxarray/grid/coordinates.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index 7c9a1fd1a..cf220fdc3 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -180,7 +180,7 @@ def teardown(self, resolution): def time_welzl(self, resolution): _construct_face_centerpoints(self.uxgrid.node_lon, self.uxgrid.node_lat, - self.uxgrid.face_nodes, + self.uxgrid.face_node_connectivity, self.uxgrid.n_nodes_per_face) @@ -188,5 +188,5 @@ def time_cartesian_averaging(self, resolution): _construct_face_centroids(self.uxgrid.node_x, self.uxgrid.node_y, self.uxgrid.node_z, - self.uxgrid.face_nodes, + self.uxgrid.face_node_connectivity, self.uxgrid.n_nodes_per_face) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 248df8e7f..2a96b1aa8 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -359,7 +359,7 @@ def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_fac Latitudes of the nodes. face_nodes : array_like Indices of nodes per face. - n_nodes_per_face : int + n_nodes_per_face : array_like Number of nodes per face. Returns From 68edd31a9653f4d92210afde9bf5f97052d50fdb Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Fri, 2 Aug 2024 13:54:03 -0500 Subject: [PATCH 10/29] fix face nodes typo --- benchmarks/mpas_ocean.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index cf220fdc3..9e81dd8ed 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -178,15 +178,15 @@ def teardown(self, resolution): del self.uxgrid def time_welzl(self, resolution): - _construct_face_centerpoints(self.uxgrid.node_lon, - self.uxgrid.node_lat, - self.uxgrid.face_node_connectivity, - self.uxgrid.n_nodes_per_face) + _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): - _construct_face_centroids(self.uxgrid.node_x, - self.uxgrid.node_y, - self.uxgrid.node_z, - self.uxgrid.face_node_connectivity, - self.uxgrid.n_nodes_per_face) + _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) From c580465734f9f907ca099be19b8dee81611922cc Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Tue, 6 Aug 2024 05:43:58 -0500 Subject: [PATCH 11/29] o Remove new properties _ctrpt infavor of using regular face lat/lon x/y/z for welzl centerpoint also: introduce construct_face_center method in Grid --- docs/user_api/index.rst | 1 + test/test_centroids.py | 12 +++--- uxarray/grid/coordinates.py | 25 ++++++------ uxarray/grid/grid.py | 76 +++++++++++-------------------------- 4 files changed, 41 insertions(+), 73 deletions(-) diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst index ac521c048..0a939e310 100644 --- a/docs/user_api/index.rst +++ b/docs/user_api/index.rst @@ -224,6 +224,7 @@ Methods Grid.get_kd_tree Grid.copy Grid.isel + Grid.compute_face_center Dimensions diff --git a/test/test_centroids.py b/test/test_centroids.py index 17556b62f..aed5ceb38 100644 --- a/test/test_centroids.py +++ b/test/test_centroids.py @@ -114,12 +114,14 @@ def test_face_centerpoint(self): 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) - _populate_face_centerpoints(uxgrid) - # the expected centerpoint should be close to centroid for this case + # 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] - # Test the values of the calculated centerpoint, giving high tolerance of one decimal place - nt.assert_array_almost_equal(ctr_lon, uxgrid.face_lon_ctrpt.values[0], decimal=1) - nt.assert_array_almost_equal(ctr_lat, uxgrid.face_lat_ctrpt.values[0], decimal=1) + # now explicitly get the centerpoints stored to face_lon/lat using welzl's centerpoint algorithm + uxgrid.compute_face_center(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) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 2a96b1aa8..a9a0e68eb 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -141,34 +141,32 @@ def _populate_face_centerpoints(grid, repopulate=False): n_nodes_per_face = grid.n_nodes_per_face.values # Check if the centerpoints are already populated - if "face_lon_ctrpt" not in grid._ds or repopulate: - # Construct the centerpoints if there are none stored - if "face_x_ctrpt" not in grid._ds: - centerpoint_lon, centerpoint_lat = _construct_face_centerpoints( - node_lon, node_lat, face_nodes, n_nodes_per_face - ) + if "face_lon" not in grid._ds or repopulate: + centerpoint_lon, centerpoint_lat = _construct_face_centerpoints( + node_lon, node_lat, face_nodes, n_nodes_per_face + ) # get the cartesian coordinates of the centerpoints ctrpt_x, ctrpt_y, ctrpt_z = _lonlat_rad_to_xyz(centerpoint_lon, centerpoint_lat) # set the grid variables for centerpoints - if "face_lon_ctrpt" not in grid._ds or repopulate: - grid._ds["face_lon_ctrpt"] = xr.DataArray( + if "face_lon" not in grid._ds or repopulate: + grid._ds["face_lon"] = xr.DataArray( centerpoint_lon, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LON_ATTRS ) - grid._ds["face_lat_ctrpt"] = xr.DataArray( + grid._ds["face_lat"] = xr.DataArray( centerpoint_lat, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LAT_ATTRS ) - if "face_x_ctrpt" not in grid._ds or repopulate: - grid._ds["face_x_ctrpt"] = xr.DataArray( + if "face_x" not in grid._ds or repopulate: + grid._ds["face_x"] = xr.DataArray( ctrpt_x, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_X_ATTRS ) - grid._ds["face_y_ctrpt"] = xr.DataArray( + grid._ds["face_y"] = xr.DataArray( ctrpt_y, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Y_ATTRS ) - grid._ds["face_z_ctrpt"] = xr.DataArray( + grid._ds["face_z"] = xr.DataArray( ctrpt_z, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Z_ATTRS ) @@ -195,7 +193,6 @@ def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_fa tuple The x, y, and z coordinates of the centroids. """ - centroid_x = np.zeros((face_nodes.shape[0]), dtype=np.float64) centroid_y = np.zeros((face_nodes.shape[0]), dtype=np.float64) centroid_z = np.zeros((face_nodes.shape[0]), dtype=np.float64) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index a9f7bc1a0..f76a96075 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -345,6 +345,28 @@ def validate(self): else: raise RuntimeError("Mesh validation failed.") + def compute_face_center(self, method="average"): + """Constructs a face_lon and face_lat. + + Parameters + ---------- + method : str, default="average" + other supported method is Welzl's algorithm, takes value "welzl" + + + Usage + ----- + >>> import uxarray as ux + >>> uxgrid = ux.open_grid("GRID_FILE_NAME") + >>> face_lat = uxgrid.construct_face_center(method="welzl") + """ + if method == "average": + _populate_face_centroids(self, repopulate=True) + elif method == "welzl": + _populate_face_centerpoints(self, repopulate=True) + else: + raise ValueError("unknown method for face center calculation") + def __repr__(self): """Constructs a string representation of the contents of a ``Grid``.""" @@ -731,60 +753,6 @@ def face_z(self) -> xr.DataArray: _populate_face_centroids(self) return self._ds["face_z"] - @property - def face_lon_ctrpt(self) -> xr.DataArray: - """Longitude of the center of each face in degrees. - - Dimensions: ``(n_face, )`` - """ - if "face_lon_ctrpt" not in self._ds: - _populate_face_centerpoints(self) - _set_desired_longitude_range(self._ds) - return self._ds["face_lon_ctrpt"] - - @property - def face_lat_ctrpt(self) -> xr.DataArray: - """Latitude of the center of each face in degrees. - - Dimensions: ``(n_face, )`` - """ - if "face_lat_ctrpt" not in self._ds: - _populate_face_centerpoints(self) - _set_desired_longitude_range(self._ds) - - return self._ds["face_lat_ctrpt"] - - @property - def face_x_ctrpt(self) -> xr.DataArray: - """Cartesian x location of the center of each face in meters. - - Dimensions: ``(n_face, )`` - """ - if "face_x_ctrpt" not in self._ds: - _populate_face_centerpoints(self) - - return self._ds["face_x_ctrpt"] - - @property - def face_y_ctrpt(self) -> xr.DataArray: - """Cartesian y location of the center of each face in meters. - - Dimensions: ``(n_face, )`` - """ - if "face_y_ctrpt" not in self._ds: - _populate_face_centerpoints(self) - return self._ds["face_y_ctrpt"] - - @property - def face_z_ctrpt(self) -> xr.DataArray: - """Cartesian z location of the center of each face in meters. - - Dimensions: ``(n_face, )`` - """ - if "face_z_ctrpt" not in self._ds: - _populate_face_centerpoints(self) - return self._ds["face_z_ctrpt"] - @property def face_node_connectivity(self) -> xr.DataArray: """Indices of the nodes that make up each face. From b73fb31eaac77e6ab6c07df1c5313d8ef1155dd8 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Mon, 26 Aug 2024 18:08:32 -0500 Subject: [PATCH 12/29] o precommit --- uxarray/grid/coordinates.py | 2 +- uxarray/grid/utils.py | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 878479b23..227a33314 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -8,12 +8,12 @@ from uxarray.grid.arcs import _angle_of_2_vectors from uxarray.grid.utils import ( _xyz_to_lonlat_rad, - _xyz_to_lonlat_rad_no_norm, _lonlat_rad_to_xyz, _xyz_to_lonlat_deg, _normalize_xyz, ) + def _populate_node_latlon(grid) -> None: """Populates the lon and lat coordinates of a Grid (`node_lon`, `node_lat`)""" diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py index c1dac9090..18bc171e9 100644 --- a/uxarray/grid/utils.py +++ b/uxarray/grid/utils.py @@ -2,7 +2,7 @@ from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE import warnings import uxarray.utils.computing as ac_utils -import math +import math from typing import Union @@ -478,6 +478,7 @@ def _xyz_to_lonlat_rad( return lon, lat + @njit def _xyz_to_lonlat_rad_no_norm( x: Union[np.ndarray, float], @@ -548,6 +549,7 @@ def _lonlat_rad_to_xyz( return x, y, z + def _xyz_to_lonlat_deg( x: Union[np.ndarray, float], y: Union[np.ndarray, float], @@ -583,6 +585,7 @@ def _xyz_to_lonlat_deg( lon = (lon + 180) % 360 - 180 return lon, lat + @njit def _normalize_xyz_scalar(x: float, y: float, z: float): denom = np.linalg.norm(np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2) From 4f65cb89f42dfcbd4ce086fd47296c3d31bcf2da Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Mon, 26 Aug 2024 18:31:08 -0500 Subject: [PATCH 13/29] o reduce some python-level loop overhead, use zip and some list comprehension in _construct_face_centerpoints --- uxarray/grid/coordinates.py | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 227a33314..c04d8c2bc 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -363,27 +363,29 @@ def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_fac ------- tuple of numpy.ndarray Two arrays containing the longitudes and latitudes of the centerpoints. - - Notes - ----- - This function calculates the centerpoints of faces defined by nodes on a sphere, using Welzl's algorithm to - find the smallest enclosing circle for each face. """ - ctrpt_lon = np.zeros(face_nodes.shape[0], dtype=np.float64) - ctrpt_lat = np.zeros(face_nodes.shape[0], dtype=np.float64) + num_faces = face_nodes.shape[0] + ctrpt_lon = np.zeros(num_faces, dtype=np.float64) + ctrpt_lat = np.zeros(num_faces, dtype=np.float64) - for face_idx, n_max_nodes in enumerate(n_nodes_per_face): - points_array = np.column_stack( + # Pre-compute all points arrays + points_arrays = [ + np.column_stack( ( node_lon[face_nodes[face_idx, :n_max_nodes]], node_lat[face_nodes[face_idx, :n_max_nodes]], ) ) - circle = smallest_enclosing_circle(points_array) - ctrpt_lon[face_idx] = circle[0][0] - ctrpt_lat[face_idx] = circle[0][1] + for face_idx, n_max_nodes in enumerate(n_nodes_per_face) + ] + + # Compute circles for all faces + circles = [smallest_enclosing_circle(points) for points in points_arrays] + + # Extract centerpoints + ctrpt_lon, ctrpt_lat = zip(*[circle[0] for circle in circles]) - return ctrpt_lon, ctrpt_lat + return np.array(ctrpt_lon), np.array(ctrpt_lat) def _populate_edge_centroids(grid, repopulate=False): From f5bd37b586df0995384fd7b4d5e8b3e44250ca12 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Tue, 27 Aug 2024 00:12:53 -0500 Subject: [PATCH 14/29] o Fix structure, numba causes issues with recursive functions, but this is much more efficient now. If need be we can later change to iterative instead of recursion for avoiding the error: Numba doesn't support the use of yield in closures. --- uxarray/grid/coordinates.py | 71 ++++++++++++++++++------------------- 1 file changed, 35 insertions(+), 36 deletions(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index c04d8c2bc..5bb32da19 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -3,7 +3,6 @@ import warnings - from uxarray.conventions import ugrid from uxarray.grid.arcs import _angle_of_2_vectors from uxarray.grid.utils import ( @@ -115,6 +114,41 @@ def _populate_face_centroids(grid, repopulate=False): ) +def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face): + """Constructs the xyz centroid coordinate for each face using Cartesian + Averaging. + + Parameters + ---------- + node_x : numpy.ndarray + X coordinates of the nodes. + node_y : numpy.ndarray + Y coordinates of the nodes. + node_z : numpy.ndarray + Z coordinates of the nodes. + face_nodes : numpy.ndarray + Indices of nodes per face. + n_nodes_per_face : numpy.ndarray + Number of nodes per face. + + Returns + ------- + tuple + The x, y, and z coordinates of the centroids. + """ + centroid_x = np.zeros((face_nodes.shape[0]), dtype=np.float64) + centroid_y = np.zeros((face_nodes.shape[0]), dtype=np.float64) + centroid_z = np.zeros((face_nodes.shape[0]), dtype=np.float64) + + for face_idx, n_max_nodes in enumerate(n_nodes_per_face): + # Compute Cartesian Average + centroid_x[face_idx] = np.mean(node_x[face_nodes[face_idx, 0:n_max_nodes]]) + centroid_y[face_idx] = np.mean(node_y[face_nodes[face_idx, 0:n_max_nodes]]) + centroid_z[face_idx] = np.mean(node_z[face_nodes[face_idx, 0:n_max_nodes]]) + + return _normalize_xyz(centroid_x, centroid_y, centroid_z) + + def _populate_face_centerpoints(grid, repopulate=False): """Calculates the face centerpoints using Welzl's algorithm. It is a randomized algorithm for finding the center and radius of the smallest @@ -171,41 +205,6 @@ def _populate_face_centerpoints(grid, repopulate=False): ) -def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face): - """Constructs the xyz centroid coordinate for each face using Cartesian - Averaging. - - Parameters - ---------- - node_x : numpy.ndarray - X coordinates of the nodes. - node_y : numpy.ndarray - Y coordinates of the nodes. - node_z : numpy.ndarray - Z coordinates of the nodes. - face_nodes : numpy.ndarray - Indices of nodes per face. - n_nodes_per_face : numpy.ndarray - Number of nodes per face. - - Returns - ------- - tuple - The x, y, and z coordinates of the centroids. - """ - centroid_x = np.zeros((face_nodes.shape[0]), dtype=np.float64) - centroid_y = np.zeros((face_nodes.shape[0]), dtype=np.float64) - centroid_z = np.zeros((face_nodes.shape[0]), dtype=np.float64) - - for face_idx, n_max_nodes in enumerate(n_nodes_per_face): - # Compute Cartesian Average - centroid_x[face_idx] = np.mean(node_x[face_nodes[face_idx, 0:n_max_nodes]]) - centroid_y[face_idx] = np.mean(node_y[face_nodes[face_idx, 0:n_max_nodes]]) - centroid_z[face_idx] = np.mean(node_z[face_nodes[face_idx, 0:n_max_nodes]]) - - return _normalize_xyz(centroid_x, centroid_y, centroid_z) - - def circle_from_two_points(p1, p2): """Calculate the smallest circle that encloses two points on a unit sphere. From 0be7cfda9625ef37cef54e22951853a482c828ae Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Tue, 27 Aug 2024 02:21:54 -0500 Subject: [PATCH 15/29] o Modify the recursion logic for Numba compatibility, works 3-10x faster --- uxarray/grid/coordinates.py | 43 ++++++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 10 deletions(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 5bb32da19..b5914c2f5 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -11,6 +11,7 @@ _xyz_to_lonlat_deg, _normalize_xyz, ) +from numba import njit def _populate_node_latlon(grid) -> None: @@ -205,6 +206,7 @@ def _populate_face_centerpoints(grid, repopulate=False): ) +@njit def circle_from_two_points(p1, p2): """Calculate the smallest circle that encloses two points on a unit sphere. @@ -223,12 +225,17 @@ def circle_from_two_points(p1, p2): center_lon = (p1[0] + p2[0]) / 2 center_lat = (p1[1] + p2[1]) / 2 center = (center_lon, center_lat) - v1, v2 = (np.array(_lonlat_rad_to_xyz(*np.radians(p))) for p in (p1, p2)) + + v1 = np.array(_lonlat_rad_to_xyz(np.radians(p1[0]), np.radians(p1[1]))) + v2 = np.array(_lonlat_rad_to_xyz(np.radians(p2[0]), np.radians(p2[1]))) + distance = _angle_of_2_vectors(v1, v2) radius = distance / 2 + return center, radius +@njit def circle_from_three_points(p1, p2, p3): """Calculate the smallest circle that encloses three points on a unit sphere. This is a placeholder implementation. @@ -247,8 +254,15 @@ def circle_from_three_points(p1, p2, p3): tuple A tuple containing the center (as a tuple of lon and lat) and the radius of the circle. """ - center = p1 # Placeholder center - v1, v2, v3 = (np.array(_lonlat_rad_to_xyz(*np.radians(p))) for p in (p1, p2, p3)) + # Placeholder implementation for three-point circle calculation + center_lon = (p1[0] + p2[0] + p3[0]) / 3 + center_lat = (p1[1] + p2[1] + p3[1]) / 3 + center = (center_lon, center_lat) + + v1 = np.array(_lonlat_rad_to_xyz(np.radians(p1[0]), np.radians(p1[1]))) + v2 = np.array(_lonlat_rad_to_xyz(np.radians(p2[0]), np.radians(p2[1]))) + v3 = np.array(_lonlat_rad_to_xyz(np.radians(p3[0]), np.radians(p3[1]))) + radius = ( max( _angle_of_2_vectors(v1, v2), @@ -257,9 +271,11 @@ def circle_from_three_points(p1, p2, p3): ) / 2 ) + return center, radius +@njit def is_inside_circle(circle, point): """Check if a point is inside a given circle on a unit sphere. @@ -276,11 +292,13 @@ def is_inside_circle(circle, point): True if the point is inside the circle, False otherwise. """ center, radius = circle - v1, v2 = (np.array(_lonlat_rad_to_xyz(*np.radians(p))) for p in (center, point)) + v1 = np.array(_lonlat_rad_to_xyz(np.radians(center[0]), np.radians(center[1]))) + v2 = np.array(_lonlat_rad_to_xyz(np.radians(point[0]), np.radians(point[1]))) distance = _angle_of_2_vectors(v1, v2) return distance <= radius +@njit def welzl_recursive(points, boundary, R): """Recursive helper function for Welzl's algorithm to find the smallest enclosing circle. @@ -303,25 +321,30 @@ def welzl_recursive(points, boundary, R): if len(points) == 0 or len(boundary) == 3: # Construct the minimal circle based on the number of boundary points if len(boundary) == 0: - return R + # Return a default circle if no boundary points are available + return ((0.0, 0.0), 0.0) elif len(boundary) == 1: - return (boundary[0], 0) + return circle_from_two_points(boundary[0], boundary[0]) elif len(boundary) == 2: return circle_from_two_points(boundary[0], boundary[1]) elif len(boundary) == 3: return circle_from_three_points(boundary[0], boundary[1], boundary[2]) - # Choose a point from the set and remove it p = points[-1] - temp_points = np.delete(points, -1, axis=0) + temp_points = points[:-1] circle = welzl_recursive(temp_points, boundary, R) # Check if the chosen point is inside the current circle if circle and is_inside_circle(circle, p): return circle + # If not, the point must be on the boundary of the minimal enclosing circle else: - # If not, the point must be on the boundary of the minimal enclosing circle - return welzl_recursive(temp_points, np.append(boundary, [p], axis=0), R) + new_boundary = np.empty( + (boundary.shape[0] + 1, boundary.shape[1]), dtype=boundary.dtype + ) + new_boundary[:-1] = boundary + new_boundary[-1] = p + return welzl_recursive(temp_points, new_boundary, R) def smallest_enclosing_circle(points): From d77e7262fb4fde593cf13e56352d1ef934a800b3 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Tue, 27 Aug 2024 09:53:25 -0500 Subject: [PATCH 16/29] o Try to get the codecov higher, add internal function tests --- benchmarks/quad_hexagon.py | 36 +++++++++++++++++++++++++++++ docs/internal_api/index.rst | 18 ++++++++++----- test/test_centroids.py | 45 +++++++++++++++++++++++++++++++++++-- uxarray/grid/coordinates.py | 26 ++++++++++----------- 4 files changed, 105 insertions(+), 20 deletions(-) diff --git a/benchmarks/quad_hexagon.py b/benchmarks/quad_hexagon.py index 4364f39b0..aeb7559bb 100644 --- a/benchmarks/quad_hexagon.py +++ b/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) diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst index b263f4ec8..e6391c5d0 100644 --- a/docs/internal_api/index.rst +++ b/docs/internal_api/index.rst @@ -147,10 +147,6 @@ 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 @@ -158,6 +154,13 @@ Coordinates 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 @@ -179,7 +182,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 diff --git a/test/test_centroids.py b/test/test_centroids.py index aed5ceb38..e97e1a30c 100644 --- a/test/test_centroids.py +++ b/test/test_centroids.py @@ -4,7 +4,7 @@ import numpy.testing as nt import uxarray as ux from pathlib import Path -from uxarray.grid.coordinates import _populate_face_centroids, _populate_edge_centroids, _populate_face_centerpoints +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__))) @@ -64,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.compute_face_center(method="average") # computed_face_x = (uxgrid.face_lon.values + 180) % 360 - 180 computed_face_x = uxgrid.face_lon.values @@ -109,6 +110,46 @@ 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.""" diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index b5914c2f5..c6f617517 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -207,7 +207,7 @@ def _populate_face_centerpoints(grid, repopulate=False): @njit -def circle_from_two_points(p1, p2): +def _circle_from_two_points(p1, p2): """Calculate the smallest circle that encloses two points on a unit sphere. Parameters @@ -236,7 +236,7 @@ def circle_from_two_points(p1, p2): @njit -def circle_from_three_points(p1, p2, p3): +def _circle_from_three_points(p1, p2, p3): """Calculate the smallest circle that encloses three points on a unit sphere. This is a placeholder implementation. @@ -276,7 +276,7 @@ def circle_from_three_points(p1, p2, p3): @njit -def is_inside_circle(circle, point): +def _is_inside_circle(circle, point): """Check if a point is inside a given circle on a unit sphere. Parameters @@ -299,7 +299,7 @@ def is_inside_circle(circle, point): @njit -def welzl_recursive(points, boundary, R): +def _welzl_recursive(points, boundary, R): """Recursive helper function for Welzl's algorithm to find the smallest enclosing circle. @@ -324,18 +324,18 @@ def welzl_recursive(points, boundary, R): # Return a default circle if no boundary points are available return ((0.0, 0.0), 0.0) elif len(boundary) == 1: - return circle_from_two_points(boundary[0], boundary[0]) + return _circle_from_two_points(boundary[0], boundary[0]) elif len(boundary) == 2: - return circle_from_two_points(boundary[0], boundary[1]) + return _circle_from_two_points(boundary[0], boundary[1]) elif len(boundary) == 3: - return circle_from_three_points(boundary[0], boundary[1], boundary[2]) + return _circle_from_three_points(boundary[0], boundary[1], boundary[2]) p = points[-1] temp_points = points[:-1] - circle = welzl_recursive(temp_points, boundary, R) + circle = _welzl_recursive(temp_points, boundary, R) # Check if the chosen point is inside the current circle - if circle and is_inside_circle(circle, p): + if circle and _is_inside_circle(circle, p): return circle # If not, the point must be on the boundary of the minimal enclosing circle else: @@ -344,10 +344,10 @@ def welzl_recursive(points, boundary, R): ) new_boundary[:-1] = boundary new_boundary[-1] = p - return welzl_recursive(temp_points, new_boundary, R) + return _welzl_recursive(temp_points, new_boundary, R) -def smallest_enclosing_circle(points): +def _smallest_enclosing_circle(points): """Find the smallest circle that encloses all given points on a unit sphere using Welzl's algorithm. @@ -364,7 +364,7 @@ def smallest_enclosing_circle(points): np.random.shuffle( points ) # Randomize the input to increase the chance of an optimal solution - return welzl_recursive(points, np.empty((0, 2)), None) + return _welzl_recursive(points, np.empty((0, 2)), None) def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_face): @@ -402,7 +402,7 @@ def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_fac ] # Compute circles for all faces - circles = [smallest_enclosing_circle(points) for points in points_arrays] + circles = [_smallest_enclosing_circle(points) for points in points_arrays] # Extract centerpoints ctrpt_lon, ctrpt_lat = zip(*[circle[0] for circle in circles]) From 1a185307acbcdc3a191591e2fc58ac60a487dd58 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Wed, 4 Sep 2024 00:24:05 -0500 Subject: [PATCH 17/29] o Try to include inside funcs --- benchmarks/mpas_ocean.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index 9e81dd8ed..b2a8da300 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -165,9 +165,9 @@ def time_nearest_neighbor_remapping(self): def time_inverse_distance_weighted_remapping(self): self.uxds_480["bottomDepth"].remap.inverse_distance_weighted(self.uxds_120.uxgrid) -from uxarray.grid.coordinates import _construct_face_centerpoints, _construct_face_centroids class ConstructFaceLatLon: + param_names = ['resolution'] params = ['480km', '120km'] @@ -178,6 +178,9 @@ 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, @@ -185,6 +188,8 @@ def time_welzl(self, resolution): 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, From 67e3de3d6884d1ca4e11a896116502d7589a5b6a Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Tue, 10 Sep 2024 11:56:50 -0500 Subject: [PATCH 18/29] o call cartesian average instead of average add more doc --- docs/user_api/index.rst | 2 +- test/test_centroids.py | 4 ++-- uxarray/grid/grid.py | 19 +++++++++++++------ 3 files changed, 16 insertions(+), 9 deletions(-) diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst index 33061a859..c35d9c01b 100644 --- a/docs/user_api/index.rst +++ b/docs/user_api/index.rst @@ -224,7 +224,7 @@ Methods Grid.get_kd_tree Grid.copy Grid.isel - Grid.compute_face_center + Grid.construct_face_centers Grid.chunk diff --git a/test/test_centroids.py b/test/test_centroids.py index e97e1a30c..eefada5c5 100644 --- a/test/test_centroids.py +++ b/test/test_centroids.py @@ -65,7 +65,7 @@ def test_centroids_from_mean_verts_scrip(self): expected_face_y = uxgrid.face_lat.values # _populate_face_centroids(uxgrid, repopulate=True) - uxgrid.compute_face_center(method="average") + uxgrid.construct_face_centers(method="cartesian average") # computed_face_x = (uxgrid.face_lon.values + 180) % 360 - 180 computed_face_x = uxgrid.face_lon.values @@ -161,7 +161,7 @@ def test_face_centerpoint(self): ctr_lat = uxgrid.face_lat.values[0] # now explicitly get the centerpoints stored to face_lon/lat using welzl's centerpoint algorithm - uxgrid.compute_face_center(method = "welzl") + 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) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 42f82507f..c47b1ead4 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -348,13 +348,18 @@ def validate(self): else: raise RuntimeError("Mesh validation failed.") - def compute_face_center(self, method="average"): - """Constructs a face_lon and face_lat. + def construct_face_centers(self, method="cartesian average"): + """Constructs face centers, this method provides users direct control + of the method for constructing the face centers, the default method is + "cartesian average", but a more efficient method is "welzl" that is + based on the recursive Welzl algorithm. It must be noted that this + method can override the parsed/recompute the original parsed face + centers. Parameters ---------- - method : str, default="average" - other supported method is Welzl's algorithm, takes value "welzl" + method : str, default="cartesian average" + Supported methods are "cartesian average" and "welzl" Usage @@ -363,12 +368,14 @@ def compute_face_center(self, method="average"): >>> uxgrid = ux.open_grid("GRID_FILE_NAME") >>> face_lat = uxgrid.construct_face_center(method="welzl") """ - if method == "average": + if method == "cartesian average": _populate_face_centroids(self, repopulate=True) elif method == "welzl": _populate_face_centerpoints(self, repopulate=True) else: - raise ValueError("unknown method for face center calculation") + raise ValueError( + f"Unknown method for face center calculation. Expected one of ['cartesian average', 'welzl'] but received {method}" + ) def __repr__(self): """Constructs a string representation of the contents of a ``Grid``.""" From fcc22c1dba26b3112f39312878bb09a667d171d0 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 12 Sep 2024 18:35:28 -0500 Subject: [PATCH 19/29] o Add return doc --- uxarray/grid/coordinates.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index c6f617517..f781f34d2 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -163,6 +163,10 @@ def _populate_face_centerpoints(grid, repopulate=False): The grid containing the nodes and faces. repopulate : bool, optional Bool used to turn on/off repopulating the face coordinates of the centerpoints, default is False. + + Returns + ------- + None, populates the grid with the face centerpoints: face_lon, face_lat """ # warnings.warn("This cannot be guaranteed to work correctly on concave polygons") From ac8e099d6f3238a7155096edf832edc705b76b94 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Fri, 13 Sep 2024 10:05:41 -0500 Subject: [PATCH 20/29] o Add return test for doc --- uxarray/grid/grid.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index abd0e113d..c5602a089 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -362,6 +362,10 @@ def construct_face_centers(self, method="cartesian average"): method : str, default="cartesian average" Supported methods are "cartesian average" and "welzl" + Returns + ------- + None + This method constructs the face_lon and face_lat attributes for the grid object. Usage ----- From dba53dc8cc6abd89571150bc06b1421554dc35e4 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Mon, 16 Sep 2024 09:44:02 -0500 Subject: [PATCH 21/29] o fix text --- uxarray/grid/grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index c5602a089..c4d134dab 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -352,7 +352,7 @@ def validate(self): def construct_face_centers(self, method="cartesian average"): """Constructs face centers, this method provides users direct control of the method for constructing the face centers, the default method is - "cartesian average", but a more efficient method is "welzl" that is + "cartesian average", but a more accurate method is "welzl" that is based on the recursive Welzl algorithm. It must be noted that this method can override the parsed/recompute the original parsed face centers. From ecc0fa032d086884afe0513837d56180aeb4c57b Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Mon, 16 Sep 2024 18:20:58 -0500 Subject: [PATCH 22/29] o Introduce circular dependency issue --- test/test_helpers.py | 2 +- uxarray/grid/arcs.py | 2 +- uxarray/grid/coordinates.py | 247 +++++++++++++++++++++++++++++++++++- uxarray/grid/utils.py | 234 ---------------------------------- uxarray/io/_exodus.py | 2 +- 5 files changed, 243 insertions(+), 244 deletions(-) diff --git a/test/test_helpers.py b/test/test_helpers.py index de6c809f7..c5b923a26 100644 --- a/test/test_helpers.py +++ b/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.utils import _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad +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 diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index d704bc1a2..63344473e 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -2,7 +2,7 @@ # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place -from uxarray.grid.utils import _xyz_to_lonlat_rad_no_norm, _normalize_xyz_scalar +from uxarray.grid.coordinates import _xyz_to_lonlat_rad_no_norm, _normalize_xyz_scalar from uxarray.constants import ERROR_TOLERANCE from uxarray.utils.computing import isclose, cross, dot diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index f781f34d2..3ecf51592 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -1,17 +1,21 @@ import xarray as xr import numpy as np - +import math import warnings from uxarray.conventions import ugrid from uxarray.grid.arcs import _angle_of_2_vectors -from uxarray.grid.utils import ( - _xyz_to_lonlat_rad, - _lonlat_rad_to_xyz, - _xyz_to_lonlat_deg, - _normalize_xyz, -) + +# from uxarray.grid.utils import ( +# _xyz_to_lonlat_rad, +# _lonlat_rad_to_xyz, +# _xyz_to_lonlat_deg, +# _normalize_xyz, +# ) from numba import njit +from uxarray.constants import ERROR_TOLERANCE + +from typing import Union def _populate_node_latlon(grid) -> None: @@ -501,3 +505,232 @@ def _set_desired_longitude_range(ds): if lon_name in ds: if ds[lon_name].max() > 180: ds[lon_name].data = (ds[lon_name].data + 180) % 360 - 180 + + +def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None): + """Replaces all instances of the current fill value (``original_fill``) in + (``grid_var``) with (``new_fill``) and converts to the dtype defined by + (``new_dtype``) + + Parameters + ---------- + grid_var : np.ndarray + grid variable to be modified + original_fill : constant + original fill value used in (``grid_var``) + new_fill : constant + new fill value to be used in (``grid_var``) + new_dtype : np.dtype, optional + new data type to convert (``grid_var``) to + + Returns + ---------- + grid_var : xarray.Dataset + Input Dataset with correct fill value and dtype + """ + + # locations of fill values + if original_fill is not None and np.isnan(original_fill): + fill_val_idx = np.isnan(grid_var) + else: + fill_val_idx = grid_var == original_fill + + # convert to new data type + if new_dtype != grid_var.dtype and new_dtype is not None: + grid_var = grid_var.astype(new_dtype) + + # ensure fill value can be represented with current integer data type + if np.issubdtype(new_dtype, np.integer): + int_min = np.iinfo(grid_var.dtype).min + int_max = np.iinfo(grid_var.dtype).max + # ensure new_fill is in range [int_min, int_max] + if new_fill < int_min or new_fill > int_max: + raise ValueError( + f"New fill value: {new_fill} not representable by" + f" integer dtype: {grid_var.dtype}" + ) + + # ensure non-nan fill value can be represented with current float data type + elif np.issubdtype(new_dtype, np.floating) and not np.isnan(new_fill): + float_min = np.finfo(grid_var.dtype).min + float_max = np.finfo(grid_var.dtype).max + # ensure new_fill is in range [float_min, float_max] + if new_fill < float_min or new_fill > float_max: + raise ValueError( + f"New fill value: {new_fill} not representable by" + f" float dtype: {grid_var.dtype}" + ) + else: + raise ValueError( + f"Data type {grid_var.dtype} not supported" f"for grid variables" + ) + + # replace all zeros with a fill value + grid_var[fill_val_idx] = new_fill + + return grid_var + + +def _xyz_to_lonlat_rad( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], + normalize: bool = True, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Converts Cartesian x, y, z coordinates in Spherical latitude and + longitude coordinates in degrees. + + Parameters + ---------- + x : Union[np.ndarray, float] + Cartesian x coordinates + y: Union[np.ndarray, float] + Cartesiain y coordinates + z: Union[np.ndarray, float] + Cartesian z coordinates + normalize: bool + Flag to select whether to normalize the coordinates + + Returns + ------- + lon : Union[np.ndarray, float] + Longitude in radians + lat: Union[np.ndarray, float] + Latitude in radians + """ + + if normalize: + x, y, z = _normalize_xyz(x, y, z) + denom = np.abs(x * x + y * y + z * z) + x /= denom + y /= denom + z /= denom + + lon = np.arctan2(y, x, dtype=np.float64) + lat = np.arcsin(z, dtype=np.float64) + + # set longitude range to [0, pi] + lon = np.mod(lon, 2 * np.pi) + + z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE + + lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) + lon = np.where(z_mask, 0.0, lon) + + return lon, lat + + +@njit +def _xyz_to_lonlat_rad_no_norm( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], +): + """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 + """ + + lon = math.atan2(y, x) + lat = math.asin(z) + + # set longitude range to [0, pi] + lon = np.mod(lon, 2 * np.pi) + + z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE + + lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) + lon = np.where(z_mask, 0.0, lon) + + return lon, lat + + +def _normalize_xyz( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Normalizes a set of Cartesiain coordinates.""" + denom = np.linalg.norm( + np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2, axis=0 + ) + + x_norm = x / denom + y_norm = y / denom + z_norm = z / denom + return x_norm, y_norm, z_norm + + +@njit(cache=True) +def _lonlat_rad_to_xyz( + lon: Union[np.ndarray, float], + lat: Union[np.ndarray, float], +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Converts Spherical lon and lat coordinates into Cartesian x, y, z + coordinates.""" + x = np.cos(lon) * np.cos(lat) + y = np.sin(lon) * np.cos(lat) + z = np.sin(lat) + + return x, y, z + + +def _xyz_to_lonlat_deg( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], + normalize: bool = True, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Converts Cartesian x, y, z coordinates in Spherical latitude and + longitude coordinates in degrees. + + Parameters + ---------- + x : Union[np.ndarray, float] + Cartesian x coordinates + y: Union[np.ndarray, float] + Cartesiain y coordinates + z: Union[np.ndarray, float] + Cartesian z coordinates + normalize: bool + Flag to select whether to normalize the coordinates + + Returns + ------- + lon : Union[np.ndarray, float] + Longitude in degrees + lat: Union[np.ndarray, float] + Latitude in degrees + """ + lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z, normalize=normalize) + + lon = np.rad2deg(lon_rad) + lat = np.rad2deg(lat_rad) + + lon = (lon + 180) % 360 - 180 + return lon, lat + + +@njit +def _normalize_xyz_scalar(x: float, y: float, z: float): + denom = np.linalg.norm(np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2) + x_norm = x / denom + y_norm = y / denom + z_norm = z / denom + return x_norm, y_norm, z_norm diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py index 18bc171e9..64acddba4 100644 --- a/uxarray/grid/utils.py +++ b/uxarray/grid/utils.py @@ -2,75 +2,6 @@ from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE import warnings import uxarray.utils.computing as ac_utils -import math - -from typing import Union - -from numba import njit - - -def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None): - """Replaces all instances of the current fill value (``original_fill``) in - (``grid_var``) with (``new_fill``) and converts to the dtype defined by - (``new_dtype``) - - Parameters - ---------- - grid_var : np.ndarray - grid variable to be modified - original_fill : constant - original fill value used in (``grid_var``) - new_fill : constant - new fill value to be used in (``grid_var``) - new_dtype : np.dtype, optional - new data type to convert (``grid_var``) to - - Returns - ---------- - grid_var : xarray.Dataset - Input Dataset with correct fill value and dtype - """ - - # locations of fill values - if original_fill is not None and np.isnan(original_fill): - fill_val_idx = np.isnan(grid_var) - else: - fill_val_idx = grid_var == original_fill - - # convert to new data type - if new_dtype != grid_var.dtype and new_dtype is not None: - grid_var = grid_var.astype(new_dtype) - - # ensure fill value can be represented with current integer data type - if np.issubdtype(new_dtype, np.integer): - int_min = np.iinfo(grid_var.dtype).min - int_max = np.iinfo(grid_var.dtype).max - # ensure new_fill is in range [int_min, int_max] - if new_fill < int_min or new_fill > int_max: - raise ValueError( - f"New fill value: {new_fill} not representable by" - f" integer dtype: {grid_var.dtype}" - ) - - # ensure non-nan fill value can be represented with current float data type - elif np.issubdtype(new_dtype, np.floating) and not np.isnan(new_fill): - float_min = np.finfo(grid_var.dtype).min - float_max = np.finfo(grid_var.dtype).max - # ensure new_fill is in range [float_min, float_max] - if new_fill < float_min or new_fill > float_max: - raise ValueError( - f"New fill value: {new_fill} not representable by" - f" float dtype: {grid_var.dtype}" - ) - else: - raise ValueError( - f"Data type {grid_var.dtype} not supported" f"for grid variables" - ) - - # replace all zeros with a fill value - grid_var[fill_val_idx] = new_fill - - return grid_var def _inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old): @@ -428,168 +359,3 @@ def _get_lonlat_rad_face_edge_nodes( face_edges_lonlat_rad[valid_mask, 1] = node_lat_rad[valid_edges] return face_edges_lonlat_rad.reshape(n_face, n_max_face_edges, 2, 2) - - -def _xyz_to_lonlat_rad( - x: Union[np.ndarray, float], - y: Union[np.ndarray, float], - z: Union[np.ndarray, float], - normalize: bool = True, -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Converts Cartesian x, y, z coordinates in Spherical latitude and - longitude coordinates in degrees. - - Parameters - ---------- - x : Union[np.ndarray, float] - Cartesian x coordinates - y: Union[np.ndarray, float] - Cartesiain y coordinates - z: Union[np.ndarray, float] - Cartesian z coordinates - normalize: bool - Flag to select whether to normalize the coordinates - - Returns - ------- - lon : Union[np.ndarray, float] - Longitude in radians - lat: Union[np.ndarray, float] - Latitude in radians - """ - - if normalize: - x, y, z = _normalize_xyz(x, y, z) - denom = np.abs(x * x + y * y + z * z) - x /= denom - y /= denom - z /= denom - - lon = np.arctan2(y, x, dtype=np.float64) - lat = np.arcsin(z, dtype=np.float64) - - # set longitude range to [0, pi] - lon = np.mod(lon, 2 * np.pi) - - z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE - - lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) - lon = np.where(z_mask, 0.0, lon) - - return lon, lat - - -@njit -def _xyz_to_lonlat_rad_no_norm( - x: Union[np.ndarray, float], - y: Union[np.ndarray, float], - z: Union[np.ndarray, float], -): - """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 - """ - - lon = math.atan2(y, x) - lat = math.asin(z) - - # set longitude range to [0, pi] - lon = np.mod(lon, 2 * np.pi) - - z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE - - lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) - lon = np.where(z_mask, 0.0, lon) - - return lon, lat - - -def _normalize_xyz( - x: Union[np.ndarray, float], - y: Union[np.ndarray, float], - z: Union[np.ndarray, float], -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Normalizes a set of Cartesiain coordinates.""" - denom = np.linalg.norm( - np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2, axis=0 - ) - - x_norm = x / denom - y_norm = y / denom - z_norm = z / denom - return x_norm, y_norm, z_norm - - -@njit(cache=True) -def _lonlat_rad_to_xyz( - lon: Union[np.ndarray, float], - lat: Union[np.ndarray, float], -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Converts Spherical lon and lat coordinates into Cartesian x, y, z - coordinates.""" - x = np.cos(lon) * np.cos(lat) - y = np.sin(lon) * np.cos(lat) - z = np.sin(lat) - - return x, y, z - - -def _xyz_to_lonlat_deg( - x: Union[np.ndarray, float], - y: Union[np.ndarray, float], - z: Union[np.ndarray, float], - normalize: bool = True, -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Converts Cartesian x, y, z coordinates in Spherical latitude and - longitude coordinates in degrees. - - Parameters - ---------- - x : Union[np.ndarray, float] - Cartesian x coordinates - y: Union[np.ndarray, float] - Cartesiain y coordinates - z: Union[np.ndarray, float] - Cartesian z coordinates - normalize: bool - Flag to select whether to normalize the coordinates - - Returns - ------- - lon : Union[np.ndarray, float] - Longitude in degrees - lat: Union[np.ndarray, float] - Latitude in degrees - """ - lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z, normalize=normalize) - - lon = np.rad2deg(lon_rad) - lat = np.rad2deg(lat_rad) - - lon = (lon + 180) % 360 - 180 - return lon, lat - - -@njit -def _normalize_xyz_scalar(x: float, y: float, z: float): - denom = np.linalg.norm(np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2) - x_norm = x / denom - y_norm = y / denom - z_norm = z / denom - return x_norm, y_norm, z_norm diff --git a/uxarray/io/_exodus.py b/uxarray/io/_exodus.py index c0c5d2778..8c91913fa 100644 --- a/uxarray/io/_exodus.py +++ b/uxarray/io/_exodus.py @@ -6,7 +6,7 @@ from uxarray.grid.connectivity import _replace_fill_values from uxarray.constants import INT_DTYPE, INT_FILL_VALUE -from uxarray.grid.utils import _lonlat_rad_to_xyz, _xyz_to_lonlat_deg +from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_deg from uxarray.conventions import ugrid From cc330d47cdcb4f83de0cb84022f14c4af45debfa Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 19 Sep 2024 11:54:50 -0500 Subject: [PATCH 23/29] Update benchmarks/mpas_ocean.py Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> --- benchmarks/mpas_ocean.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index baf86c9f6..e07d867ab 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -146,7 +146,7 @@ def time_construct_hole_edge_indices(self, resolution): ux.grid.geometry._construct_hole_edge_indices(self.uxds.uxgrid.edge_face_connectivity) -class ConstructFaceLatLon: +class ConstructFaceLatLon(GridBenchmark): param_names = ['resolution'] params = ['480km', '120km'] From 1744a308d2df8c242fb2ef05a9e5cf295a909270 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 19 Sep 2024 11:54:58 -0500 Subject: [PATCH 24/29] Update benchmarks/mpas_ocean.py Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> --- benchmarks/mpas_ocean.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index e07d867ab..3b81c2226 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -148,14 +148,6 @@ def time_construct_hole_edge_indices(self, resolution): class ConstructFaceLatLon(GridBenchmark): - 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 From 57e3112d313f14f665a96a3b44c55e73c71ee2f2 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 19 Sep 2024 12:47:09 -0500 Subject: [PATCH 25/29] o Fix imports and benchmarks --- benchmarks/mpas_ocean.py | 1 - test/test_centroids.py | 2 +- uxarray/grid/area.py | 2 +- uxarray/grid/coordinates.py | 65 ------------------------------------- 4 files changed, 2 insertions(+), 68 deletions(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index b887129a2..3343c969f 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -173,4 +173,3 @@ def teardown(self, resolution): def time_check_norm(self, resolution): from uxarray.grid.validation import _check_normalization _check_normalization(self.uxgrid) - diff --git a/test/test_centroids.py b/test/test_centroids.py index eefada5c5..70c5cb208 100644 --- a/test/test_centroids.py +++ b/test/test_centroids.py @@ -5,7 +5,7 @@ import uxarray as ux from pathlib import Path 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 +from uxarray.grid.coordinates import _normalize_xyz current_path = Path(os.path.dirname(os.path.realpath(__file__))) diff --git a/uxarray/grid/area.py b/uxarray/grid/area.py index e9ca6096e..13720e786 100644 --- a/uxarray/grid/area.py +++ b/uxarray/grid/area.py @@ -1,6 +1,6 @@ import numpy as np -from uxarray.grid.utils import _lonlat_rad_to_xyz +from uxarray.grid.coordinates import _lonlat_rad_to_xyz from numba import njit, config from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 2e2824ab3..cc8b84318 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -9,7 +9,6 @@ from numba import njit from uxarray.constants import ERROR_TOLERANCE from typing import Union -import math @njit(cache=True) @@ -715,70 +714,6 @@ def _set_desired_longitude_range(ds): ds[lon_name].data = (ds[lon_name].data + 180) % 360 - 180 -def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None): - """Replaces all instances of the current fill value (``original_fill``) in - (``grid_var``) with (``new_fill``) and converts to the dtype defined by - (``new_dtype``) - - Parameters - ---------- - grid_var : np.ndarray - grid variable to be modified - original_fill : constant - original fill value used in (``grid_var``) - new_fill : constant - new fill value to be used in (``grid_var``) - new_dtype : np.dtype, optional - new data type to convert (``grid_var``) to - - Returns - ---------- - grid_var : xarray.Dataset - Input Dataset with correct fill value and dtype - """ - - # locations of fill values - if original_fill is not None and np.isnan(original_fill): - fill_val_idx = np.isnan(grid_var) - else: - fill_val_idx = grid_var == original_fill - - # convert to new data type - if new_dtype != grid_var.dtype and new_dtype is not None: - grid_var = grid_var.astype(new_dtype) - - # ensure fill value can be represented with current integer data type - if np.issubdtype(new_dtype, np.integer): - int_min = np.iinfo(grid_var.dtype).min - int_max = np.iinfo(grid_var.dtype).max - # ensure new_fill is in range [int_min, int_max] - if new_fill < int_min or new_fill > int_max: - raise ValueError( - f"New fill value: {new_fill} not representable by" - f" integer dtype: {grid_var.dtype}" - ) - - # ensure non-nan fill value can be represented with current float data type - elif np.issubdtype(new_dtype, np.floating) and not np.isnan(new_fill): - float_min = np.finfo(grid_var.dtype).min - float_max = np.finfo(grid_var.dtype).max - # ensure new_fill is in range [float_min, float_max] - if new_fill < float_min or new_fill > float_max: - raise ValueError( - f"New fill value: {new_fill} not representable by" - f" float dtype: {grid_var.dtype}" - ) - else: - raise ValueError( - f"Data type {grid_var.dtype} not supported" f"for grid variables" - ) - - # replace all zeros with a fill value - grid_var[fill_val_idx] = new_fill - - return grid_var - - def _xyz_to_lonlat_rad( x: Union[np.ndarray, float], y: Union[np.ndarray, float], From 83cdc6796605b65d2383f960d2052a37a958fdf9 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 14:45:38 -0500 Subject: [PATCH 26/29] resolve circuluar import --- uxarray/grid/arcs.py | 28 +----- uxarray/grid/coordinates.py | 196 ++++++++++++++++++------------------ uxarray/grid/utils.py | 27 +++++ 3 files changed, 128 insertions(+), 123 deletions(-) diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 2bffa2177..c2ef5f009 100644 --- a/uxarray/grid/arcs.py +++ b/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. diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index cc8b84318..415fa2fe1 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -4,12 +4,13 @@ import warnings from uxarray.conventions import ugrid -from uxarray.grid.arcs import _angle_of_2_vectors from numba import njit from uxarray.constants import ERROR_TOLERANCE from typing import Union +from uxarray.grid.utils import _angle_of_2_vectors + @njit(cache=True) def _lonlat_rad_to_xyz( @@ -361,64 +362,72 @@ def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_fa return _normalize_xyz(centroid_x, centroid_y, centroid_z) -def _populate_face_centerpoints(grid, repopulate=False): - """Calculates the face centerpoints using Welzl's algorithm. It is a - randomized algorithm for finding the center and radius of the smallest - circle that encloses a set of points. It is here adapted to work on a unit - sphere. Also, this algorithm cannot be guaranteed to work on concave - polygons. +def _welzl_recursive(points, boundary, R): + """Recursive helper function for Welzl's algorithm to find the smallest + enclosing circle. Parameters ---------- - grid : Grid - The grid containing the nodes and faces. - repopulate : bool, optional - Bool used to turn on/off repopulating the face coordinates of the centerpoints, default is False. + points : numpy.ndarray + The set of points to consider. + boundary : numpy.ndarray + The current boundary points of the minimal enclosing circle. + R : tuple + The current minimal enclosing circle. Returns ------- - None, populates the grid with the face centerpoints: face_lon, face_lat + tuple + The smallest enclosing circle as a tuple of center and radius. """ - # warnings.warn("This cannot be guaranteed to work correctly on concave polygons") - - node_lon = grid.node_lon.values - node_lat = grid.node_lat.values - - centerpoint_lat = [] - centerpoint_lon = [] + # Base case: no points or boundary has 3 points + if len(points) == 0 or len(boundary) == 3: + # Construct the minimal circle based on the number of boundary points + if len(boundary) == 0: + # Return a default circle if no boundary points are available + return ((0.0, 0.0), 0.0) + elif len(boundary) == 1: + return _circle_from_two_points(boundary[0], boundary[0]) + elif len(boundary) == 2: + return _circle_from_two_points(boundary[0], boundary[1]) + elif len(boundary) == 3: + return _circle_from_three_points(boundary[0], boundary[1], boundary[2]) - face_nodes = grid.face_node_connectivity.values - n_nodes_per_face = grid.n_nodes_per_face.values + p = points[-1] + temp_points = points[:-1] + circle = _welzl_recursive(temp_points, boundary, R) - # Check if the centerpoints are already populated - if "face_lon" not in grid._ds or repopulate: - centerpoint_lon, centerpoint_lat = _construct_face_centerpoints( - node_lon, node_lat, face_nodes, n_nodes_per_face + # Check if the chosen point is inside the current circle + if circle and _is_inside_circle(circle, p): + return circle + # If not, the point must be on the boundary of the minimal enclosing circle + else: + new_boundary = np.empty( + (boundary.shape[0] + 1, boundary.shape[1]), dtype=boundary.dtype ) - # get the cartesian coordinates of the centerpoints - ctrpt_x, ctrpt_y, ctrpt_z = _lonlat_rad_to_xyz(centerpoint_lon, centerpoint_lat) + new_boundary[:-1] = boundary + new_boundary[-1] = p + return _welzl_recursive(temp_points, new_boundary, R) - # set the grid variables for centerpoints - if "face_lon" not in grid._ds or repopulate: - grid._ds["face_lon"] = xr.DataArray( - centerpoint_lon, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LON_ATTRS - ) - grid._ds["face_lat"] = xr.DataArray( - centerpoint_lat, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LAT_ATTRS - ) - if "face_x" not in grid._ds or repopulate: - grid._ds["face_x"] = xr.DataArray( - ctrpt_x, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_X_ATTRS - ) +def _smallest_enclosing_circle(points): + """Find the smallest circle that encloses all given points on a unit sphere + using Welzl's algorithm. - grid._ds["face_y"] = xr.DataArray( - ctrpt_y, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Y_ATTRS - ) + Parameters + ---------- + points : numpy.ndarray + An array of points as tuples of (lon, lat). - grid._ds["face_z"] = xr.DataArray( - ctrpt_z, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Z_ATTRS - ) + Returns + ------- + tuple + The smallest enclosing circle as a tuple of center and radius. + """ + np.random.shuffle( + points + ) # Randomize the input to increase the chance of an optimal solution + return _welzl_recursive(points, np.empty((0, 2)), None) @njit @@ -513,73 +522,64 @@ def _is_inside_circle(circle, point): return distance <= radius -@njit -def _welzl_recursive(points, boundary, R): - """Recursive helper function for Welzl's algorithm to find the smallest - enclosing circle. +def _populate_face_centerpoints(grid, repopulate=False): + """Calculates the face centerpoints using Welzl's algorithm. It is a + randomized algorithm for finding the center and radius of the smallest + circle that encloses a set of points. It is here adapted to work on a unit + sphere. Also, this algorithm cannot be guaranteed to work on concave + polygons. Parameters ---------- - points : numpy.ndarray - The set of points to consider. - boundary : numpy.ndarray - The current boundary points of the minimal enclosing circle. - R : tuple - The current minimal enclosing circle. + grid : Grid + The grid containing the nodes and faces. + repopulate : bool, optional + Bool used to turn on/off repopulating the face coordinates of the centerpoints, default is False. Returns ------- - tuple - The smallest enclosing circle as a tuple of center and radius. + None, populates the grid with the face centerpoints: face_lon, face_lat """ - # Base case: no points or boundary has 3 points - if len(points) == 0 or len(boundary) == 3: - # Construct the minimal circle based on the number of boundary points - if len(boundary) == 0: - # Return a default circle if no boundary points are available - return ((0.0, 0.0), 0.0) - elif len(boundary) == 1: - return _circle_from_two_points(boundary[0], boundary[0]) - elif len(boundary) == 2: - return _circle_from_two_points(boundary[0], boundary[1]) - elif len(boundary) == 3: - return _circle_from_three_points(boundary[0], boundary[1], boundary[2]) + # warnings.warn("This cannot be guaranteed to work correctly on concave polygons") - p = points[-1] - temp_points = points[:-1] - circle = _welzl_recursive(temp_points, boundary, R) + node_lon = grid.node_lon.values + node_lat = grid.node_lat.values - # Check if the chosen point is inside the current circle - if circle and _is_inside_circle(circle, p): - return circle - # If not, the point must be on the boundary of the minimal enclosing circle - else: - new_boundary = np.empty( - (boundary.shape[0] + 1, boundary.shape[1]), dtype=boundary.dtype + centerpoint_lat = [] + centerpoint_lon = [] + + face_nodes = grid.face_node_connectivity.values + n_nodes_per_face = grid.n_nodes_per_face.values + + # Check if the centerpoints are already populated + if "face_lon" not in grid._ds or repopulate: + centerpoint_lon, centerpoint_lat = _construct_face_centerpoints( + node_lon, node_lat, face_nodes, n_nodes_per_face ) - new_boundary[:-1] = boundary - new_boundary[-1] = p - return _welzl_recursive(temp_points, new_boundary, R) + # get the cartesian coordinates of the centerpoints + ctrpt_x, ctrpt_y, ctrpt_z = _lonlat_rad_to_xyz(centerpoint_lon, centerpoint_lat) + # set the grid variables for centerpoints + if "face_lon" not in grid._ds or repopulate: + grid._ds["face_lon"] = xr.DataArray( + centerpoint_lon, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LON_ATTRS + ) + grid._ds["face_lat"] = xr.DataArray( + centerpoint_lat, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LAT_ATTRS + ) -def _smallest_enclosing_circle(points): - """Find the smallest circle that encloses all given points on a unit sphere - using Welzl's algorithm. + if "face_x" not in grid._ds or repopulate: + grid._ds["face_x"] = xr.DataArray( + ctrpt_x, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_X_ATTRS + ) - Parameters - ---------- - points : numpy.ndarray - An array of points as tuples of (lon, lat). + grid._ds["face_y"] = xr.DataArray( + ctrpt_y, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Y_ATTRS + ) - Returns - ------- - tuple - The smallest enclosing circle as a tuple of center and radius. - """ - np.random.shuffle( - points - ) # Randomize the input to increase the chance of an optimal solution - return _welzl_recursive(points, np.empty((0, 2)), None) + grid._ds["face_z"] = xr.DataArray( + ctrpt_z, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Z_ATTRS + ) def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_face): diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py index 9701cef0b..63cb60213 100644 --- a/uxarray/grid/utils.py +++ b/uxarray/grid/utils.py @@ -3,6 +3,33 @@ import warnings import uxarray.utils.computing as ac_utils +from numba import njit + + +@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 _inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old): """Calculate the inverse Jacobian matrix for a given set of parameters. From 795e5274a2d594d1290f3dcd7828a05da6774e4d Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 14:46:55 -0500 Subject: [PATCH 27/29] clean up quad hex benchmark --- benchmarks/quad_hexagon.py | 36 ------------------------------------ 1 file changed, 36 deletions(-) diff --git a/benchmarks/quad_hexagon.py b/benchmarks/quad_hexagon.py index aeb7559bb..4364f39b0 100644 --- a/benchmarks/quad_hexagon.py +++ b/benchmarks/quad_hexagon.py @@ -34,39 +34,3 @@ 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) From b54dfca91ac3b01fa0abd6b29e517fd5113da768 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 14:50:10 -0500 Subject: [PATCH 28/29] update internal API --- docs/internal_api/index.rst | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst index 2b42a56f5..34586a506 100644 --- a/docs/internal_api/index.rst +++ b/docs/internal_api/index.rst @@ -162,6 +162,12 @@ Coordinates 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 @@ -180,15 +186,10 @@ 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 - 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 From f1058afe74a9e2b9bf706272cb4ae3b04434e057 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 15:20:34 -0500 Subject: [PATCH 29/29] update benchmark --- benchmarks/mpas_ocean.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index 1cfe73d02..0c761434c 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -139,19 +139,10 @@ def time_construct_hole_edge_indices(self, resolution): class ConstructFaceLatLon(GridBenchmark): def time_welzl(self, resolution): - from uxarray.grid.coordinates import _construct_face_centerpoint - _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) + self.uxgrid.construct_face_centers(method='welzl') 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) + self.uxgrid.construct_face_centers(method='cartesian average') class CheckNorm: param_names = ['resolution']