From a031ea16558a6b1f983475757125f09fda20beaa Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Tue, 11 Mar 2025 19:00:31 -0700 Subject: [PATCH 1/9] initial commit --- uxarray/core/zonal.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py index e57975b87..cdfd1f69c 100644 --- a/uxarray/core/zonal.py +++ b/uxarray/core/zonal.py @@ -6,6 +6,8 @@ from uxarray.grid.utils import _get_cartesian_face_edge_nodes +#TODO: Optimize the entire zonal average function to avoid unnecessary computations + def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=False): """Computes the non-conservative zonal mean across one or more latitudes.""" uxgrid = uxda.uxgrid From 3fc6358dca0335d8dc3cf98df64c4a288ec27e98 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Tue, 11 Mar 2025 20:15:42 -0700 Subject: [PATCH 2/9] Finish all easyckean without touching the latlon bounds --- uxarray/core/zonal.py | 2 +- uxarray/grid/integrate.py | 36 ++++++++++++++++++++---------------- 2 files changed, 21 insertions(+), 17 deletions(-) diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py index cdfd1f69c..80a9f6c55 100644 --- a/uxarray/core/zonal.py +++ b/uxarray/core/zonal.py @@ -6,7 +6,6 @@ from uxarray.grid.utils import _get_cartesian_face_edge_nodes -#TODO: Optimize the entire zonal average function to avoid unnecessary computations def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=False): """Computes the non-conservative zonal mean across one or more latitudes.""" @@ -60,3 +59,4 @@ def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=Fal ) return result + diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index d0beade57..c1ebcb9f9 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -215,19 +215,20 @@ def _get_zonal_face_interval( face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1] try: - unique_intersections, pt_lon_min, pt_lon_max = ( + + unique_intersections_cart,unique_intersection_lonlat, pt_lon_min, pt_lon_max = ( _get_faces_constLat_intersection_info( face_edges_cart, latitude_cart, is_GCA_list, is_latlonface ) ) # If there's exactly one intersection, the face is only "touched" - if len(unique_intersections) == 1: + if len(unique_intersections_cart) == 1: return pl.DataFrame({"start": [0.0], "end": [0.0]}) - # Convert intersection points to (lon, lat) in radians + # Extract unique longitudes from intersection points longitudes = np.array( - [_xyz_to_lonlat_rad(*pt.tolist())[0] for pt in unique_intersections] + [pt[0] for pt in unique_intersection_lonlat] ) # Handle special wrap-around cases (crossing anti-meridian, etc.) @@ -378,7 +379,8 @@ def _get_faces_constLat_intersection_info( ------- tuple A tuple containing: - - intersections_pts_list_cart (list): A list of intersection points where each point is where an edge intersects with the latitude. + - intersections_pts_list_cart (list): A list of intersection points in cartesian where each point is where an edge intersects with the latitude. + - intersections_pts_list_lon (list): A list of intersection points in lonlat where each point is where an edge intersects with the latitude. - pt_lon_min (float): The min longnitude of the interseted intercal in radian if any; otherwise, None.. - pt_lon_max (float): The max longnitude of the interseted intercal in radian, if any; otherwise, None. """ @@ -418,28 +420,30 @@ def _get_faces_constLat_intersection_info( intersections_pts_list_cart.append(intersections[0]) # Find the unique intersection points - unique_intersections = np.unique(intersections_pts_list_cart, axis=0) + unique_intersections_cart = np.unique(intersections_pts_list_cart, axis=0) - if len(unique_intersections) == 2: + if len(unique_intersections_cart) == 2: unique_intersection_lonlat = np.array( - [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections] + [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections_cart] ) sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0) pt_lon_min, pt_lon_max = sorted_lonlat[:, 0] - return unique_intersections, pt_lon_min, pt_lon_max - elif len(unique_intersections) == 1: - return unique_intersections, None, None - elif len(unique_intersections) != 0 and len(unique_intersections) != 1: + return unique_intersections_cart,unique_intersection_lonlat, pt_lon_min, pt_lon_max + elif len(unique_intersections_cart) == 1: + return unique_intersections_cart, np.array(_xyz_to_lonlat_rad(unique_intersections_cart[0][0], + unique_intersections_cart[0][1], + unique_intersections_cart[0][2])), None, None + elif len(unique_intersections_cart) != 0 and len(unique_intersections_cart) != 1: # If the unique intersections numbers is larger than n_edges * 2, then it means the face is concave - if len(unique_intersections) > len(valid_edges) * 2: + if len(unique_intersections_cart) > len(valid_edges) * 2: raise ValueError( "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly" ) else: # Now return all the intersections points and the pt_lon_min, pt_lon_max unique_intersection_lonlat = np.array( - [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections] + [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections_cart] ) sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0) @@ -449,8 +453,8 @@ def _get_faces_constLat_intersection_info( np.max(sorted_lonlat[:, 0]), ) - return unique_intersections, pt_lon_min, pt_lon_max - elif len(unique_intersections) == 0: + return unique_intersections_cart, unique_intersection_lonlat, pt_lon_min, pt_lon_max + elif len(unique_intersections_cart) == 0: raise ValueError( "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" ) From d2a3aa70ce0243d7b06726b6058b5e0a2a332080 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Wed, 12 Mar 2025 23:36:16 -0700 Subject: [PATCH 3/9] Try adding `face_edges_spherical` and `face_edges_cartesian` as properties to grid --- uxarray/core/zonal.py | 25 ++++-- uxarray/grid/geometry.py | 168 ++++++++++++++++++++++----------------- uxarray/grid/grid.py | 16 ++++ 3 files changed, 132 insertions(+), 77 deletions(-) diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py index 80a9f6c55..a9fcdc358 100644 --- a/uxarray/core/zonal.py +++ b/uxarray/core/zonal.py @@ -3,7 +3,10 @@ from uxarray.grid.integrate import _zonal_face_weights, _zonal_face_weights_robust -from uxarray.grid.utils import _get_cartesian_face_edge_nodes +from uxarray.grid.utils import ( + _get_cartesian_face_edge_nodes, + _get_lonlat_rad_face_edge_nodes, +) @@ -19,23 +22,35 @@ def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=Fal # Create a NumPy array for storing results result = np.zeros(shape, dtype=uxda.dtype) - faces_edge_nodes_xyz = _get_cartesian_face_edge_nodes( + + faces_edges_cartesian = _get_cartesian_face_edge_nodes( uxgrid.face_node_connectivity.values, uxgrid.n_face, - uxgrid.n_max_face_nodes, + uxgrid.n_max_face_edges, uxgrid.node_x.values, uxgrid.node_y.values, uxgrid.node_z.values, ) - bounds = uxgrid.bounds.values + + faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( + uxgrid.face_node_connectivity.values, + uxgrid.n_face, + uxgrid.n_max_face_edges, + uxgrid.node_lon.values, + uxgrid.node_lat.values, + ) + + # Check if the bounds are available + if "bounds" not in uxgrid._ds: + uxgrid._populate_bounds() for i, lat in enumerate(latitudes): face_indices = uxda.uxgrid.get_faces_at_constant_latitude(lat) z = np.sin(np.deg2rad(lat)) - faces_edge_nodes_xyz_candidate = faces_edge_nodes_xyz[face_indices, :, :, :] + faces_edge_nodes_xyz_candidate = faces_edges_cartesian[face_indices, :, :, :] n_nodes_per_face_candidate = n_nodes_per_face[face_indices] diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 1c5660da6..8870bf573 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -103,7 +103,7 @@ def _unique_points(points, tolerance=ERROR_TOLERANCE): @njit(cache=True) def _pad_closed_face_nodes( - face_node_connectivity, n_face, n_max_face_nodes, n_nodes_per_face + face_node_connectivity, n_face, n_max_face_nodes, n_nodes_per_face ): """Pads a closed array of face nodes by inserting the first element at any point a fill value is encountered. @@ -123,14 +123,14 @@ def _pad_closed_face_nodes( def _build_polygon_shells( - node_lon, - node_lat, - face_node_connectivity, - n_face, - n_max_face_nodes, - n_nodes_per_face, - projection=None, - central_longitude=0.0, + node_lon, + node_lat, + face_node_connectivity, + n_face, + n_max_face_nodes, + n_nodes_per_face, + projection=None, + central_longitude=0.0, ): """Builds an array of polygon shells, which can be used with Shapely to construct polygons.""" @@ -266,7 +266,7 @@ def _grid_to_polygon_geodataframe(grid, periodic_elements, projection, project, def _build_geodataframe_without_antimeridian( - polygon_shells, projected_polygon_shells, antimeridian_face_indices, engine + polygon_shells, projected_polygon_shells, antimeridian_face_indices, engine ): """Builds a ``spatialpandas.GeoDataFrame`` or ``geopandas.GeoDataFrame``excluding any faces that cross the @@ -295,10 +295,10 @@ def _build_geodataframe_without_antimeridian( def _build_geodataframe_with_antimeridian( - polygon_shells, - projected_polygon_shells, - antimeridian_face_indices, - engine, + polygon_shells, + projected_polygon_shells, + antimeridian_face_indices, + engine, ): """Builds a ``spatialpandas.GeoDataFrame`` or ``geopandas.GeoDataFrame`` including any faces that cross the antimeridian.""" @@ -317,9 +317,9 @@ def _build_geodataframe_with_antimeridian( def _build_corrected_shapely_polygons( - polygon_shells, - projected_polygon_shells, - antimeridian_face_indices, + polygon_shells, + projected_polygon_shells, + antimeridian_face_indices, ): if projected_polygon_shells is not None: # use projected shells if a projection is applied @@ -438,7 +438,7 @@ def _build_corrected_polygon_shells(polygon_shells): def _grid_to_matplotlib_polycollection( - grid, periodic_elements, projection=None, **kwargs + grid, periodic_elements, projection=None, **kwargs ): """Constructs and returns a ``matplotlib.collections.PolyCollection``""" @@ -642,7 +642,7 @@ def _get_polygons(grid, periodic_elements, projection=None, apply_projection=Tru def _grid_to_matplotlib_linecollection( - grid, periodic_elements, projection=None, **kwargs + grid, periodic_elements, projection=None, **kwargs ): """Constructs and returns a ``matplotlib.collections.LineCollection``""" @@ -913,7 +913,7 @@ def _check_intersection(ref_edge_xyz, edges_xyz): for i in range(n_edges): edge_xyz = edges_xyz[i] if allclose( - intersection_point, edge_xyz[0], atol=ERROR_TOLERANCE + intersection_point, edge_xyz[0], atol=ERROR_TOLERANCE ) or allclose(intersection_point, edge_xyz[1], atol=ERROR_TOLERANCE): return 0 @@ -1140,10 +1140,10 @@ def insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True): @njit(cache=True) def _populate_face_latlon_bound( - face_edges_xyz, - face_edges_lonlat, - is_latlonface=False, - is_GCA_list=None, + face_edges_xyz, + face_edges_lonlat, + is_latlonface=False, + is_GCA_list=None, ): # Check if face_edges contains pole points has_north_pole = pole_point_inside_polygon(1, face_edges_xyz, face_edges_lonlat) @@ -1198,9 +1198,9 @@ def _populate_face_latlon_bound( # Check if the node matches the pole point or if the pole point is within the edge max_abs_diff = np.max(np.abs(n1_cart - pole_point_xyz)) if max_abs_diff <= ERROR_TOLERANCE or point_within_gca( - pole_point_xyz, - n1_cart, - n2_cart, + pole_point_xyz, + n1_cart, + n2_cart, ): is_center_pole = False face_latlon_array = insert_pt_in_latlonbox( @@ -1292,15 +1292,15 @@ def _populate_face_latlon_bound( # Insert extreme latitude points into the latlonbox if ( - abs(node1_lat_rad - lat_max) > ERROR_TOLERANCE - and abs(node2_lat_rad - lat_max) > ERROR_TOLERANCE + abs(node1_lat_rad - lat_max) > ERROR_TOLERANCE + and abs(node2_lat_rad - lat_max) > ERROR_TOLERANCE ): face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([lat_max, node1_lon_rad]) ) elif ( - abs(node1_lat_rad - lat_min) > ERROR_TOLERANCE - and abs(node2_lat_rad - lat_min) > ERROR_TOLERANCE + abs(node1_lat_rad - lat_min) > ERROR_TOLERANCE + and abs(node2_lat_rad - lat_min) > ERROR_TOLERANCE ): face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([lat_min, node1_lon_rad]) @@ -1315,23 +1315,23 @@ def _populate_face_latlon_bound( @njit(cache=True, parallel=True) def compute_temp_latlon_array( - face_node_connectivity, - faces_edges_cartesian, - faces_edges_lonlat_rad, - n_nodes_per_face, - is_latlonface, - is_face_GCA_list, - INT_FILL_VALUE, + face_node_connectivity, + faces_edges_cartesian, + faces_edges_lonlat_rad, + n_nodes_per_face, + is_latlonface, + is_face_GCA_list, + INT_FILL_VALUE, ): n_face = face_node_connectivity.shape[0] temp_latlon_array = np.full((n_face, 2, 2), INT_FILL_VALUE, dtype=np.float64) for face_idx in prange(n_face): cur_face_edges_cartesian = faces_edges_cartesian[ - face_idx, 0 : n_nodes_per_face[face_idx] - ] + face_idx, 0: n_nodes_per_face[face_idx] + ] cur_faces_edges_lonlat_rad = faces_edges_lonlat_rad[ - face_idx, 0 : n_nodes_per_face[face_idx] - ] + face_idx, 0: n_nodes_per_face[face_idx] + ] if is_face_GCA_list is not None: is_GCA_list = is_face_GCA_list[face_idx] else: @@ -1347,7 +1347,12 @@ def compute_temp_latlon_array( def _populate_bounds( - grid, is_latlonface: bool = False, is_face_GCA_list=None, return_array=False + grid, + faces_edges_cartesian: np.ndarray = None, + faces_edges_lonlat_rad: np.ndarray = None, + is_latlonface: bool = False, + is_face_GCA_list=None, + return_array=False ): """Populates the bounds of the grid based on the geometry of its faces, taking into account special conditions such as faces crossing the @@ -1357,6 +1362,21 @@ def _populate_bounds( Parameters ---------- + grid : object + The grid object whose internal representation will be updated with the computed + face bounds. + + faces_edges_cartesian : np.ndarray, optional + An array of shape (n_face, n_max_face_edges, 2, 3) containing the Cartesian coordinates + of each face's edges. This array may include dummy values for grids with holes. + Default is None. + + faces_edges_lonlat_rad : np.ndarray, optional + An array of shape (n_face, n_max_face_edges, 2, 2) containing the longitude and latitude + (in radians) coordinates of each face's edges. This array may include dummy values for + grids with holes. + Default is None. + is_latlonface : bool, optional A global flag that indicates if faces are latlon faces. If True, all faces are treated as latlon faces, meaning that all edges are either longitude or @@ -1407,25 +1427,29 @@ def _populate_bounds( """ # Ensure grid's cartesian coordinates are normalized + # TODO: Is it possible to have a more flexible normalization? (Avoid duplicated normalizations) grid.normalize_cartesian_coordinates() # Prepare data for Numba functions - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + # TODO: remove the duplicate call for this function in both zonal average and bounds + if faces_edges_cartesian is None: + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) - faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_lon.values, - grid.node_lat.values, - ) + if faces_edges_lonlat_rad is None: + faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_lon.values, + grid.node_lat.values, + ) n_nodes_per_face = grid.n_nodes_per_face.values @@ -1512,16 +1536,16 @@ def stereographic_projection(lon, lat, central_lon, central_lat): # Calculate constant used for calculation k = 2.0 / ( - 1.0 - + np.sin(central_lat) * np.sin(lat) - + np.cos(central_lat) * np.cos(lat) * np.cos(lon - central_lon) + 1.0 + + np.sin(central_lat) * np.sin(lat) + + np.cos(central_lat) * np.cos(lat) * np.cos(lon - central_lon) ) # Calculate the x and y coordinates x = k * np.cos(lat) * np.sin(lon - central_lon) y = k * ( - np.cos(central_lat) * np.sin(lat) - - np.sin(central_lat) * np.cos(lat) * np.cos(lon - central_lon) + np.cos(central_lat) * np.sin(lat) + - np.sin(central_lat) * np.cos(lat) * np.cos(lon - central_lon) ) return x, y @@ -1557,7 +1581,7 @@ def inverse_stereographic_projection(x, y, central_lon, central_lat): central_lat = np.deg2rad(central_lat) # Calculate constants used for calculation - p = np.sqrt(x**2 + y**2) + p = np.sqrt(x ** 2 + y ** 2) c = 2 * np.arctan(p / 2) @@ -1576,9 +1600,9 @@ def inverse_stereographic_projection(x, y, central_lon, central_lat): @njit(cache=True) def point_in_face( - edges_xyz, - point_xyz, - inclusive=True, + edges_xyz, + point_xyz, + inclusive=True, ): """Determines if a point lies inside a face. @@ -1628,9 +1652,9 @@ def point_in_face( for ind in range(len(edges_xyz)): # If the point lies on an edge, return True if inclusive if point_within_gca( - point_xyz, - edges_xyz[ind][0], - edges_xyz[ind][1], + point_xyz, + edges_xyz[ind][0], + edges_xyz[ind][1], ): if inclusive: return True @@ -1723,7 +1747,7 @@ def _populate_max_face_radius(self): @njit(cache=True) def calculate_max_face_radius( - face_node_connectivity, node_lats_rad, node_lons_rad, face_lats_rad, face_lons_rad + face_node_connectivity, node_lats_rad, node_lons_rad, face_lats_rad, face_lons_rad ): """Finds the max face radius in the mesh. Parameters @@ -1797,7 +1821,7 @@ def haversine_distance(lon_a, lat_a, lon_b, lat_b): # Haversine formula equation_in_sqrt = (np.sin(dlat / 2) ** 2) + np.cos(lat_a) * np.cos(lat_b) * ( - np.sin(dlon / 2) ** 2 + np.sin(dlon / 2) ** 2 ) distance = 2 * np.arcsin(np.sqrt(equation_in_sqrt)) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 9cae36bac..baaa03469 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -1443,6 +1443,22 @@ def face_areas(self, value): assert isinstance(value, xr.DataArray) self._ds["face_areas"] = value + @property + def face_edges_cartesian(self): + """Cartesian Coordinates for each Face Edge. + + Dimensions ``(n_face, n_max_face_edges, 2, 3)`` + """ + pass + + @property + def face_edges_spherical(self): + """Latitude Longitude Bounds for each Face in radians. + + Dimensions ``(n_face, n_max_face_edges, 2, 2)`` + """ + pass + @property def bounds(self): """Latitude Longitude Bounds for each Face in radians. From bbf52c34d184136c6481b356bbb902f6bb9d76e8 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Thu, 13 Mar 2025 17:22:57 -0700 Subject: [PATCH 4/9] Apply the initial changes --- uxarray/core/zonal.py | 28 +--------- uxarray/grid/geometry.py | 113 ++++++++++++++++++++++++++++++--------- uxarray/grid/grid.py | 20 +++++-- 3 files changed, 105 insertions(+), 56 deletions(-) diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py index a9fcdc358..eb0c6f312 100644 --- a/uxarray/core/zonal.py +++ b/uxarray/core/zonal.py @@ -3,10 +3,6 @@ from uxarray.grid.integrate import _zonal_face_weights, _zonal_face_weights_robust -from uxarray.grid.utils import ( - _get_cartesian_face_edge_nodes, - _get_lonlat_rad_face_edge_nodes, -) @@ -22,28 +18,8 @@ def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=Fal # Create a NumPy array for storing results result = np.zeros(shape, dtype=uxda.dtype) - - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - uxgrid.face_node_connectivity.values, - uxgrid.n_face, - uxgrid.n_max_face_edges, - uxgrid.node_x.values, - uxgrid.node_y.values, - uxgrid.node_z.values, - ) - - - faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( - uxgrid.face_node_connectivity.values, - uxgrid.n_face, - uxgrid.n_max_face_edges, - uxgrid.node_lon.values, - uxgrid.node_lat.values, - ) - - # Check if the bounds are available - if "bounds" not in uxgrid._ds: - uxgrid._populate_bounds() + bounds = uxgrid.bounds.values + faces_edges_cartesian = uxgrid.face_edges_cartesian.values for i, lat in enumerate(latitudes): face_indices = uxda.uxgrid.get_faces_at_constant_latitude(lat) diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 8870bf573..77b1cff4b 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -1346,13 +1346,91 @@ def compute_temp_latlon_array( return temp_latlon_array +def _populate_face_edges_cartesian(grid, return_array=False): + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + faces_edges_cartesian_xarray = xr.DataArray( + faces_edges_cartesian, + dims=["n_face", "n_max_face_edges", "two", "three"], + attrs={ + "cf_role": "face_edges_cartesian", + "_FillValue": INT_FILL_VALUE, + "long_name": "Provide the Cartesian coordinates of the edge for each face", + "start_index": INT_DTYPE(0) + }, + ) + + if return_array: + return faces_edges_cartesian_xarray + else: + grid._ds["face_edges_cartesian"] = faces_edges_cartesian_xarray + + +def _populate_faces_edges_cartesian(grid, return_array=False): + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + faces_edges_cartesian_xarray = xr.DataArray( + faces_edges_cartesian, + dims=["n_face", "n_max_face_edges", "two", "three"], + attrs={ + "cf_role": "faces_edges_cartesian", + "_FillValue": INT_FILL_VALUE, + "long_name": "Provide the Cartesian coordinates of the edge for each face", + "start_index": INT_DTYPE(0) + }, + ) + + if return_array: + return faces_edges_cartesian_xarray + else: + grid._ds["faces_edges_cartesian"] = faces_edges_cartesian_xarray + + +def _populate_faces_edges_spherical(grid, return_array=False): + faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_lon.values, + grid.node_lat.values, + ) + + faces_edges_lonlat_rad_xarray = xr.DataArray( + faces_edges_lonlat_rad, + dims=["n_face", "n_max_face_edges", "two", "two"], + attrs={ + "cf_role": "faces_edges_spherical", + "_FillValue": INT_FILL_VALUE, + "long_name": "Provide the Spherical coordinates (lon lat) in radians of the edge for each face", + "start_index": INT_DTYPE(0) + }, + ) + + if return_array: + return faces_edges_lonlat_rad_xarray + else: + grid._ds["faces_edges_spherical"] = faces_edges_lonlat_rad_xarray + + def _populate_bounds( - grid, - faces_edges_cartesian: np.ndarray = None, - faces_edges_lonlat_rad: np.ndarray = None, - is_latlonface: bool = False, - is_face_GCA_list=None, - return_array=False + grid, + is_latlonface: bool = False, + is_face_GCA_list=None, + return_array=False ): """Populates the bounds of the grid based on the geometry of its faces, taking into account special conditions such as faces crossing the @@ -1430,26 +1508,11 @@ def _populate_bounds( # TODO: Is it possible to have a more flexible normalization? (Avoid duplicated normalizations) grid.normalize_cartesian_coordinates() - # Prepare data for Numba functions - # TODO: remove the duplicate call for this function in both zonal average and bounds - if faces_edges_cartesian is None: - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + # If the face_edges_cartesian and face_edges_lonlat_rad didn't exist, + # The following call will automatically populate them + faces_edges_cartesian = grid.faces_edges_cartesian.values - if faces_edges_lonlat_rad is None: - faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_lon.values, - grid.node_lat.values, - ) + faces_edges_lonlat_rad = grid.faces_edges_spherical.values n_nodes_per_face = grid.n_nodes_per_face.values diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index baaa03469..44c72028e 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -71,6 +71,8 @@ _grid_to_matplotlib_polycollection, _grid_to_matplotlib_linecollection, _populate_bounds, + _populate_faces_edges_cartesian, + _populate_faces_edges_spherical, _construct_boundary_edge_indices, compute_temp_latlon_array, _find_faces, @@ -1444,20 +1446,24 @@ def face_areas(self, value): self._ds["face_areas"] = value @property - def face_edges_cartesian(self): + def faces_edges_cartesian(self): """Cartesian Coordinates for each Face Edge. Dimensions ``(n_face, n_max_face_edges, 2, 3)`` """ - pass + if "faces_edges_cartesian" not in self._ds: + _populate_faces_edges_cartesian(self) + return self._ds["faces_edges_cartesian"] @property - def face_edges_spherical(self): - """Latitude Longitude Bounds for each Face in radians. + def faces_edges_spherical(self): + """Latitude Longitude Coordinates for each Face in radians. Dimensions ``(n_face, n_max_face_edges, 2, 2)`` """ - pass + if "faces_edges_spherical" not in self._ds: + _populate_faces_edges_spherical(self) + return self._ds["faces_edges_spherical"] @property def bounds(self): @@ -1472,6 +1478,10 @@ def bounds(self): "This initial execution will be significantly longer.", RuntimeWarning, ) + if "faces_edges_cartesian" not in self._ds: + _populate_faces_edges_cartesian(self) + if "faces_edges_spherical" not in self._ds: + _populate_faces_edges_spherical(self) _populate_bounds(self) return self._ds["bounds"] From f994a84fc2e5e33eadedac3dfcdebe423e24b10d Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Thu, 13 Mar 2025 17:40:38 -0700 Subject: [PATCH 5/9] remove the duplicate call for `face_edges_cartesian` calculation --- test/test_geometry.py | 90 +++++++++++--------------------------- test/test_helpers.py | 32 ++++++-------- test/test_integrate.py | 13 ++---- test/test_intersections.py | 12 ++--- uxarray/core/zonal.py | 2 +- uxarray/grid/geometry.py | 52 ++++------------------ uxarray/grid/grid.py | 14 +++--- uxarray/grid/utils.py | 8 ++-- 8 files changed, 63 insertions(+), 160 deletions(-) diff --git a/test/test_geometry.py b/test/test_geometry.py index c066f612a..0c729374e 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -12,7 +12,7 @@ from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad, \ _xyz_to_lonlat_deg, _xyz_to_lonlat_rad_scalar from uxarray.grid.arcs import extreme_gca_latitude, extreme_gca_z -from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes +from uxarray.grid.utils import _get_cartesian_faces_edge_nodes, _get_lonlat_rad_faces_edge_nodes from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, _pole_point_inside_polygon_cartesian, \ stereographic_projection, inverse_stereographic_projection, point_in_face, haversine_distance @@ -1184,14 +1184,9 @@ def test_point_inside(): grid = ux.open_grid(grid_mpas_2) # Get the face edges of all faces in the grid - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) # Loop through each face for i in range(grid.n_face): @@ -1209,14 +1204,9 @@ def test_point_outside(): grid = ux.open_grid(grid_mpas_2) # Get the face edges of all faces in the grid - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) # Set the point as the face center of a different face than the face tested point_xyz = np.array([grid.face_x[1].values, grid.face_y[1].values, grid.face_z[1].values]) @@ -1232,14 +1222,9 @@ def test_point_on_node(): grid = ux.open_grid(grid_mpas_2) # Get the face edges of all faces in the grid - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) # Set the point as a node point_xyz = np.array([*faces_edges_cartesian[0][0][0]]) @@ -1263,14 +1248,9 @@ def test_point_inside_close(): # Create the grid and face edges grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) # Use point in face to determine if the point is inside or out of the face assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=False) @@ -1288,14 +1268,9 @@ def test_point_outside_close(): # Create the grid and face edges grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) # Use point in face to determine if the point is inside or out of the face assert not point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=False) @@ -1312,14 +1287,9 @@ def test_face_at_pole(): # Create the grid and face edges grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True) @@ -1334,14 +1304,9 @@ def test_face_at_antimeridian(): # Create the grid and face edges grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True) @@ -1357,14 +1322,9 @@ def test_face_normal_face(): # Create the grid and face edges grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True) diff --git a/test/test_helpers.py b/test/test_helpers.py index 76e352487..e8a4adbcd 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -11,7 +11,7 @@ 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.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.utils import _get_cartesian_faces_edge_nodes, _get_lonlat_rad_faces_edge_nodes from uxarray.grid.geometry import pole_point_inside_polygon, _pole_point_inside_polygon_cartesian try: @@ -273,9 +273,8 @@ def test_get_cartesian_face_edge_nodes_pipeline(): node_y = grid.node_y.values node_z = grid.node_z.values - face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes( - face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z - ) + face_edges_connectivity_cartesian = _get_cartesian_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges, + node_x, node_y, node_z) result = _pole_point_inside_polygon_cartesian( 'North', face_edges_connectivity_cartesian[0] @@ -298,9 +297,8 @@ def test_get_cartesian_face_edge_nodes_filled_value(): node_y = grid.node_y.values node_z = grid.node_z.values - face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes( - face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z - ) + face_edges_connectivity_cartesian = _get_cartesian_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges, + node_x, node_y, node_z) result = _pole_point_inside_polygon_cartesian( 'North', face_edges_connectivity_cartesian[0] @@ -335,9 +333,8 @@ def test_get_cartesian_face_edge_nodes_filled_value2(): node_y = np.array([v0_cart[1],v1_cart[1],v2_cart[1],v3_cart[1],v4_cart[1]]) node_z = np.array([v0_cart[2],v1_cart[2],v2_cart[2],v3_cart[2],v4_cart[2]]) - face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes( - face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z - ) + face_edges_connectivity_cartesian = _get_cartesian_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges, + node_x, node_y, node_z) correct_result = np.array([ [ @@ -368,9 +365,8 @@ def test_get_lonlat_face_edge_nodes_pipeline(): node_lon = grid.node_lon.values node_lat = grid.node_lat.values - face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes( - face_node_conn, n_face, n_max_face_edges, node_lon, node_lat - ) + face_edges_connectivity_lonlat = _get_lonlat_rad_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges, + node_lon, node_lat) face_edges_connectivity_lonlat = face_edges_connectivity_lonlat[0] face_edges_connectivity_cartesian = [] @@ -398,9 +394,8 @@ def test_get_lonlat_face_edge_nodes_filled_value(): node_lon = grid.node_lon.values node_lat = grid.node_lat.values - face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes( - face_node_conn, n_face, n_max_face_edges, node_lon, node_lat - ) + face_edges_connectivity_lonlat = _get_lonlat_rad_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges, + node_lon, node_lat) face_edges_connectivity_lonlat = face_edges_connectivity_lonlat[0] face_edges_connectivity_cartesian = [] @@ -434,9 +429,8 @@ def test_get_lonlat_face_edge_nodes_filled_value2(): node_lon = np.array([v0_rad[0],v1_rad[0],v2_rad[0],v3_rad[0],v4_rad[0]]) node_lat = np.array([v0_rad[1],v1_rad[1],v2_rad[1],v3_rad[1],v4_rad[1]]) - face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes( - face_node_conn, n_face, n_max_face_edges, node_lon, node_lat - ) + face_edges_connectivity_lonlat = _get_lonlat_rad_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges, + node_lon, node_lat) correct_result = np.array([ [ diff --git a/test/test_integrate.py b/test/test_integrate.py index ad9881c04..18289e0bf 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -19,7 +19,7 @@ _get_faces_constLat_intersection_info, _zonal_face_weights, \ _zonal_face_weights_robust -from uxarray.grid.utils import _get_cartesian_face_edge_nodes +from uxarray.grid.utils import _get_cartesian_faces_edge_nodes current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -1037,14 +1037,9 @@ def test_compare_zonal_weights(): for gridfile in gridfiles: uxgrid = ux.open_grid(gridfile) n_nodes_per_face = uxgrid.n_nodes_per_face.values - face_edge_nodes_xyz = _get_cartesian_face_edge_nodes( - uxgrid.face_node_connectivity.values, - uxgrid.n_face, - uxgrid.n_max_face_edges, - uxgrid.node_x.values, - uxgrid.node_y.values, - uxgrid.node_z.values, - ) + face_edge_nodes_xyz = _get_cartesian_faces_edge_nodes(uxgrid.face_node_connectivity.values, uxgrid.n_face, + uxgrid.n_max_face_edges, uxgrid.node_x.values, + uxgrid.node_y.values, uxgrid.node_z.values) bounds = uxgrid.bounds.values for i, lat in enumerate(latitudes): diff --git a/test/test_intersections.py b/test/test_intersections.py index f2fbfe1ee..6d52e4d46 100644 --- a/test/test_intersections.py +++ b/test/test_intersections.py @@ -161,20 +161,16 @@ def test_GCA_GCA_north_pole_angled(): def test_GCA_edge_intersection_count(): - from uxarray.grid.utils import _get_cartesian_face_edge_nodes + from uxarray.grid.utils import _get_cartesian_faces_edge_nodes # Generate a normal face that is not crossing the antimeridian or the poles vertices_lonlat = [[29.5, 11.0], [29.5, 10.0], [30.5, 10.0], [30.5, 11.0]] vertices_lonlat = np.array(vertices_lonlat) grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) - face_edge_nodes_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values) + face_edge_nodes_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) face_center_xyz = np.array([grid.face_x.values[0], grid.face_y.values[0], grid.face_z.values[0]], dtype=np.float64) north_pole_xyz = np.array([0.0, 0.0, 1.0], dtype=np.float64) diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py index eb0c6f312..02e63e787 100644 --- a/uxarray/core/zonal.py +++ b/uxarray/core/zonal.py @@ -19,7 +19,7 @@ def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=Fal result = np.zeros(shape, dtype=uxda.dtype) bounds = uxgrid.bounds.values - faces_edges_cartesian = uxgrid.face_edges_cartesian.values + faces_edges_cartesian = uxgrid.faces_edges_cartesian.values for i, lat in enumerate(latitudes): face_indices = uxda.uxgrid.get_faces_at_constant_latitude(lat) diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 77b1cff4b..13027b2d3 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -29,8 +29,8 @@ gca_gca_intersection, ) from uxarray.grid.utils import ( - _get_cartesian_face_edge_nodes, - _get_lonlat_rad_face_edge_nodes, + _get_cartesian_faces_edge_nodes, + _get_lonlat_rad_faces_edge_nodes, ) from uxarray.utils.computing import allclose, isclose @@ -1346,42 +1346,10 @@ def compute_temp_latlon_array( return temp_latlon_array -def _populate_face_edges_cartesian(grid, return_array=False): - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) - - faces_edges_cartesian_xarray = xr.DataArray( - faces_edges_cartesian, - dims=["n_face", "n_max_face_edges", "two", "three"], - attrs={ - "cf_role": "face_edges_cartesian", - "_FillValue": INT_FILL_VALUE, - "long_name": "Provide the Cartesian coordinates of the edge for each face", - "start_index": INT_DTYPE(0) - }, - ) - - if return_array: - return faces_edges_cartesian_xarray - else: - grid._ds["face_edges_cartesian"] = faces_edges_cartesian_xarray - - def _populate_faces_edges_cartesian(grid, return_array=False): - faces_edges_cartesian = _get_cartesian_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_x.values, - grid.node_y.values, - grid.node_z.values, - ) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_x.values, + grid.node_y.values, grid.node_z.values) faces_edges_cartesian_xarray = xr.DataArray( faces_edges_cartesian, @@ -1401,13 +1369,9 @@ def _populate_faces_edges_cartesian(grid, return_array=False): def _populate_faces_edges_spherical(grid, return_array=False): - faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( - grid.face_node_connectivity.values, - grid.n_face, - grid.n_max_face_edges, - grid.node_lon.values, - grid.node_lat.values, - ) + faces_edges_lonlat_rad = _get_lonlat_rad_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, + grid.n_max_face_edges, grid.node_lon.values, + grid.node_lat.values) faces_edges_lonlat_rad_xarray = xr.DataArray( faces_edges_lonlat_rad, diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 44c72028e..5596bfe15 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -14,7 +14,7 @@ Tuple, ) -from uxarray.grid.utils import _get_cartesian_face_edge_nodes +from uxarray.grid.utils import _get_cartesian_faces_edge_nodes # reader and writer imports from uxarray.io._exodus import _read_exodus, _encode_exodus @@ -2664,14 +2664,10 @@ def get_faces_containing_point( return np.empty(0, dtype=np.int64) # Get the faces in terms of their edges - face_edge_nodes_xyz = _get_cartesian_face_edge_nodes( - subset.face_node_connectivity.values, - subset.n_face, - subset.n_max_face_nodes, - subset.node_x.values, - subset.node_y.values, - subset.node_z.values, - ) + # Since this is a new subset, the cartesian_faces_edge_nodes need to be recalculated anyway + face_edge_nodes_xyz = _get_cartesian_faces_edge_nodes(subset.face_node_connectivity.values, subset.n_face, + subset.n_max_face_nodes, subset.node_x.values, + subset.node_y.values, subset.node_z.values) # Get the original face indices from the subset inverse_indices = subset.inverse_indices.face.values diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py index d9a6621c9..068479dad 100644 --- a/uxarray/grid/utils.py +++ b/uxarray/grid/utils.py @@ -138,7 +138,7 @@ def _swap_first_fill_value_with_last(arr): return arr -def _get_cartesian_face_edge_nodes( +def _get_cartesian_faces_edge_nodes( face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z ): """Construct an array to hold the edge Cartesian coordinates connectivity @@ -179,9 +179,7 @@ def _get_cartesian_face_edge_nodes( >>> node_x = np.array([0, 1, 1, 0, 1, 0]) >>> node_y = np.array([0, 0, 1, 1, 2, 2]) >>> node_z = np.array([0, 0, 0, 0, 1, 1]) - >>> _get_cartesian_face_edge_nodes( - ... face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z - ... ) + >>> _get_cartesian_faces_edge_nodes(face_node_conn,n_face,n_max_face_edges,node_x,node_y,node_z) array([[[[ 0, 0, 0], [ 1, 0, 0]], @@ -259,7 +257,7 @@ def _get_cartesian_face_edge_nodes( return face_edges_cartesian.reshape(n_face, n_max_face_edges, 2, 3) -def _get_lonlat_rad_face_edge_nodes( +def _get_lonlat_rad_faces_edge_nodes( face_node_conn, n_face, n_max_face_edges, node_lon, node_lat ): """Construct an array to hold the edge latitude and longitude in radians From b392046ab7847eb8f0081b2e47087dfe6055265d Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Thu, 13 Mar 2025 18:04:38 -0700 Subject: [PATCH 6/9] Finish all necessary changes --- test/test_integrate.py | 8 +- uxarray/core/zonal.py | 2 - uxarray/grid/geometry.py | 154 +++++++++++++++++++------------------- uxarray/grid/grid.py | 11 ++- uxarray/grid/integrate.py | 55 ++++++++++---- uxarray/grid/utils.py | 4 +- 6 files changed, 130 insertions(+), 104 deletions(-) diff --git a/test/test_integrate.py b/test/test_integrate.py index 18289e0bf..e7450d2b0 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -67,7 +67,7 @@ def test_get_faces_constLat_intersection_info_one_intersection(): latitude_cart = -0.8660254037844386 is_latlonface = False is_GCA_list = None - unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, + unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface) assert len(unique_intersections) == 1 @@ -97,7 +97,7 @@ def test_get_faces_constLat_intersection_info_encompass_pole(): is_latlonface = False is_GCA_list = None - unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, + unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface) assert len(unique_intersections) <= 2 * len(face_edges_cart) @@ -119,7 +119,7 @@ def test_get_faces_constLat_intersection_info_on_pole(): latitude_cart = -0.9998476951563913 is_latlonface = False is_GCA_list = None - unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, + unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface) assert len(unique_intersections) == 2 @@ -141,7 +141,7 @@ def test_get_faces_constLat_intersection_info_near_pole(): latitude_deg = np.rad2deg(latitude_rad) is_latlonface = False is_GCA_list = None - unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, + unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface) assert len(unique_intersections) == 1 diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py index 02e63e787..271d3d8a6 100644 --- a/uxarray/core/zonal.py +++ b/uxarray/core/zonal.py @@ -5,7 +5,6 @@ from uxarray.grid.integrate import _zonal_face_weights, _zonal_face_weights_robust - def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=False): """Computes the non-conservative zonal mean across one or more latitudes.""" uxgrid = uxda.uxgrid @@ -50,4 +49,3 @@ def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=Fal ) return result - diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 13027b2d3..e7ae36403 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -103,7 +103,7 @@ def _unique_points(points, tolerance=ERROR_TOLERANCE): @njit(cache=True) def _pad_closed_face_nodes( - face_node_connectivity, n_face, n_max_face_nodes, n_nodes_per_face + face_node_connectivity, n_face, n_max_face_nodes, n_nodes_per_face ): """Pads a closed array of face nodes by inserting the first element at any point a fill value is encountered. @@ -123,14 +123,14 @@ def _pad_closed_face_nodes( def _build_polygon_shells( - node_lon, - node_lat, - face_node_connectivity, - n_face, - n_max_face_nodes, - n_nodes_per_face, - projection=None, - central_longitude=0.0, + node_lon, + node_lat, + face_node_connectivity, + n_face, + n_max_face_nodes, + n_nodes_per_face, + projection=None, + central_longitude=0.0, ): """Builds an array of polygon shells, which can be used with Shapely to construct polygons.""" @@ -266,7 +266,7 @@ def _grid_to_polygon_geodataframe(grid, periodic_elements, projection, project, def _build_geodataframe_without_antimeridian( - polygon_shells, projected_polygon_shells, antimeridian_face_indices, engine + polygon_shells, projected_polygon_shells, antimeridian_face_indices, engine ): """Builds a ``spatialpandas.GeoDataFrame`` or ``geopandas.GeoDataFrame``excluding any faces that cross the @@ -295,10 +295,10 @@ def _build_geodataframe_without_antimeridian( def _build_geodataframe_with_antimeridian( - polygon_shells, - projected_polygon_shells, - antimeridian_face_indices, - engine, + polygon_shells, + projected_polygon_shells, + antimeridian_face_indices, + engine, ): """Builds a ``spatialpandas.GeoDataFrame`` or ``geopandas.GeoDataFrame`` including any faces that cross the antimeridian.""" @@ -317,9 +317,9 @@ def _build_geodataframe_with_antimeridian( def _build_corrected_shapely_polygons( - polygon_shells, - projected_polygon_shells, - antimeridian_face_indices, + polygon_shells, + projected_polygon_shells, + antimeridian_face_indices, ): if projected_polygon_shells is not None: # use projected shells if a projection is applied @@ -438,7 +438,7 @@ def _build_corrected_polygon_shells(polygon_shells): def _grid_to_matplotlib_polycollection( - grid, periodic_elements, projection=None, **kwargs + grid, periodic_elements, projection=None, **kwargs ): """Constructs and returns a ``matplotlib.collections.PolyCollection``""" @@ -642,7 +642,7 @@ def _get_polygons(grid, periodic_elements, projection=None, apply_projection=Tru def _grid_to_matplotlib_linecollection( - grid, periodic_elements, projection=None, **kwargs + grid, periodic_elements, projection=None, **kwargs ): """Constructs and returns a ``matplotlib.collections.LineCollection``""" @@ -688,6 +688,7 @@ def _convert_shells_to_polygons(shells): return polygons +# TODO: If we don't actually use this one, remove it def _pole_point_inside_polygon_cartesian(pole, face_edges_xyz): if isinstance(pole, str): pole = POLE_NAME_TO_INT[pole] @@ -913,7 +914,7 @@ def _check_intersection(ref_edge_xyz, edges_xyz): for i in range(n_edges): edge_xyz = edges_xyz[i] if allclose( - intersection_point, edge_xyz[0], atol=ERROR_TOLERANCE + intersection_point, edge_xyz[0], atol=ERROR_TOLERANCE ) or allclose(intersection_point, edge_xyz[1], atol=ERROR_TOLERANCE): return 0 @@ -1140,10 +1141,10 @@ def insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True): @njit(cache=True) def _populate_face_latlon_bound( - face_edges_xyz, - face_edges_lonlat, - is_latlonface=False, - is_GCA_list=None, + face_edges_xyz, + face_edges_lonlat, + is_latlonface=False, + is_GCA_list=None, ): # Check if face_edges contains pole points has_north_pole = pole_point_inside_polygon(1, face_edges_xyz, face_edges_lonlat) @@ -1198,9 +1199,9 @@ def _populate_face_latlon_bound( # Check if the node matches the pole point or if the pole point is within the edge max_abs_diff = np.max(np.abs(n1_cart - pole_point_xyz)) if max_abs_diff <= ERROR_TOLERANCE or point_within_gca( - pole_point_xyz, - n1_cart, - n2_cart, + pole_point_xyz, + n1_cart, + n2_cart, ): is_center_pole = False face_latlon_array = insert_pt_in_latlonbox( @@ -1292,15 +1293,15 @@ def _populate_face_latlon_bound( # Insert extreme latitude points into the latlonbox if ( - abs(node1_lat_rad - lat_max) > ERROR_TOLERANCE - and abs(node2_lat_rad - lat_max) > ERROR_TOLERANCE + abs(node1_lat_rad - lat_max) > ERROR_TOLERANCE + and abs(node2_lat_rad - lat_max) > ERROR_TOLERANCE ): face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([lat_max, node1_lon_rad]) ) elif ( - abs(node1_lat_rad - lat_min) > ERROR_TOLERANCE - and abs(node2_lat_rad - lat_min) > ERROR_TOLERANCE + abs(node1_lat_rad - lat_min) > ERROR_TOLERANCE + and abs(node2_lat_rad - lat_min) > ERROR_TOLERANCE ): face_latlon_array = insert_pt_in_latlonbox( face_latlon_array, np.array([lat_min, node1_lon_rad]) @@ -1315,23 +1316,23 @@ def _populate_face_latlon_bound( @njit(cache=True, parallel=True) def compute_temp_latlon_array( - face_node_connectivity, - faces_edges_cartesian, - faces_edges_lonlat_rad, - n_nodes_per_face, - is_latlonface, - is_face_GCA_list, - INT_FILL_VALUE, + face_node_connectivity, + faces_edges_cartesian, + faces_edges_lonlat_rad, + n_nodes_per_face, + is_latlonface, + is_face_GCA_list, + INT_FILL_VALUE, ): n_face = face_node_connectivity.shape[0] temp_latlon_array = np.full((n_face, 2, 2), INT_FILL_VALUE, dtype=np.float64) for face_idx in prange(n_face): cur_face_edges_cartesian = faces_edges_cartesian[ - face_idx, 0: n_nodes_per_face[face_idx] - ] + face_idx, 0 : n_nodes_per_face[face_idx] + ] cur_faces_edges_lonlat_rad = faces_edges_lonlat_rad[ - face_idx, 0: n_nodes_per_face[face_idx] - ] + face_idx, 0 : n_nodes_per_face[face_idx] + ] if is_face_GCA_list is not None: is_GCA_list = is_face_GCA_list[face_idx] else: @@ -1347,9 +1348,14 @@ def compute_temp_latlon_array( def _populate_faces_edges_cartesian(grid, return_array=False): - faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, - grid.n_max_face_edges, grid.node_x.values, - grid.node_y.values, grid.node_z.values) + faces_edges_cartesian = _get_cartesian_faces_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) faces_edges_cartesian_xarray = xr.DataArray( faces_edges_cartesian, @@ -1358,7 +1364,7 @@ def _populate_faces_edges_cartesian(grid, return_array=False): "cf_role": "faces_edges_cartesian", "_FillValue": INT_FILL_VALUE, "long_name": "Provide the Cartesian coordinates of the edge for each face", - "start_index": INT_DTYPE(0) + "start_index": INT_DTYPE(0), }, ) @@ -1369,9 +1375,13 @@ def _populate_faces_edges_cartesian(grid, return_array=False): def _populate_faces_edges_spherical(grid, return_array=False): - faces_edges_lonlat_rad = _get_lonlat_rad_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face, - grid.n_max_face_edges, grid.node_lon.values, - grid.node_lat.values) + faces_edges_lonlat_rad = _get_lonlat_rad_faces_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_lon.values, + grid.node_lat.values, + ) faces_edges_lonlat_rad_xarray = xr.DataArray( faces_edges_lonlat_rad, @@ -1380,7 +1390,7 @@ def _populate_faces_edges_spherical(grid, return_array=False): "cf_role": "faces_edges_spherical", "_FillValue": INT_FILL_VALUE, "long_name": "Provide the Spherical coordinates (lon lat) in radians of the edge for each face", - "start_index": INT_DTYPE(0) + "start_index": INT_DTYPE(0), }, ) @@ -1391,10 +1401,7 @@ def _populate_faces_edges_spherical(grid, return_array=False): def _populate_bounds( - grid, - is_latlonface: bool = False, - is_face_GCA_list=None, - return_array=False + grid, is_latlonface: bool = False, is_face_GCA_list=None, return_array=False ): """Populates the bounds of the grid based on the geometry of its faces, taking into account special conditions such as faces crossing the @@ -1408,17 +1415,6 @@ def _populate_bounds( The grid object whose internal representation will be updated with the computed face bounds. - faces_edges_cartesian : np.ndarray, optional - An array of shape (n_face, n_max_face_edges, 2, 3) containing the Cartesian coordinates - of each face's edges. This array may include dummy values for grids with holes. - Default is None. - - faces_edges_lonlat_rad : np.ndarray, optional - An array of shape (n_face, n_max_face_edges, 2, 2) containing the longitude and latitude - (in radians) coordinates of each face's edges. This array may include dummy values for - grids with holes. - Default is None. - is_latlonface : bool, optional A global flag that indicates if faces are latlon faces. If True, all faces are treated as latlon faces, meaning that all edges are either longitude or @@ -1563,16 +1559,16 @@ def stereographic_projection(lon, lat, central_lon, central_lat): # Calculate constant used for calculation k = 2.0 / ( - 1.0 - + np.sin(central_lat) * np.sin(lat) - + np.cos(central_lat) * np.cos(lat) * np.cos(lon - central_lon) + 1.0 + + np.sin(central_lat) * np.sin(lat) + + np.cos(central_lat) * np.cos(lat) * np.cos(lon - central_lon) ) # Calculate the x and y coordinates x = k * np.cos(lat) * np.sin(lon - central_lon) y = k * ( - np.cos(central_lat) * np.sin(lat) - - np.sin(central_lat) * np.cos(lat) * np.cos(lon - central_lon) + np.cos(central_lat) * np.sin(lat) + - np.sin(central_lat) * np.cos(lat) * np.cos(lon - central_lon) ) return x, y @@ -1608,7 +1604,7 @@ def inverse_stereographic_projection(x, y, central_lon, central_lat): central_lat = np.deg2rad(central_lat) # Calculate constants used for calculation - p = np.sqrt(x ** 2 + y ** 2) + p = np.sqrt(x**2 + y**2) c = 2 * np.arctan(p / 2) @@ -1627,9 +1623,9 @@ def inverse_stereographic_projection(x, y, central_lon, central_lat): @njit(cache=True) def point_in_face( - edges_xyz, - point_xyz, - inclusive=True, + edges_xyz, + point_xyz, + inclusive=True, ): """Determines if a point lies inside a face. @@ -1679,9 +1675,9 @@ def point_in_face( for ind in range(len(edges_xyz)): # If the point lies on an edge, return True if inclusive if point_within_gca( - point_xyz, - edges_xyz[ind][0], - edges_xyz[ind][1], + point_xyz, + edges_xyz[ind][0], + edges_xyz[ind][1], ): if inclusive: return True @@ -1774,7 +1770,7 @@ def _populate_max_face_radius(self): @njit(cache=True) def calculate_max_face_radius( - face_node_connectivity, node_lats_rad, node_lons_rad, face_lats_rad, face_lons_rad + face_node_connectivity, node_lats_rad, node_lons_rad, face_lats_rad, face_lons_rad ): """Finds the max face radius in the mesh. Parameters @@ -1848,7 +1844,7 @@ def haversine_distance(lon_a, lat_a, lon_b, lat_b): # Haversine formula equation_in_sqrt = (np.sin(dlat / 2) ** 2) + np.cos(lat_a) * np.cos(lat_b) * ( - np.sin(dlon / 2) ** 2 + np.sin(dlon / 2) ** 2 ) distance = 2 * np.arcsin(np.sqrt(equation_in_sqrt)) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 5596bfe15..3d0acaa0a 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -2665,9 +2665,14 @@ def get_faces_containing_point( # Get the faces in terms of their edges # Since this is a new subset, the cartesian_faces_edge_nodes need to be recalculated anyway - face_edge_nodes_xyz = _get_cartesian_faces_edge_nodes(subset.face_node_connectivity.values, subset.n_face, - subset.n_max_face_nodes, subset.node_x.values, - subset.node_y.values, subset.node_z.values) + face_edge_nodes_xyz = _get_cartesian_faces_edge_nodes( + subset.face_node_connectivity.values, + subset.n_face, + subset.n_max_face_nodes, + subset.node_x.values, + subset.node_y.values, + subset.node_z.values, + ) # Get the original face indices from the subset inverse_indices = subset.inverse_indices.face.values diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index c1ebcb9f9..4ac4dba11 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -215,11 +215,13 @@ def _get_zonal_face_interval( face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1] try: - - unique_intersections_cart,unique_intersection_lonlat, pt_lon_min, pt_lon_max = ( - _get_faces_constLat_intersection_info( - face_edges_cart, latitude_cart, is_GCA_list, is_latlonface - ) + ( + unique_intersections_cart, + unique_intersection_lonlat, + pt_lon_min, + pt_lon_max, + ) = _get_faces_constLat_intersection_info( + face_edges_cart, latitude_cart, is_GCA_list, is_latlonface ) # If there's exactly one intersection, the face is only "touched" @@ -227,9 +229,7 @@ def _get_zonal_face_interval( return pl.DataFrame({"start": [0.0], "end": [0.0]}) # Extract unique longitudes from intersection points - longitudes = np.array( - [pt[0] for pt in unique_intersection_lonlat] - ) + longitudes = np.array([pt[0] for pt in unique_intersection_lonlat]) # Handle special wrap-around cases (crossing anti-meridian, etc.) if face_lon_bound_left >= face_lon_bound_right or ( @@ -424,16 +424,33 @@ def _get_faces_constLat_intersection_info( if len(unique_intersections_cart) == 2: unique_intersection_lonlat = np.array( - [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections_cart] + [ + _xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) + for pt in unique_intersections_cart + ] ) sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0) pt_lon_min, pt_lon_max = sorted_lonlat[:, 0] - return unique_intersections_cart,unique_intersection_lonlat, pt_lon_min, pt_lon_max + return ( + unique_intersections_cart, + unique_intersection_lonlat, + pt_lon_min, + pt_lon_max, + ) elif len(unique_intersections_cart) == 1: - return unique_intersections_cart, np.array(_xyz_to_lonlat_rad(unique_intersections_cart[0][0], - unique_intersections_cart[0][1], - unique_intersections_cart[0][2])), None, None + return ( + unique_intersections_cart, + np.array( + _xyz_to_lonlat_rad( + unique_intersections_cart[0][0], + unique_intersections_cart[0][1], + unique_intersections_cart[0][2], + ) + ), + None, + None, + ) elif len(unique_intersections_cart) != 0 and len(unique_intersections_cart) != 1: # If the unique intersections numbers is larger than n_edges * 2, then it means the face is concave if len(unique_intersections_cart) > len(valid_edges) * 2: @@ -443,7 +460,10 @@ def _get_faces_constLat_intersection_info( else: # Now return all the intersections points and the pt_lon_min, pt_lon_max unique_intersection_lonlat = np.array( - [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections_cart] + [ + _xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) + for pt in unique_intersections_cart + ] ) sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0) @@ -453,7 +473,12 @@ def _get_faces_constLat_intersection_info( np.max(sorted_lonlat[:, 0]), ) - return unique_intersections_cart, unique_intersection_lonlat, pt_lon_min, pt_lon_max + return ( + unique_intersections_cart, + unique_intersection_lonlat, + pt_lon_min, + pt_lon_max, + ) elif len(unique_intersections_cart) == 0: raise ValueError( "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py index 068479dad..13a053081 100644 --- a/uxarray/grid/utils.py +++ b/uxarray/grid/utils.py @@ -179,7 +179,9 @@ def _get_cartesian_faces_edge_nodes( >>> node_x = np.array([0, 1, 1, 0, 1, 0]) >>> node_y = np.array([0, 0, 1, 1, 2, 2]) >>> node_z = np.array([0, 0, 0, 0, 1, 1]) - >>> _get_cartesian_faces_edge_nodes(face_node_conn,n_face,n_max_face_edges,node_x,node_y,node_z) + >>> _get_cartesian_faces_edge_nodes( + ... face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z + ... ) array([[[[ 0, 0, 0], [ 1, 0, 0]], From 8621fa52601f691b404948ce743fd9819ad7257c Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Mon, 7 Apr 2025 12:09:38 -0700 Subject: [PATCH 7/9] Avoid unnecessary co-planer check and created `_intersection_within_interval` --- test/test_integrate.py | 31 +++++++++++++--------- uxarray/grid/arcs.py | 50 +++++++++++++++++++++++++++++++++++ uxarray/grid/intersections.py | 22 ++++++++------- 3 files changed, 81 insertions(+), 22 deletions(-) diff --git a/test/test_integrate.py b/test/test_integrate.py index e7450d2b0..8657a06cb 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -19,7 +19,7 @@ _get_faces_constLat_intersection_info, _zonal_face_weights, \ _zonal_face_weights_robust -from uxarray.grid.utils import _get_cartesian_faces_edge_nodes +from uxarray.grid.utils import _get_cartesian_face_edge_nodes current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -67,9 +67,9 @@ def test_get_faces_constLat_intersection_info_one_intersection(): latitude_cart = -0.8660254037844386 is_latlonface = False is_GCA_list = None - unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, + unique_intersections_cart,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface) - assert len(unique_intersections) == 1 + assert len(unique_intersections_cart) == 1 def test_get_faces_constLat_intersection_info_encompass_pole(): @@ -97,9 +97,9 @@ def test_get_faces_constLat_intersection_info_encompass_pole(): is_latlonface = False is_GCA_list = None - unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, - is_GCA_list, is_latlonface) - assert len(unique_intersections) <= 2 * len(face_edges_cart) + unique_intersections_cart,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, + is_GCA_list, is_latlonface) + assert len(unique_intersections_cart) <= 2 * len(face_edges_cart) def test_get_faces_constLat_intersection_info_on_pole(): @@ -119,9 +119,9 @@ def test_get_faces_constLat_intersection_info_on_pole(): latitude_cart = -0.9998476951563913 is_latlonface = False is_GCA_list = None - unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, + unique_intersections_cart,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface) - assert len(unique_intersections) == 2 + assert len(unique_intersections_cart) == 2 def test_get_faces_constLat_intersection_info_near_pole(): @@ -141,9 +141,9 @@ def test_get_faces_constLat_intersection_info_near_pole(): latitude_deg = np.rad2deg(latitude_rad) is_latlonface = False is_GCA_list = None - unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, + unique_intersections_cart,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface) - assert len(unique_intersections) == 1 + assert len(unique_intersections_cart) == 1 def test_get_zonal_face_interval(): @@ -1037,9 +1037,14 @@ def test_compare_zonal_weights(): for gridfile in gridfiles: uxgrid = ux.open_grid(gridfile) n_nodes_per_face = uxgrid.n_nodes_per_face.values - face_edge_nodes_xyz = _get_cartesian_faces_edge_nodes(uxgrid.face_node_connectivity.values, uxgrid.n_face, - uxgrid.n_max_face_edges, uxgrid.node_x.values, - uxgrid.node_y.values, uxgrid.node_z.values) + face_edge_nodes_xyz = _get_cartesian_face_edge_nodes( + uxgrid.face_node_connectivity.values, + uxgrid.n_face, + uxgrid.n_max_face_edges, + uxgrid.node_x.values, + uxgrid.node_y.values, + uxgrid.node_z.values, + ) bounds = uxgrid.bounds.values for i, lat in enumerate(latitudes): diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index dc53a2bd7..c4e42bcec 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -90,6 +90,55 @@ def point_within_gca(pt_xyz, gca_a_xyz, gca_b_xyz): return False +@njit(cache=True) +def _intersection_within_interval(pt_xyz, gca_a_xyz, gca_b_xyz): + """ + Check if a an intersection point, which is certain to on the same plane as the CGA, lies on a given Great Circle Arc (GCA) interval, considering the smaller arc of the circle. + Handles the anti-meridian case as well. + + Parameters + ---------- + pt_xyz : numpy.ndarray + Cartesian coordinates of the point. + gca_a_xyz : numpy.ndarray + Cartesian coordinates of the first endpoint of the Great Circle Arc. + gca_b_xyz : numpy.ndarray + Cartesian coordinates of the second endpoint of the Great Circle Arc. + + Returns + ------- + bool + True if the point lies within the specified GCA interval, False otherwise. + + Raises + ------ + ValueError + In such cases, consider breaking the GCA into two separate arcs. + + Notes + ----- + - The function assumes that the point lies on the same plane as the GCA and it also assumes that the input will not be exactly 180 degree. + - It assumes the input represents the smaller arc of the Great Circle. + - The `_angle_of_2_vectors` and `_xyz_to_lonlat_rad_scalar` functions are used for calculations. + """ + + # Check if the point lies within the Great Circle Arc interval + pt_a = gca_a_xyz - pt_xyz + pt_b = gca_b_xyz - pt_xyz + + # Use the dot product to determine the sign of the angle between pt_a and pt_b + cos_theta = np.dot(pt_a, pt_b) + + # Return True if the point lies within the interval (smaller arc) + if cos_theta < 0: + return True + elif isclose(cos_theta, 0.0, atol=MACHINE_EPSILON): + # set error tolerance to 0.0 + return True + else: + return False + + @njit(cache=True) def in_between(p, q, r) -> bool: """Determines whether the number q is between p and r. @@ -369,3 +418,4 @@ def compute_arc_length(pt_a, pt_b): delta_theta = np.atan2(cross_2d, dot_2d) return rho * abs(delta_theta) + diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 0eabb7b26..97ec9a245 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -7,6 +7,7 @@ in_between, extreme_gca_z, point_within_gca, + _intersection_within_interval ) from uxarray.utils.computing import allclose, cross, norm @@ -347,13 +348,13 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz): x2_xyz = -x1_xyz # Check intersection points - if point_within_gca(x1_xyz, w0_xyz, w1_xyz) and point_within_gca( + if _intersection_within_interval(x1_xyz, w0_xyz, w1_xyz) and _intersection_within_interval( x1_xyz, v0_xyz, v1_xyz ): res[count, :] = x1_xyz count += 1 - if point_within_gca(x2_xyz, w0_xyz, w1_xyz) and point_within_gca( + if _intersection_within_interval(x2_xyz, w0_xyz, w1_xyz) and _intersection_within_interval( x2_xyz, v0_xyz, v1_xyz ): res[count, :] = x2_xyz @@ -417,17 +418,20 @@ def gca_const_lat_intersection(gca_cart, const_z): nx, ny, nz = n - s_tilde = np.sqrt(nx**2 + ny**2 - (nx**2 + ny**2 + nz**2) * const_z**2) - p1_x = -(1.0 / (nx**2 + ny**2)) * (const_z * nx * nz + s_tilde * ny) - p2_x = -(1.0 / (nx**2 + ny**2)) * (const_z * nx * nz - s_tilde * ny) - p1_y = -(1.0 / (nx**2 + ny**2)) * (const_z * ny * nz - s_tilde * nx) - p2_y = -(1.0 / (nx**2 + ny**2)) * (const_z * ny * nz + s_tilde * nx) + nx_square_ny_square = nx**2 + ny**2 + n_square = nx_square_ny_square + nz**2 + + s_tilde = np.sqrt(nx_square_ny_square - n_square * const_z**2) + p1_x = -(1.0 / (nx_square_ny_square)) * (const_z * nx * nz + s_tilde * ny) + p2_x = -(1.0 / (nx_square_ny_square)) * (const_z * nx * nz - s_tilde * ny) + p1_y = -(1.0 / (nx_square_ny_square)) * (const_z * ny * nz - s_tilde * nx) + p2_y = -(1.0 / (nx_square_ny_square)) * (const_z * ny * nz + s_tilde * nx) p1 = np.array([p1_x, p1_y, const_z]) p2 = np.array([p2_x, p2_y, const_z]) - p1_intersects_gca = point_within_gca(p1, gca_cart[0], gca_cart[1]) - p2_intersects_gca = point_within_gca(p2, gca_cart[0], gca_cart[1]) + p1_intersects_gca = _intersection_within_interval(p1, gca_cart[0], gca_cart[1]) + p2_intersects_gca = _intersection_within_interval(p2, gca_cart[0], gca_cart[1]) if p1_intersects_gca and p2_intersects_gca: res[0] = p1 From 0b87406daba0d9c537d38487c4688fde53381fa9 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Mon, 7 Apr 2025 12:10:52 -0700 Subject: [PATCH 8/9] Run precommit --- uxarray/grid/arcs.py | 1 - uxarray/grid/intersections.py | 14 +++++++------- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index c4e42bcec..8bb12a04a 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -418,4 +418,3 @@ def compute_arc_length(pt_a, pt_b): delta_theta = np.atan2(cross_2d, dot_2d) return rho * abs(delta_theta) - diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 97ec9a245..dccf0ef06 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -7,7 +7,7 @@ in_between, extreme_gca_z, point_within_gca, - _intersection_within_interval + _intersection_within_interval, ) from uxarray.utils.computing import allclose, cross, norm @@ -348,15 +348,15 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz): x2_xyz = -x1_xyz # Check intersection points - if _intersection_within_interval(x1_xyz, w0_xyz, w1_xyz) and _intersection_within_interval( - x1_xyz, v0_xyz, v1_xyz - ): + if _intersection_within_interval( + x1_xyz, w0_xyz, w1_xyz + ) and _intersection_within_interval(x1_xyz, v0_xyz, v1_xyz): res[count, :] = x1_xyz count += 1 - if _intersection_within_interval(x2_xyz, w0_xyz, w1_xyz) and _intersection_within_interval( - x2_xyz, v0_xyz, v1_xyz - ): + if _intersection_within_interval( + x2_xyz, w0_xyz, w1_xyz + ) and _intersection_within_interval(x2_xyz, v0_xyz, v1_xyz): res[count, :] = x2_xyz count += 1 From c6e58028029eb0ade1c4c646a88b6c3d78433149 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Mon, 7 Apr 2025 12:59:44 -0700 Subject: [PATCH 9/9] fix bugs --- test/test_integrate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_integrate.py b/test/test_integrate.py index 8657a06cb..2f17c39ae 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -19,7 +19,7 @@ _get_faces_constLat_intersection_info, _zonal_face_weights, \ _zonal_face_weights_robust -from uxarray.grid.utils import _get_cartesian_face_edge_nodes +from uxarray.grid.utils import _get_cartesian_faces_edge_nodes current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -1037,7 +1037,7 @@ def test_compare_zonal_weights(): for gridfile in gridfiles: uxgrid = ux.open_grid(gridfile) n_nodes_per_face = uxgrid.n_nodes_per_face.values - face_edge_nodes_xyz = _get_cartesian_face_edge_nodes( + face_edge_nodes_xyz = _get_cartesian_faces_edge_nodes( uxgrid.face_node_connectivity.values, uxgrid.n_face, uxgrid.n_max_face_edges,