From db7a7ccd80266b6b33d81a4ec1bfa5d4966d7db8 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Tue, 2 Apr 2024 22:41:51 -0500 Subject: [PATCH 01/14] initial work on vectorization --- test/test_coords.py | 31 ++++++++ uxarray/grid/veccoordinates.py | 131 +++++++++++++++++++++++++++++++++ 2 files changed, 162 insertions(+) create mode 100644 test/test_coords.py create mode 100644 uxarray/grid/veccoordinates.py diff --git a/test/test_coords.py b/test/test_coords.py new file mode 100644 index 000000000..8e3fcc5ec --- /dev/null +++ b/test/test_coords.py @@ -0,0 +1,31 @@ +import uxarray as ux + +from uxarray.constants import INT_FILL_VALUE +import numpy.testing as nt +import os + +import pytest + +from pathlib import Path + +current_path = Path(os.path.dirname(os.path.realpath(__file__))) + +GRID_PATHS = [ + current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc", + current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" +] + + + + +def test_lonlat_to_xyz(): + from uxarray.grid.veccoordinates import _populate_node_xyz + + for grid_path in GRID_PATHS: + uxgrid = ux.open_grid(grid_path) + + _populate_node_xyz(uxgrid) + + +def test_xyz_to_lonlat(): + pass diff --git a/uxarray/grid/veccoordinates.py b/uxarray/grid/veccoordinates.py new file mode 100644 index 000000000..b40d8d644 --- /dev/null +++ b/uxarray/grid/veccoordinates.py @@ -0,0 +1,131 @@ +import xarray as xr +import numpy as np + +from numba import config + +from uxarray.constants import ENABLE_JIT, ERROR_TOLERANCE +from uxarray.conventions import ugrid + +config.DISABLE_JIT = not ENABLE_JIT + + +# ====================================================================================================================== +# Conversion Functions +# ====================================================================================================================== +def _lonlat_rad_to_xyz( + lon: np.ndarray, lat: np.ndarray +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Docstring TODO.""" + 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: np.ndarray, y: np.ndarray, z: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """Docstring TODO.""" + denom = np.abs(x * x + y * y + z * z) + + x /= denom + y /= denom + z /= denom + + lon = np.empty_like(x) + lat = np.empty_like(x) + + # is there a way to vectorize this loop? + for i in range(len(x)): + if np.absolute(z[i]) < (1.0 - ERROR_TOLERANCE): + lon[i] = np.arctan2(y[i], x[i]) + lat[i] = np.arcsin(z[i]) + + if lon[i] < 0.0: + lon[i] += 2.0 * np.pi + elif z[i] > 0.0: + lon[i] = 0.0 + lat[i] = 0.5 * np.pi + else: + lon[i] = 0.0 + lat[i] = -0.5 * np.pi + + return lon, lat + + +# don't use numba here unless we figure out ord=2 +def _normalize_xyz( + x: np.ndarray, y: np.ndarray, z: np.ndarray +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Docstring TODO.""" + x_norm = x / np.linalg.norm(x, ord=2) + y_norm = y / np.linalg.norm(y, ord=2) + z_norm = z / np.linalg.norm(z, ord=2) + return x_norm, y_norm, z_norm + + +# ====================================================================================================================== +# Population Functions +# ====================================================================================================================== +def _populate_node_latlon(grid) -> None: + """Docstring TODO.""" + x_norm, y_norm, z_norm = _normalize_xyz( + grid.node_x.values, grid.node_y.values, grid.node_z.values + ) + lon_rad, lat_rad = _xyz_to_lonlat_rad(x_norm, y_norm, z_norm) + + # Convert to degrees + lon = np.rad2deg(lon_rad) + lat = np.rad2deg(lat_rad) + + grid._ds["node_lon"] = xr.DataArray( + data=lon, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LON_ATTRS + ) + grid._ds["node_lat"] = xr.DataArray( + data=lat, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LAT_ATTRS + ) + + +def _populate_edge_latlon(grid) -> None: + pass + + +def _populate_face_latlon(grid) -> None: + pass + + +def _populate_node_xyz(grid) -> None: + """Docstring TODO.""" + + node_lon_rad = np.deg2rad(grid.node_lon.values) + node_lat_rad = np.deg2rad(grid.node_lat.values) + x, y, z = _lonlat_rad_to_xyz(node_lon_rad, node_lat_rad) + + grid._ds["node_x"] = xr.DataArray( + data=x, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_X_ATTRS + ) + grid._ds["node_y"] = xr.DataArray( + data=y, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_Y_ATTRS + ) + grid._ds["node_z"] = xr.DataArray( + data=z, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_Z_ATTRS + ) + + +def _populate_edge_xyz(grid) -> None: + pass + + +def _populate_face_xyz(grid) -> None: + pass + + +# ====================================================================================================================== +# Build Functions +# ====================================================================================================================== +def _build_face_centroids(): + pass + + +def _build_edge_centroids(): + pass From cc3c81c3ab7d719726347669fc930bf11551fd11 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 8 May 2024 02:42:51 -0500 Subject: [PATCH 02/14] start overhaul of code to support new functions --- test/test_grid.py | 4 +- test/test_intersections.py | 8 +- test/test_neighbors.py | 2 - uxarray/grid/arcs.py | 615 ++++++++++++++-------------- uxarray/grid/area.py | 10 +- uxarray/grid/coordinates.py | 357 ++++++----------- uxarray/grid/grid.py | 14 +- uxarray/grid/integrate.py | 705 +++++++++++++++++---------------- uxarray/grid/veccoordinates.py | 131 ------ uxarray/io/_exodus.py | 7 +- 10 files changed, 809 insertions(+), 1044 deletions(-) delete mode 100644 uxarray/grid/veccoordinates.py diff --git a/test/test_grid.py b/test/test_grid.py index 48aad8e86..b9be85816 100644 --- a/test/test_grid.py +++ b/test/test_grid.py @@ -10,7 +10,7 @@ from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity, _build_edge_node_connectivity -from uxarray.grid.coordinates import _populate_lonlat_coord, node_lonlat_rad_to_xyz +from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz from uxarray.constants import INT_FILL_VALUE, ERROR_TOLERANCE @@ -459,7 +459,7 @@ def test_populate_lonlat_coord(self): verts_cart = np.stack((cart_x, cart_y, cart_z), axis=1) vgrid = ux.open_grid(verts_cart, latlon=False) - _populate_lonlat_coord(vgrid) + _populate_node_latlon(vgrid) # The connectivity in `__from_vert__()` will be formed in a reverse order lon_deg, lat_deg = zip(*reversed(list(zip(lon_deg, lat_deg)))) for i in range(0, vgrid.n_node): diff --git a/test/test_intersections.py b/test/test_intersections.py index f17b3c359..6291b93b0 100644 --- a/test/test_intersections.py +++ b/test/test_intersections.py @@ -3,7 +3,9 @@ import uxarray as ux from uxarray.constants import ERROR_TOLERANCE -from uxarray.grid.coordinates import node_lonlat_rad_to_xyz, node_xyz_to_lonlat_rad +# from uxarray.grid.coordinates import node_lonlat_rad_to_xyz, node_xyz_to_lonlat_rad + +from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad from uxarray.grid.intersections import gca_gca_intersection, gca_constLat_intersection @@ -11,6 +13,10 @@ class TestGCAGCAIntersection(TestCase): def test_get_GCA_GCA_intersections_antimeridian(self): # Test the case where the two GCAs are on the antimeridian + + + + GCA1 = node_lonlat_rad_to_xyz([np.deg2rad(170.0), np.deg2rad(89.99)]) GCR1_cart = np.array([ node_lonlat_rad_to_xyz([np.deg2rad(170.0), diff --git a/test/test_neighbors.py b/test/test_neighbors.py index 87736fae4..74321c247 100644 --- a/test/test_neighbors.py +++ b/test/test_neighbors.py @@ -10,8 +10,6 @@ from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity -from uxarray.grid.coordinates import _populate_lonlat_coord - from uxarray.constants import INT_FILL_VALUE try: diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index e2c4bddfb..df4a9dc7f 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -1,303 +1,312 @@ -import numpy as np - -from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place -from uxarray.constants import ERROR_TOLERANCE - - -def _to_list(obj): - if not isinstance(obj, list): - if isinstance(obj, np.ndarray): - # Convert the NumPy array to a list using .tolist() - obj = obj.tolist() - else: - # If not a list or NumPy array, return the object as-is - obj = [obj] - return obj - - -def point_within_gca(pt, gca_cart, is_directed=False): - """Check if a point lies on a given Great Circle Arc (GCA). The anti- - meridian case is also considered. - - Parameters - ---------- - pt : numpy.ndarray (float) - Cartesian coordinates of the point. - gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr) - Cartesian coordinates of the Great Circle Arc (GCR). - is_directed : bool, optional, default = False - If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected, - and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False. - For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1 - - Returns - ------- - bool - True if the point lies between the two endpoints of the GCR, False otherwise. - - Raises - ------ - ValueError - If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes. - In such cases, consider breaking the GCR into two separate GCRs. - - ValueError - If the input GCR spans more than 180 degrees (π radians). - In such cases, consider breaking the GCR into two separate GCRs. - - Notes - ----- - The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and - accounting for the anti-meridian case. - - The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude). - In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation. - - The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations. - - Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. - """ - # Convert the cartesian coordinates to lonlat coordinates - pt_lonlat = node_xyz_to_lonlat_rad(_to_list(pt)) - GCRv0_lonlat = node_xyz_to_lonlat_rad(_to_list(gca_cart[0])) - GCRv1_lonlat = node_xyz_to_lonlat_rad(_to_list(gca_cart[1])) - - # Convert the list to np.float64 - gca_cart[0] = np.array(gca_cart[0], dtype=np.float64) - gca_cart[1] = np.array(gca_cart[1], dtype=np.float64) - - # First if the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes - angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1]) - if np.allclose(angle, np.pi, rtol=0, atol=ERROR_TOLERANCE): - raise ValueError( - "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. " - "Consider breaking the Great Circle Arc" - "into two Great Circle Arcs" - ) - - if not np.allclose( - np.dot(np.cross(gca_cart[0], gca_cart[1]), pt), 0, rtol=0, atol=ERROR_TOLERANCE - ): - return False - - if np.isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): - # If the pt and the GCA are on the same longitude (the y coordinates are the same) - if np.isclose(GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): - # Now use the latitude to determine if the pt falls between the interval - return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1]) - else: - # If the pt and the GCA are not on the same longitude when the GCA is a longnitude arc, then the pt is not on the GCA - return False - - # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point - if np.isclose( - abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0, atol=ERROR_TOLERANCE - ): - # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0] - if np.isclose(abs(pt_lonlat[1]), np.pi / 2, rtol=0, atol=ERROR_TOLERANCE): - pt_lonlat[0] = GCRv0_lonlat[0] - if not np.isclose( - GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE - ) and not np.isclose( - GCRv1_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE - ): - return False - else: - # Determine the pole latitude and latitude extension - if (GCRv0_lonlat[1] > 0 and GCRv1_lonlat[1] > 0) or ( - GCRv0_lonlat[1] < 0 and GCRv1_lonlat[1] < 0 - ): - pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0 else -np.pi / 2 - lat_extend = ( - abs(np.pi / 2 - abs(GCRv0_lonlat[1])) - + np.pi / 2 - + abs(GCRv1_lonlat[1]) - ) - else: - pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1]) - lat_extend = ( - abs(np.pi / 2 - abs(GCRv0_lonlat[1])) - + np.pi / 2 - + abs(GCRv1_lonlat[1]) - ) - - if is_directed and lat_extend >= np.pi: - raise ValueError( - "The input GCA spans more than 180 degrees. Please divide it into two." - ) - - return in_between(GCRv0_lonlat[1], pt_lonlat[1], pole_lat) or in_between( - pole_lat, pt_lonlat[1], GCRv1_lonlat[1] - ) - - if is_directed: - # The anti-meridian case Sufficient condition: absolute difference between the longitudes of the two - # vertices is greater than 180 degrees (π radians): abs(GCRv1_lon - GCRv0_lon) > π - if abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]) > np.pi: - # The necessary condition: the pt longitude is on the opposite side of the anti-meridian - # Case 1: where 0 --> x0--> 180 -->x1 -->0 case is lager than the 180degrees (pi radians) - if GCRv0_lonlat[0] <= np.pi <= GCRv1_lonlat[0]: - raise ValueError( - "The input Great Circle Arc span is larger than 180 degree, please break it into two." - ) - - # The necessary condition: the pt longitude is on the opposite side of the anti-meridian - # Case 2: The anti-meridian case where 180 -->x0 --> 0 lon --> x1 --> 180 lon - elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0: - return in_between( - GCRv0_lonlat[0], pt_lonlat[0], 2 * np.pi - ) or in_between(0, pt_lonlat[0], GCRv1_lonlat[0]) - - # The non-anti-meridian case. - else: - return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) - else: - # The undirected case - # sort the longitude - GCRv0_lonlat_min, GCRv1_lonlat_max = sorted([GCRv0_lonlat[0], GCRv1_lonlat[0]]) - if np.pi > GCRv1_lonlat_max - GCRv0_lonlat_min >= 0.0: - return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) - else: - return in_between(GCRv1_lonlat_max, pt_lonlat[0], 2 * np.pi) or in_between( - 0.0, pt_lonlat[0], GCRv0_lonlat_min - ) - - -def in_between(p, q, r) -> bool: - """Determines whether the number q is between p and r. - - Parameters - ---------- - p : float - The lower bound. - q : float - The number to check. - r : float - The upper bound. - - Returns - ------- - bool - True if q is between p and r, False otherwise. - """ - - return p <= q <= r or r <= q <= p - - -def _decide_pole_latitude(lat1, lat2): - """Determine the pole latitude based on the latitudes of two points on a - Great Circle Arc (GCA). - - This function calculates the combined latitude span from each point to its nearest pole - and decides which pole (North or South) the smaller GCA will pass. This decision is crucial - for handling GCAs that span exactly or more than 180 degrees in longitude, indicating - the arc might pass through or close to one of the Earth's poles. - - Parameters - ---------- - lat1 : float - Latitude of the first point in radians. Positive for the Northern Hemisphere, negative for the Southern. - lat2 : float - Latitude of the second point in radians. Positive for the Northern Hemisphere, negative for the Southern. - - Returns - ------- - float - The latitude of the pole (np.pi/2 for the North Pole or -np.pi/2 for the South Pole) the GCA is closer to. - - Notes - ----- - The function assumes the input latitudes are valid (i.e., between -np.pi/2 and np.pi/2) and expressed in radians. - The determination of which pole a GCA is closer to is based on the sum of the latitudinal spans from each point - to its nearest pole, considering the shortest path on the sphere. - """ - # Calculate the total latitudinal span to the nearest poles - lat_extend = abs(np.pi / 2 - abs(lat1)) + np.pi / 2 + abs(lat2) - - # Determine the closest pole based on the latitudinal span - if lat_extend < np.pi: - closest_pole = np.pi / 2 if lat1 > 0 else -np.pi / 2 - else: - closest_pole = -np.pi / 2 if lat1 > 0 else np.pi / 2 - - return closest_pole - - -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. - - Parameters - ---------- - gca_cart : numpy.ndarray - An array containing two 3D vectors that define a great circle arc. - - extreme_type : str - The type of extreme latitude to calculate. Must be either 'max' or 'min'. - - Returns - ------- - float - The maximum or minimum latitude of the great circle arc in radians. - - Raises - ------ - ValueError - If `extreme_type` is not 'max' or 'min'. - """ - extreme_type = extreme_type.lower() - - if extreme_type not in ("max", "min"): - raise ValueError("extreme_type must be either 'max' or 'min'") - - n1, n2 = gca_cart - dot_n1_n2 = np.dot(n1, n2) - denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0) - d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom - - d_a_max = ( - np.clip(d_a_max, 0, 1) - if np.isclose(d_a_max, [0, 1], atol=ERROR_TOLERANCE).any() - else d_a_max - ) - lat_n1, lat_n2 = ( - node_xyz_to_lonlat_rad(n1.tolist())[1], - node_xyz_to_lonlat_rad(n2.tolist())[1], - ) - - if 0 < d_a_max < 1: - node3 = (1 - d_a_max) * n1 + d_a_max * n2 - node3 = np.array(normalize_in_place(node3.tolist())) - d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1)) - - return ( - max(d_lat_rad, lat_n1, lat_n2) - if extreme_type == "max" - else min(d_lat_rad, lat_n1, lat_n2) - ) - else: - return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2) +# import numpy as np +# +# # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place +# from uxarray.grid.coordinates import _xyz_to_lonlat_rad +# +# +# from uxarray.constants import ERROR_TOLERANCE +# +# +# def _to_list(obj): +# if not isinstance(obj, list): +# if isinstance(obj, np.ndarray): +# # Convert the NumPy array to a list using .tolist() +# obj = obj.tolist() +# else: +# # If not a list or NumPy array, return the object as-is +# obj = [obj] +# return obj +# +# +# def point_within_gca(pt, gca_cart, is_directed=False): +# """Check if a point lies on a given Great Circle Arc (GCA). The anti- +# meridian case is also considered. +# +# Parameters +# ---------- +# pt : numpy.ndarray (float) +# Cartesian coordinates of the point. +# gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr) +# Cartesian coordinates of the Great Circle Arc (GCR). +# is_directed : bool, optional, default = False +# If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected, +# and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False. +# For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1 +# +# Returns +# ------- +# bool +# True if the point lies between the two endpoints of the GCR, False otherwise. +# +# Raises +# ------ +# ValueError +# If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes. +# In such cases, consider breaking the GCR into two separate GCRs. +# +# ValueError +# If the input GCR spans more than 180 degrees (π radians). +# In such cases, consider breaking the GCR into two separate GCRs. +# +# Notes +# ----- +# The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and +# accounting for the anti-meridian case. +# +# The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude). +# In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation. +# +# The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations. +# +# Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. +# """ +# # Convert the cartesian coordinates to lonlat coordinates +# +# pt_lon, pt_lat = _xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) +# +# GCRv0_lon, GCRv0_lat = _xyz_to_lonlat_rad( +# gca_cart[0][0], gca_cart[0][1], gca_cart[0][2] +# ) +# GCRv1_lon, GCRv1_lat = _xyz_to_lonlat_rad( +# gca_cart[1][0], gca_cart[1][1], gca_cart[1][2] +# ) +# +# # Convert the list to np.float64 +# gca_cart[0] = np.array(gca_cart[0], dtype=np.float64) +# gca_cart[1] = np.array(gca_cart[1], dtype=np.float64) +# +# # First if the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes +# angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1]) +# if np.allclose(angle, np.pi, rtol=0, atol=ERROR_TOLERANCE): +# raise ValueError( +# "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. " +# "Consider breaking the Great Circle Arc" +# "into two Great Circle Arcs" +# ) +# +# if not np.allclose( +# np.dot(np.cross(gca_cart[0], gca_cart[1]), pt), 0, rtol=0, atol=ERROR_TOLERANCE +# ): +# return False +# +# if np.isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): +# # If the pt and the GCA are on the same longitude (the y coordinates are the same) +# if np.isclose(GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): +# # Now use the latitude to determine if the pt falls between the interval +# return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1]) +# else: +# # If the pt and the GCA are not on the same longitude when the GCA is a longnitude arc, then the pt is not on the GCA +# return False +# +# # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point +# if np.isclose( +# abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0, atol=ERROR_TOLERANCE +# ): +# # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0] +# if np.isclose(abs(pt_lonlat[1]), np.pi / 2, rtol=0, atol=ERROR_TOLERANCE): +# pt_lonlat[0] = GCRv0_lonlat[0] +# if not np.isclose( +# GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE +# ) and not np.isclose( +# GCRv1_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE +# ): +# return False +# else: +# # Determine the pole latitude and latitude extension +# if (GCRv0_lonlat[1] > 0 and GCRv1_lonlat[1] > 0) or ( +# GCRv0_lonlat[1] < 0 and GCRv1_lonlat[1] < 0 +# ): +# pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0 else -np.pi / 2 +# lat_extend = ( +# abs(np.pi / 2 - abs(GCRv0_lonlat[1])) +# + np.pi / 2 +# + abs(GCRv1_lonlat[1]) +# ) +# else: +# pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1]) +# lat_extend = ( +# abs(np.pi / 2 - abs(GCRv0_lonlat[1])) +# + np.pi / 2 +# + abs(GCRv1_lonlat[1]) +# ) +# +# if is_directed and lat_extend >= np.pi: +# raise ValueError( +# "The input GCA spans more than 180 degrees. Please divide it into two." +# ) +# +# return in_between(GCRv0_lonlat[1], pt_lonlat[1], pole_lat) or in_between( +# pole_lat, pt_lonlat[1], GCRv1_lonlat[1] +# ) +# +# if is_directed: +# # The anti-meridian case Sufficient condition: absolute difference between the longitudes of the two +# # vertices is greater than 180 degrees (π radians): abs(GCRv1_lon - GCRv0_lon) > π +# if abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]) > np.pi: +# # The necessary condition: the pt longitude is on the opposite side of the anti-meridian +# # Case 1: where 0 --> x0--> 180 -->x1 -->0 case is lager than the 180degrees (pi radians) +# if GCRv0_lonlat[0] <= np.pi <= GCRv1_lonlat[0]: +# raise ValueError( +# "The input Great Circle Arc span is larger than 180 degree, please break it into two." +# ) +# +# # The necessary condition: the pt longitude is on the opposite side of the anti-meridian +# # Case 2: The anti-meridian case where 180 -->x0 --> 0 lon --> x1 --> 180 lon +# elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0: +# return in_between( +# GCRv0_lonlat[0], pt_lonlat[0], 2 * np.pi +# ) or in_between(0, pt_lonlat[0], GCRv1_lonlat[0]) +# +# # The non-anti-meridian case. +# else: +# return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) +# else: +# # The undirected case +# # sort the longitude +# GCRv0_lonlat_min, GCRv1_lonlat_max = sorted([GCRv0_lonlat[0], GCRv1_lonlat[0]]) +# if np.pi > GCRv1_lonlat_max - GCRv0_lonlat_min >= 0.0: +# return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) +# else: +# return in_between(GCRv1_lonlat_max, pt_lonlat[0], 2 * np.pi) or in_between( +# 0.0, pt_lonlat[0], GCRv0_lonlat_min +# ) +# +# +# def in_between(p, q, r) -> bool: +# """Determines whether the number q is between p and r. +# +# Parameters +# ---------- +# p : float +# The lower bound. +# q : float +# The number to check. +# r : float +# The upper bound. +# +# Returns +# ------- +# bool +# True if q is between p and r, False otherwise. +# """ +# +# return p <= q <= r or r <= q <= p +# +# +# def _decide_pole_latitude(lat1, lat2): +# """Determine the pole latitude based on the latitudes of two points on a +# Great Circle Arc (GCA). +# +# This function calculates the combined latitude span from each point to its nearest pole +# and decides which pole (North or South) the smaller GCA will pass. This decision is crucial +# for handling GCAs that span exactly or more than 180 degrees in longitude, indicating +# the arc might pass through or close to one of the Earth's poles. +# +# Parameters +# ---------- +# lat1 : float +# Latitude of the first point in radians. Positive for the Northern Hemisphere, negative for the Southern. +# lat2 : float +# Latitude of the second point in radians. Positive for the Northern Hemisphere, negative for the Southern. +# +# Returns +# ------- +# float +# The latitude of the pole (np.pi/2 for the North Pole or -np.pi/2 for the South Pole) the GCA is closer to. +# +# Notes +# ----- +# The function assumes the input latitudes are valid (i.e., between -np.pi/2 and np.pi/2) and expressed in radians. +# The determination of which pole a GCA is closer to is based on the sum of the latitudinal spans from each point +# to its nearest pole, considering the shortest path on the sphere. +# """ +# # Calculate the total latitudinal span to the nearest poles +# lat_extend = abs(np.pi / 2 - abs(lat1)) + np.pi / 2 + abs(lat2) +# +# # Determine the closest pole based on the latitudinal span +# if lat_extend < np.pi: +# closest_pole = np.pi / 2 if lat1 > 0 else -np.pi / 2 +# else: +# closest_pole = -np.pi / 2 if lat1 > 0 else np.pi / 2 +# +# return closest_pole +# +# +# 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. +# +# Parameters +# ---------- +# gca_cart : numpy.ndarray +# An array containing two 3D vectors that define a great circle arc. +# +# extreme_type : str +# The type of extreme latitude to calculate. Must be either 'max' or 'min'. +# +# Returns +# ------- +# float +# The maximum or minimum latitude of the great circle arc in radians. +# +# Raises +# ------ +# ValueError +# If `extreme_type` is not 'max' or 'min'. +# """ +# extreme_type = extreme_type.lower() +# +# if extreme_type not in ("max", "min"): +# raise ValueError("extreme_type must be either 'max' or 'min'") +# +# n1, n2 = gca_cart +# dot_n1_n2 = np.dot(n1, n2) +# denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0) +# d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom +# +# d_a_max = ( +# np.clip(d_a_max, 0, 1) +# if np.isclose(d_a_max, [0, 1], atol=ERROR_TOLERANCE).any() +# else d_a_max +# ) +# _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2]) +# _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2]) +# +# if 0 < d_a_max < 1: +# node3 = (1 - d_a_max) * n1 + d_a_max * n2 +# # TODO: +# node3 = np.array(normalize_in_place(node3.tolist())) +# node3 = np.array(normalize) +# d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1)) +# +# return ( +# max(d_lat_rad, lat_n1, lat_n2) +# if extreme_type == "max" +# else min(d_lat_rad, lat_n1, lat_n2) +# ) +# else: +# return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2) diff --git a/uxarray/grid/area.py b/uxarray/grid/area.py index 086a22334..5afa72e69 100644 --- a/uxarray/grid/area.py +++ b/uxarray/grid/area.py @@ -2,7 +2,7 @@ from numba import njit, config from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT -from uxarray.grid.coordinates import node_lonlat_rad_to_xyz +from uxarray.grid.coordinates import _lonlat_rad_to_xyz config.DISABLE_JIT = not ENABLE_JIT @@ -67,14 +67,12 @@ def calculate_face_area( node3 = np.array([x[j + 2], y[j + 2], z[j + 2]], dtype=np.float64) if coords_type == "spherical": - node1 = np.array( - node_lonlat_rad_to_xyz([np.deg2rad(x[0]), np.deg2rad(y[0])]) - ) + node1 = np.array(_lonlat_rad_to_xyz(np.deg2rad(x[0]), np.deg2rad(y[0]))) node2 = np.array( - node_lonlat_rad_to_xyz([np.deg2rad(x[j + 1]), np.deg2rad(y[j + 1])]) + _lonlat_rad_to_xyz(np.deg2rad(x[j + 1]), np.deg2rad(y[j + 1])) ) node3 = np.array( - node_lonlat_rad_to_xyz([np.deg2rad(x[j + 2]), np.deg2rad(y[j + 2])]) + _lonlat_rad_to_xyz(np.deg2rad(x[j + 2]), np.deg2rad(y[j + 2])) ) for p in range(len(dW)): diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 2aba490ea..228c206f6 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -1,154 +1,123 @@ -import warnings - import xarray as xr import numpy as np -import math - -from numba import njit, config -from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT, ERROR_TOLERANCE +from uxarray.constants import ERROR_TOLERANCE from uxarray.conventions import ugrid -config.DISABLE_JIT = not ENABLE_JIT +from typing import Union + + +def _lonlat_rad_to_xyz( + lon: Union[np.ndarray, float], + lat: Union[np.ndarray, float], +) -> tuple[ + Union[np.ndarray, float], Union[np.ndarray, float], Union[np.ndarray, float] +]: + """Converts Spherical latitude and longitude 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], + scalar: bool = False, +) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: + """Converts Cartesian x, y, z coordinates in Spherical latitude and + longitude coordinates.""" + denom = np.abs(x * x + y * y + z * z) + + x /= denom + y /= denom + z /= denom + + lon = np.arctan2(y, x) + lat = np.arcsin(z) + + # set longitude range to [0, pi] + lon = np.mod(lon, 2 * np.pi) + + if scalar: + if np.abs(z) > 1.0 - ERROR_TOLERANCE: + lat = np.sign(z) * np.pi / 2 + else: + # adjust z values near +- 1 + lat = np.where(np.abs(z) > 1.0 - ERROR_TOLERANCE, np.sign(z) * np.pi / 2, lat) + return lon, lat -@njit(cache=ENABLE_JIT_CACHE) -def node_lonlat_rad_to_xyz(node_coord): - """Helper function to Convert the node coordinate from 2D - longitude/latitude to normalized 3D xyz. - Parameters - ---------- - node: float list - 2D coordinates[longitude, latitude] in radiance +def _xyz_to_lonlat_deg( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], +) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: + lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z) - Returns - ---------- - float list - the result array of the unit 3D coordinates [x, y, z] vector where :math:`x^2 + y^2 + z^2 = 1` + # set lon range to [-pi, pi] + lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi - Raises - ---------- - RuntimeError - The input array doesn't have the size of 3. - """ - if len(node_coord) != 2: - raise RuntimeError( - "Input array should have a length of 2: [longitude, latitude]" - ) - lon = node_coord[0] - lat = node_coord[1] - return [np.cos(lon) * np.cos(lat), np.sin(lon) * np.cos(lat), np.sin(lat)] + lon = np.rad2deg(lon_rad) + lat = np.rad2deg(lat_rad) + return lon, lat -@njit(cache=ENABLE_JIT_CACHE) -def node_xyz_to_lonlat_rad(node_coord): - """Calculate the latitude and longitude in radiance for a node represented - in the [x, y, z] 3D Cartesian coordinates. +def _normalize_xyz( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], + scalar: bool = False, +) -> tuple[ + Union[np.ndarray, float], Union[np.ndarray, float], Union[np.ndarray, float] +]: + denom = np.linalg.norm([x, y, z], ord=2) - Parameters - ---------- - node_coord: float list - 3D Cartesian Coordinates [x, y, z] of the node + x_norm = x / denom + y_norm = y / denom + z_norm = z / denom + return x_norm, y_norm, z_norm - Returns - ---------- - float list - the result array of longitude and latitude in radian [longitude_rad, latitude_rad] - Raises - ---------- - RuntimeError - The input array doesn't have the size of 3. - """ - if len(node_coord) != 3: - raise RuntimeError("Input array should have a length of 3: [x, y, z]") - - [dx, dy, dz] = normalize_in_place(node_coord) - dx /= np.absolute(dx * dx + dy * dy + dz * dz) - dy /= np.absolute(dx * dx + dy * dy + dz * dz) - dz /= np.absolute(dx * dx + dy * dy + dz * dz) - - if np.absolute(dz) < (1.0 - ERROR_TOLERANCE): - d_lon_rad = math.atan2(dy, dx) - d_lat_rad = np.arcsin(dz) - - if d_lon_rad < 0.0: - d_lon_rad += 2.0 * np.pi - elif dz > 0.0: - d_lon_rad = 0.0 - d_lat_rad = 0.5 * np.pi - else: - d_lon_rad = 0.0 - d_lat_rad = -0.5 * np.pi +def _populate_node_latlon(grid) -> None: + """Docstring TODO.""" + x_norm, y_norm, z_norm = _normalize_xyz( + grid.node_x.values, grid.node_y.values, grid.node_z.values + ) + lon_rad, lat_rad = _xyz_to_lonlat_rad(x_norm, y_norm, z_norm) - return [d_lon_rad, d_lat_rad] + # set longitude range to [-pi, pi] + lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi + # convert to degrees + lon = np.rad2deg(lon_rad) + lat = np.rad2deg(lat_rad) -@njit(cache=ENABLE_JIT_CACHE) -def normalize_in_place(node): - """Helper function to project an arbitrary node in 3D coordinates [x, y, z] - on the unit sphere. It uses the `np.linalg.norm` internally to calculate - the magnitude. + grid._ds["node_lon"] = xr.DataArray( + data=lon, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LON_ATTRS + ) + grid._ds["node_lat"] = xr.DataArray( + data=lat, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LAT_ATTRS + ) - Parameters - ---------- - node: float list - 3D Cartesian Coordinates [x, y, z] - Returns - ---------- - float list - the result unit vector [x, y, z] where :math:`x^2 + y^2 + z^2 = 1` +def _populate_edge_latlon(grid) -> None: + pass - Raises - ---------- - RuntimeError - The input array doesn't have the size of 3. - """ - if len(node) != 3: - raise RuntimeError("Input array should have a length of 3: [x, y, z]") - return list(np.array(node) / np.linalg.norm(np.array(node), ord=2)) - - -def _get_xyz_from_lonlat(node_lon, node_lat): - nodes_lon_rad = np.deg2rad(node_lon) - nodes_lat_rad = np.deg2rad(node_lat) - nodes_rad = np.stack((nodes_lon_rad, nodes_lat_rad), axis=1) - nodes_cart = np.asarray(list(map(node_lonlat_rad_to_xyz, list(nodes_rad)))) - - return nodes_cart[:, 0], nodes_cart[:, 1], nodes_cart[:, 2] - - -def _populate_cartesian_xyz_coord(grid): - """A helper function that populates the xyz attribute in UXarray.Grid._ds. - This function is called when we need to use the cartesian coordinates for - each node to do the calculation but the input data only has the - ``node_lon`` and ``node_lat`` in degree. - - Note - ---- - In the UXarray, we abide the UGRID convention and make sure the following attributes will always have its - corresponding units as stated below: - - node_lon - unit: "degree_east" for longitude - node_lat - unit: "degrees_north" for latitude - node_x - unit: "m" - node_y - unit: "m" - node_z - unit: "m" - """ - # Check if the cartesian coordinates are already populated - if "node_x" in grid._ds.keys(): - return +def _populate_face_latlon(grid) -> None: + pass + - # get Cartesian (x, y, z) coordinates from lon/lat - x, y, z = _get_xyz_from_lonlat(grid.node_lon.values, grid.node_lat.values) +def _populate_node_xyz(grid) -> None: + """Docstring TODO.""" + + node_lon_rad = np.deg2rad(grid.node_lon.values) + node_lat_rad = np.deg2rad(grid.node_lat.values) + x, y, z = _lonlat_rad_to_xyz(node_lon_rad, node_lat_rad) grid._ds["node_x"] = xr.DataArray( data=x, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_X_ATTRS @@ -161,65 +130,29 @@ def _populate_cartesian_xyz_coord(grid): ) -def _get_lonlat_from_xyz(x, y, z): - nodes_cart = np.stack((x, y, z), axis=1).tolist() - nodes_rad = list(map(node_xyz_to_lonlat_rad, nodes_cart)) - nodes_degree = np.rad2deg(nodes_rad) - - return nodes_degree[:, 0], nodes_degree[:, 1] - - -def _populate_lonlat_coord(grid): - """Helper function that populates the longitude and latitude and store it - into the ``node_lon`` and ``node_lat``. This is called when the input data - has ``node_x``, ``node_y``, ``node_z`` in meters. Since we want - ``node_lon`` and ``node_lat`` to always have the "degree" units. For more - details, please read the following. - - Note - ---- - In the UXarray, we abide the UGRID convention and make sure the following attributes will always have its - corresponding units as stated below: - - node_lon - unit: "degree_east" for longitude - node_lat - unit: "degrees_north" for latitude - node_x - unit: "m" - node_y - unit: "m" - node_z - unit: "m" - """ +def _build_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face): + 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) - # get lon/lat coordinates from Cartesian (x, y, z) - lon, lat = _get_lonlat_from_xyz( - grid.node_x.values, grid.node_y.values, grid.node_z.values - ) + 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]]) - # populate dataset - grid._ds["node_lon"] = xr.DataArray( - data=lon, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LON_ATTRS - ) - grid._ds["node_lat"] = xr.DataArray( - data=lat, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LAT_ATTRS - ) + return _normalize_xyz(centroid_x, centroid_y, centroid_z) -def _populate_face_centroids(grid, repopulate=False): - """Finds the centroids of faces using cartesian averaging based off the - vertices. The centroid is defined as the average of the x, y, z - coordinates, normalized. This cannot be guaranteed to work on concave - polygons. +def _build_edge_centroids(node_x, node_y, node_z, edge_nodes): + centroid_x = np.mean(node_x[edge_nodes]) + centroid_y = np.mean(node_y[edge_nodes]) + centroid_z = np.mean(node_z[edge_nodes]) + + return _normalize_xyz(centroid_x, centroid_y, centroid_z) - Parameters - ---------- - repopulate : bool, optional - Bool used to turn on/off repopulating the face coordinates of the centroids - """ - warnings.warn("This cannot be guaranteed to work correctly on concave polygons") +def _populate_face_centroids(grid, repopulate=False): node_x = grid.node_x.values node_y = grid.node_y.values node_z = grid.node_z.values @@ -229,7 +162,7 @@ def _populate_face_centroids(grid, repopulate=False): if "face_lon" not in grid._ds or repopulate: # Construct the centroids if there are none stored if "face_x" not in grid._ds: - centroid_x, centroid_y, centroid_z = _construct_face_centroids( + centroid_x, centroid_y, centroid_z = _build_face_centroids( node_x, node_y, node_z, face_nodes, n_nodes_per_face ) @@ -238,13 +171,13 @@ def _populate_face_centroids(grid, repopulate=False): centroid_x, centroid_y, centroid_z = grid.face_x, grid.face_y, grid.face_z # Convert from xyz to latlon - centroid_lon, centroid_lat = _get_lonlat_from_xyz( + centroid_lon, centroid_lat = _xyz_to_lonlat_rad( centroid_x, centroid_y, centroid_z ) else: # Convert to xyz if there are latlon centroids already stored centroid_lon, centroid_lat = grid.face_lon.values, grid.face_lat.values - centroid_x, centroid_y, centroid_z = _get_xyz_from_lonlat( + centroid_x, centroid_y, centroid_z = _lonlat_rad_to_xyz( centroid_lon, centroid_lat ) @@ -271,31 +204,6 @@ def _populate_face_centroids(grid, repopulate=False): ) -@njit() -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.""" - centroids = np.zeros((3, face_nodes.shape[0]), dtype=np.float64) - - for face_idx, n_max_nodes in enumerate(n_nodes_per_face): - # compute cartesian average - centroid_x = np.mean(node_x[face_nodes[face_idx, 0:n_max_nodes]]) - centroid_y = np.mean(node_y[face_nodes[face_idx, 0:n_max_nodes]]) - centroid_z = np.mean(node_z[face_nodes[face_idx, 0:n_max_nodes]]) - - # normalize coordinates - centroid_normalized_xyz = normalize_in_place( - [centroid_x, centroid_y, centroid_z] - ) - - # store xyz - centroids[0, face_idx] = centroid_normalized_xyz[0] - centroids[1, face_idx] = centroid_normalized_xyz[1] - centroids[2, face_idx] = centroid_normalized_xyz[2] - - return centroids[0, :], centroids[1, :], centroids[2, :] - - 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 @@ -315,7 +223,7 @@ def _populate_edge_centroids(grid, repopulate=False): if "edge_lon" not in grid._ds or repopulate: # Construct the centroids if there are none stored if "edge_x" not in grid._ds: - centroid_x, centroid_y, centroid_z = _construct_edge_centroids( + centroid_x, centroid_y, centroid_z = _build_edge_centroids( node_x, node_y, node_z, edge_nodes_con ) @@ -324,13 +232,13 @@ def _populate_edge_centroids(grid, repopulate=False): centroid_x, centroid_y, centroid_z = grid.edge_x, grid.edge_y, grid.edge_z # Convert from xyz to latlon - centroid_lon, centroid_lat = _get_lonlat_from_xyz( + centroid_lon, centroid_lat = _xyz_to_lonlat_rad( centroid_x, centroid_y, centroid_z ) else: # Convert to xyz if there are latlon centroids already stored centroid_lon, centroid_lat = grid.edge_lon.values, grid.edge_lat.values - centroid_x, centroid_y, centroid_z = _get_xyz_from_lonlat( + centroid_x, centroid_y, centroid_z = _lonlat_rad_to_xyz( centroid_lon, centroid_lat ) @@ -365,31 +273,6 @@ def _populate_edge_centroids(grid, repopulate=False): ) -@njit() -def _construct_edge_centroids(node_x, node_y, node_z, edge_nodes_con): - """Constructs the xyz centroid coordinate for each edge using Cartesian - Averaging.""" - centroids = np.zeros((3, edge_nodes_con.shape[0]), dtype=np.float64) - - for edge_idx, connectivity in enumerate(edge_nodes_con): - # compute cartesian average - centroid_x = np.mean(node_x[connectivity[0:]]) - centroid_y = np.mean(node_y[connectivity[0:]]) - centroid_z = np.mean(node_z[connectivity[0:]]) - - # normalize coordinates - centroid_normalized_xyz = normalize_in_place( - [centroid_x, centroid_y, centroid_z] - ) - - # store xyz - centroids[0, edge_idx] = centroid_normalized_xyz[0] - centroids[1, edge_idx] = centroid_normalized_xyz[1] - centroids[2, edge_idx] = centroid_normalized_xyz[2] - - return centroids[0, :], centroids[1, :], centroids[2, :] - - def _set_desired_longitude_range(ds): """Sets the longitude range to [-180, 180] for all longitude variables.""" diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index ffb0fdae5..716a1983b 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -28,8 +28,8 @@ _populate_face_centroids, _populate_edge_centroids, _set_desired_longitude_range, - _populate_lonlat_coord, - _populate_cartesian_xyz_coord, + _populate_node_latlon, + _populate_node_xyz, ) from uxarray.grid.connectivity import ( _populate_edge_node_connectivity, @@ -554,7 +554,7 @@ def node_lon(self) -> xr.DataArray: """ if "node_lon" not in self._ds: _set_desired_longitude_range(self._ds) - _populate_lonlat_coord(self) + _populate_node_latlon(self) return self._ds["node_lon"] @property @@ -565,7 +565,7 @@ def node_lat(self) -> xr.DataArray: """ if "node_lat" not in self._ds: _set_desired_longitude_range(self._ds) - _populate_lonlat_coord(self) + _populate_node_latlon(self) return self._ds["node_lat"] @property @@ -575,7 +575,7 @@ def node_x(self) -> xr.DataArray: Dimensions: ``(n_node, )`` """ if "node_x" not in self._ds: - _populate_cartesian_xyz_coord(self) + _populate_node_xyz(self) return self._ds["node_x"] @@ -586,7 +586,7 @@ def node_y(self) -> xr.DataArray: Dimensions: ``(n_node, )`` """ if "node_y" not in self._ds: - _populate_cartesian_xyz_coord(self) + _populate_node_xyz(self) return self._ds["node_y"] @property @@ -596,7 +596,7 @@ def node_z(self) -> xr.DataArray: Dimensions: ``(n_node, )`` """ if "node_z" not in self._ds: - _populate_cartesian_xyz_coord(self) + _populate_node_xyz(self) return self._ds["node_z"] @property diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index 7542fdbbb..0930fa6dc 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -1,352 +1,353 @@ -import numpy as np -from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE -from uxarray.grid.intersections import gca_constLat_intersection -from uxarray.grid.coordinates import node_xyz_to_lonlat_rad -import pandas as pd - -DUMMY_EDGE_VALUE = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE] - - -def _get_zonal_faces_weight_at_constLat( - faces_edges_cart, - latitude_cart, - face_latlon_bound, - is_directed=False, - is_latlonface=False, - is_face_GCA_list=None, -): - """Utilize the sweep line algorithm to calculate the weight of each face at - a constant latitude. - - Parameters - ---------- - face_edges_cart : np.ndarray - A list of face polygon represented by edges in Cartesian coordinates. Shape: (n_faces, n_edges, 2, 3) - - latitude_cart : float - The latitude in Cartesian coordinates (The normalized z coordinate) - - face_latlon_bound : np.ndarray - The list of latitude and longitude bounds of faces. Shape: (n_faces,2, 2), - [...,[lat_min, lat_max], [lon_min, lon_max],...] - - is_directed : bool, optional (default=False) - If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected, - and we will always assume the small circle (The one less than 180 degree) side is the GCA. - - is_latlonface : bool, optional, default=False - A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means - That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. - Default is False. - - is_face_GCA_list : np.ndarray, optional (default=None) - A list of boolean values that indicates if the edge in that face is a GCA. Shape: (n_faces,n_edges). - True means edge face is a GCA. - False mean this edge is a constant latitude. - If None, all edges are considered as GCA. This attribute will overwrite the is_latlonface attribute. - - Returns - ------- - weights_df : pandas.DataFrame - A DataFrame with the calculated weights of each face. The DataFrame has two columns: - - 'face_index': The index of the face (integer). - - 'weight': The calculated weight of the face in radian (float). - The DataFrame is indexed by the face indices, providing a mapping from each face to its corresponding weight. - """ - intervals_list = [] - - # Iterate through all faces and their edges - for face_index, face_edges in enumerate(faces_edges_cart): - if is_face_GCA_list is not None: - is_GCA_list = is_face_GCA_list[face_index] - else: - is_GCA_list = None - face_interval_df = _get_zonal_face_interval( - face_edges, - latitude_cart, - face_latlon_bound[face_index], - is_directed=is_directed, - is_latlonface=is_latlonface, - is_GCA_list=is_GCA_list, - ) - for _, row in face_interval_df.iterrows(): - intervals_list.append( - {"start": row["start"], "end": row["end"], "face_index": face_index} - ) - - intervals_df = pd.DataFrame(intervals_list) - overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) - - # Calculate weights for each face - weights = { - face_index: overlap_contributions.get(face_index, 0.0) / total_length - for face_index in range(len(faces_edges_cart)) - } - - # Convert weights to DataFrame - weights_df = pd.DataFrame(list(weights.items()), columns=["face_index", "weight"]) - return weights_df - - -def _is_edge_gca(is_GCA_list, is_latlonface, edges_z): - """Determine if each edge is a Great Circle Arc (GCA) or a constant - latitude line in a vectorized manner. - - Parameters: - ---------- - is_GCA_list : np.ndarray or None - An array indicating whether each edge is a GCA (True) or a constant latitude line (False). - Shape: (n_edges). If None, edge types are determined based on `is_latlonface` and the z-coordinates. - is_latlonface : bool - Flag indicating if all edges should be considered as lat-lon faces, which implies all edges - are either constant latitude or longitude lines. - edges_z : np.ndarray - Array containing the z-coordinates for each vertex of the edges. This is used to determine - whether edges are on the equator or if they are aligned in latitude when `is_GCA_list` is None. - Shape should be (n_edges, 2). - - Returns: - ------- - np.ndarray - A boolean array where each element indicates whether the corresponding edge is considered a GCA. - True for GCA, False for constant latitude line. - """ - if is_GCA_list is not None: - return is_GCA_list - if is_latlonface: - return ~np.isclose(edges_z[:, 0], edges_z[:, 1], atol=ERROR_TOLERANCE) - return ~( - np.isclose(edges_z[:, 0], 0, atol=ERROR_TOLERANCE) - & np.isclose(edges_z[:, 1], 0, atol=ERROR_TOLERANCE) - ) - - -def _get_faces_constLat_intersection_info( - face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed -): - """Processes each edge of a face polygon in a vectorized manner to - determine overlaps and calculate the intersections for a given latitude and - the faces. - - Parameters: - ---------- - face_edges_cart : np.ndarray - A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3). - latitude_cart : float - The latitude in Cartesian coordinates to which intersections or overlaps are calculated. - is_GCA_list : np.ndarray or None - An array indicating whether each edge is a GCA (True) or a constant latitude line (False). - Shape: (n_edges). If None, the function will determine edge types based on `is_latlonface`. - is_latlonface : bool - Flag indicating if all faces are considered as lat-lon faces, meaning all edges are either - constant latitude or longitude lines. This parameter overwrites the `is_GCA_list` if set to True. - is_directed : bool - Flag indicating if the GCA should be considered as directed (from v0 to v1). If False, - the smaller circle (less than 180 degrees) side of the GCA is used. - - Returns: - ------- - 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. - - overlap_flag (bool): A boolean indicating if any overlap with the latitude was detected. - - overlap_edge (np.ndarray or None): The edge that overlaps with the latitude, if any; otherwise, None. - """ - valid_edges_mask = ~(np.any(face_edges_cart == DUMMY_EDGE_VALUE, axis=(1, 2))) - - # Apply mask to filter out dummy edges - valid_edges = face_edges_cart[valid_edges_mask] - - # Extract Z coordinates for edge determination - edges_z = valid_edges[:, :, 2] - - # Determine if each edge is GCA or constant latitude - is_GCA = _is_edge_gca(is_GCA_list, is_latlonface, edges_z) - - # Check overlap with latitude - overlaps_with_latitude = np.all( - np.isclose(edges_z, latitude_cart, atol=ERROR_TOLERANCE), axis=1 - ) - overlap_flag = np.any(overlaps_with_latitude & ~is_GCA) - - # Identify overlap edges if needed - intersections_pts_list_cart = [] - if overlap_flag: - overlap_index = np.where(overlaps_with_latitude & ~is_GCA)[0][0] - intersections_pts_list_cart.extend(valid_edges[overlap_index]) - else: - # Calculate intersections (assuming a batch-capable intersection function) - for idx, edge in enumerate(valid_edges): - if is_GCA[idx]: - intersections = gca_constLat_intersection( - edge, latitude_cart, is_directed=is_directed - ) - if intersections.size == 0: - continue - elif intersections.shape[0] == 2: - intersections_pts_list_cart.extend(intersections) - else: - intersections_pts_list_cart.append(intersections[0]) - - # Handle unique intersections and check for convex or concave cases - unique_intersections = np.unique(intersections_pts_list_cart, axis=0) - if len(unique_intersections) == 2: - unique_intersection_lonlat = np.array( - [node_xyz_to_lonlat_rad(pt.tolist()) for pt in unique_intersections] - ) - 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) != 0: - raise ValueError( - "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly" - ) - elif len(unique_intersections) == 0: - raise ValueError( - "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" - ) - - -def _get_zonal_face_interval( - face_edges_cart, - latitude_cart, - face_latlon_bound, - is_directed=False, - is_latlonface=False, - is_GCA_list=None, -): - """Processes a face polygon represented by edges in Cartesian coordinates - to find intervals where the face intersects with a given latitude. This - function handles directed and undirected Great Circle Arcs (GCAs) and edges - at constant latitude. - - Requires the face edges to be sorted in counter-clockwise order, and the span of the - face in longitude should be less than pi. Also, all arcs/edges length should be within pi. - - Users can specify which edges are GCAs and which are constant latitude using `is_GCA_list`. - However, edges on the equator are always treated as constant latitude edges regardless of - `is_GCA_list`. - - Parameters - ---------- - face_edges_cart : np.ndarray - A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3) - latitude_cart : float - The latitude in cartesian, the normalized Z coordinates. - face_latlon_bound : np.ndarray - The latitude and longitude bounds of the face. Shape: (2, 2), [[lat_min, lat_max], [lon_min, lon_max]] - is_directed : bool, optional - If True, the GCA is considered to be directed (from v0 to v1). If False, the GCA is undirected, - and the smaller circle (less than 180 degrees) side of the GCA is used. Default is False. - is_latlonface : bool, optional, default=False - A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means - That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. - Default is False. This attribute will overwrite the is_latlonface attribute. - is_GCA_list : np.ndarray, optional, default=False - An array indicating if each edge is a GCA (True) or a constant latitude (False). Shape: (n_edges,). - If None, all edges are considered as GCAs. Default is None. - - Returns - ------- - Intervals_df : pd.DataFrame - A DataFrame containing the intervals stored as pandas Intervals for the face. - The columns of the DataFrame are: ['start', 'end'] - """ - face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1] - - # Call the vectorized function to process all edges - ( - unique_intersections, - pt_lon_min, - pt_lon_max, - ) = _get_faces_constLat_intersection_info( - face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed - ) - - # Convert intersection points to longitude-latitude - longitudes = np.array( - [node_xyz_to_lonlat_rad(pt.tolist())[0] for pt in unique_intersections] - ) - - # Handle special wrap-around cases by checking the face bounds - if face_lon_bound_left >= face_lon_bound_right: - if not ( - (pt_lon_max >= np.pi and pt_lon_min >= np.pi) - or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi) - ): - # They're at different sides of the 0-lon, adding wrap-around points - longitudes = np.append(longitudes, [0.0, 2 * np.pi]) - - # Ensure longitudes are sorted - longitudes.sort() - - # Split the sorted longitudes into start and end points of intervals - starts = longitudes[::2] # Start points - ends = longitudes[1::2] # End points - - # Create the intervals DataFrame - Intervals_df = pd.DataFrame({"start": starts, "end": ends}) - # For consistency, sort the intervals by start - interval_df_sorted = Intervals_df.sort_values(by="start").reset_index(drop=True) - return interval_df_sorted - - -def _process_overlapped_intervals(intervals_df): - """Process the overlapped intervals using the sweep line algorithm, - considering multiple intervals per face. - - Parameters - ---------- - intervals_df : pd.DataFrame - The DataFrame that contains the intervals and the corresponding face index. - The columns of the DataFrame are: ['start', 'end', 'face_index'] - - Returns - ------- - dict - A dictionary where keys are face indices and values are their respective contributions to the total length. - float - The total length of all intervals considering their overlaps. - - Example - ------- - >>> intervals_data = [ - ... {'start': 0.0, 'end': 100.0, 'face_index': 0}, - ... {'start': 50.0, 'end': 150.0, 'face_index': 1}, - ... {'start': 140.0, 'end': 150.0, 'face_index': 2} - ... ] - >>> intervals_df = pd.DataFrame(intervals_data) - >>> overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) - >>> print(overlap_contributions) - >>> print(total_length) - """ - events = [] - for idx, row in intervals_df.iterrows(): - events.append((row["start"], "start", row["face_index"])) - events.append((row["end"], "end", row["face_index"])) - - events.sort() # Sort by position and then by start/end - - active_faces = set() - last_position = None - total_length = 0.0 - overlap_contributions = {} - - for position, event_type, face_idx in events: - if last_position is not None and active_faces: - segment_length = position - last_position - segment_weight = segment_length / len(active_faces) if active_faces else 0 - for active_face in active_faces: - overlap_contributions[active_face] = ( - overlap_contributions.get(active_face, 0) + segment_weight - ) - total_length += segment_length - - if event_type == "start": - active_faces.add(face_idx) - elif event_type == "end": - active_faces.remove(face_idx) - - last_position = position - - return overlap_contributions, total_length +# import numpy as np +# from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE +# from uxarray.grid.intersections import gca_constLat_intersection +# +# # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad +# import pandas as pd +# +# DUMMY_EDGE_VALUE = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE] +# +# +# def _get_zonal_faces_weight_at_constLat( +# faces_edges_cart, +# latitude_cart, +# face_latlon_bound, +# is_directed=False, +# is_latlonface=False, +# is_face_GCA_list=None, +# ): +# """Utilize the sweep line algorithm to calculate the weight of each face at +# a constant latitude. +# +# Parameters +# ---------- +# face_edges_cart : np.ndarray +# A list of face polygon represented by edges in Cartesian coordinates. Shape: (n_faces, n_edges, 2, 3) +# +# latitude_cart : float +# The latitude in Cartesian coordinates (The normalized z coordinate) +# +# face_latlon_bound : np.ndarray +# The list of latitude and longitude bounds of faces. Shape: (n_faces,2, 2), +# [...,[lat_min, lat_max], [lon_min, lon_max],...] +# +# is_directed : bool, optional (default=False) +# If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected, +# and we will always assume the small circle (The one less than 180 degree) side is the GCA. +# +# is_latlonface : bool, optional, default=False +# A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means +# That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. +# Default is False. +# +# is_face_GCA_list : np.ndarray, optional (default=None) +# A list of boolean values that indicates if the edge in that face is a GCA. Shape: (n_faces,n_edges). +# True means edge face is a GCA. +# False mean this edge is a constant latitude. +# If None, all edges are considered as GCA. This attribute will overwrite the is_latlonface attribute. +# +# Returns +# ------- +# weights_df : pandas.DataFrame +# A DataFrame with the calculated weights of each face. The DataFrame has two columns: +# - 'face_index': The index of the face (integer). +# - 'weight': The calculated weight of the face in radian (float). +# The DataFrame is indexed by the face indices, providing a mapping from each face to its corresponding weight. +# """ +# intervals_list = [] +# +# # Iterate through all faces and their edges +# for face_index, face_edges in enumerate(faces_edges_cart): +# if is_face_GCA_list is not None: +# is_GCA_list = is_face_GCA_list[face_index] +# else: +# is_GCA_list = None +# face_interval_df = _get_zonal_face_interval( +# face_edges, +# latitude_cart, +# face_latlon_bound[face_index], +# is_directed=is_directed, +# is_latlonface=is_latlonface, +# is_GCA_list=is_GCA_list, +# ) +# for _, row in face_interval_df.iterrows(): +# intervals_list.append( +# {"start": row["start"], "end": row["end"], "face_index": face_index} +# ) +# +# intervals_df = pd.DataFrame(intervals_list) +# overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) +# +# # Calculate weights for each face +# weights = { +# face_index: overlap_contributions.get(face_index, 0.0) / total_length +# for face_index in range(len(faces_edges_cart)) +# } +# +# # Convert weights to DataFrame +# weights_df = pd.DataFrame(list(weights.items()), columns=["face_index", "weight"]) +# return weights_df +# +# +# def _is_edge_gca(is_GCA_list, is_latlonface, edges_z): +# """Determine if each edge is a Great Circle Arc (GCA) or a constant +# latitude line in a vectorized manner. +# +# Parameters: +# ---------- +# is_GCA_list : np.ndarray or None +# An array indicating whether each edge is a GCA (True) or a constant latitude line (False). +# Shape: (n_edges). If None, edge types are determined based on `is_latlonface` and the z-coordinates. +# is_latlonface : bool +# Flag indicating if all edges should be considered as lat-lon faces, which implies all edges +# are either constant latitude or longitude lines. +# edges_z : np.ndarray +# Array containing the z-coordinates for each vertex of the edges. This is used to determine +# whether edges are on the equator or if they are aligned in latitude when `is_GCA_list` is None. +# Shape should be (n_edges, 2). +# +# Returns: +# ------- +# np.ndarray +# A boolean array where each element indicates whether the corresponding edge is considered a GCA. +# True for GCA, False for constant latitude line. +# """ +# if is_GCA_list is not None: +# return is_GCA_list +# if is_latlonface: +# return ~np.isclose(edges_z[:, 0], edges_z[:, 1], atol=ERROR_TOLERANCE) +# return ~( +# np.isclose(edges_z[:, 0], 0, atol=ERROR_TOLERANCE) +# & np.isclose(edges_z[:, 1], 0, atol=ERROR_TOLERANCE) +# ) +# +# +# def _get_faces_constLat_intersection_info( +# face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed +# ): +# """Processes each edge of a face polygon in a vectorized manner to +# determine overlaps and calculate the intersections for a given latitude and +# the faces. +# +# Parameters: +# ---------- +# face_edges_cart : np.ndarray +# A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3). +# latitude_cart : float +# The latitude in Cartesian coordinates to which intersections or overlaps are calculated. +# is_GCA_list : np.ndarray or None +# An array indicating whether each edge is a GCA (True) or a constant latitude line (False). +# Shape: (n_edges). If None, the function will determine edge types based on `is_latlonface`. +# is_latlonface : bool +# Flag indicating if all faces are considered as lat-lon faces, meaning all edges are either +# constant latitude or longitude lines. This parameter overwrites the `is_GCA_list` if set to True. +# is_directed : bool +# Flag indicating if the GCA should be considered as directed (from v0 to v1). If False, +# the smaller circle (less than 180 degrees) side of the GCA is used. +# +# Returns: +# ------- +# 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. +# - overlap_flag (bool): A boolean indicating if any overlap with the latitude was detected. +# - overlap_edge (np.ndarray or None): The edge that overlaps with the latitude, if any; otherwise, None. +# """ +# valid_edges_mask = ~(np.any(face_edges_cart == DUMMY_EDGE_VALUE, axis=(1, 2))) +# +# # Apply mask to filter out dummy edges +# valid_edges = face_edges_cart[valid_edges_mask] +# +# # Extract Z coordinates for edge determination +# edges_z = valid_edges[:, :, 2] +# +# # Determine if each edge is GCA or constant latitude +# is_GCA = _is_edge_gca(is_GCA_list, is_latlonface, edges_z) +# +# # Check overlap with latitude +# overlaps_with_latitude = np.all( +# np.isclose(edges_z, latitude_cart, atol=ERROR_TOLERANCE), axis=1 +# ) +# overlap_flag = np.any(overlaps_with_latitude & ~is_GCA) +# +# # Identify overlap edges if needed +# intersections_pts_list_cart = [] +# if overlap_flag: +# overlap_index = np.where(overlaps_with_latitude & ~is_GCA)[0][0] +# intersections_pts_list_cart.extend(valid_edges[overlap_index]) +# else: +# # Calculate intersections (assuming a batch-capable intersection function) +# for idx, edge in enumerate(valid_edges): +# if is_GCA[idx]: +# intersections = gca_constLat_intersection( +# edge, latitude_cart, is_directed=is_directed +# ) +# if intersections.size == 0: +# continue +# elif intersections.shape[0] == 2: +# intersections_pts_list_cart.extend(intersections) +# else: +# intersections_pts_list_cart.append(intersections[0]) +# +# # Handle unique intersections and check for convex or concave cases +# unique_intersections = np.unique(intersections_pts_list_cart, axis=0) +# if len(unique_intersections) == 2: +# unique_intersection_lonlat = np.array( +# [node_xyz_to_lonlat_rad(pt.tolist()) for pt in unique_intersections] +# ) +# 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) != 0: +# raise ValueError( +# "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly" +# ) +# elif len(unique_intersections) == 0: +# raise ValueError( +# "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" +# ) +# +# +# def _get_zonal_face_interval( +# face_edges_cart, +# latitude_cart, +# face_latlon_bound, +# is_directed=False, +# is_latlonface=False, +# is_GCA_list=None, +# ): +# """Processes a face polygon represented by edges in Cartesian coordinates +# to find intervals where the face intersects with a given latitude. This +# function handles directed and undirected Great Circle Arcs (GCAs) and edges +# at constant latitude. +# +# Requires the face edges to be sorted in counter-clockwise order, and the span of the +# face in longitude should be less than pi. Also, all arcs/edges length should be within pi. +# +# Users can specify which edges are GCAs and which are constant latitude using `is_GCA_list`. +# However, edges on the equator are always treated as constant latitude edges regardless of +# `is_GCA_list`. +# +# Parameters +# ---------- +# face_edges_cart : np.ndarray +# A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3) +# latitude_cart : float +# The latitude in cartesian, the normalized Z coordinates. +# face_latlon_bound : np.ndarray +# The latitude and longitude bounds of the face. Shape: (2, 2), [[lat_min, lat_max], [lon_min, lon_max]] +# is_directed : bool, optional +# If True, the GCA is considered to be directed (from v0 to v1). If False, the GCA is undirected, +# and the smaller circle (less than 180 degrees) side of the GCA is used. Default is False. +# is_latlonface : bool, optional, default=False +# A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means +# That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. +# Default is False. This attribute will overwrite the is_latlonface attribute. +# is_GCA_list : np.ndarray, optional, default=False +# An array indicating if each edge is a GCA (True) or a constant latitude (False). Shape: (n_edges,). +# If None, all edges are considered as GCAs. Default is None. +# +# Returns +# ------- +# Intervals_df : pd.DataFrame +# A DataFrame containing the intervals stored as pandas Intervals for the face. +# The columns of the DataFrame are: ['start', 'end'] +# """ +# face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1] +# +# # Call the vectorized function to process all edges +# ( +# unique_intersections, +# pt_lon_min, +# pt_lon_max, +# ) = _get_faces_constLat_intersection_info( +# face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed +# ) +# +# # Convert intersection points to longitude-latitude +# longitudes = np.array( +# [node_xyz_to_lonlat_rad(pt.tolist())[0] for pt in unique_intersections] +# ) +# +# # Handle special wrap-around cases by checking the face bounds +# if face_lon_bound_left >= face_lon_bound_right: +# if not ( +# (pt_lon_max >= np.pi and pt_lon_min >= np.pi) +# or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi) +# ): +# # They're at different sides of the 0-lon, adding wrap-around points +# longitudes = np.append(longitudes, [0.0, 2 * np.pi]) +# +# # Ensure longitudes are sorted +# longitudes.sort() +# +# # Split the sorted longitudes into start and end points of intervals +# starts = longitudes[::2] # Start points +# ends = longitudes[1::2] # End points +# +# # Create the intervals DataFrame +# Intervals_df = pd.DataFrame({"start": starts, "end": ends}) +# # For consistency, sort the intervals by start +# interval_df_sorted = Intervals_df.sort_values(by="start").reset_index(drop=True) +# return interval_df_sorted +# +# +# def _process_overlapped_intervals(intervals_df): +# """Process the overlapped intervals using the sweep line algorithm, +# considering multiple intervals per face. +# +# Parameters +# ---------- +# intervals_df : pd.DataFrame +# The DataFrame that contains the intervals and the corresponding face index. +# The columns of the DataFrame are: ['start', 'end', 'face_index'] +# +# Returns +# ------- +# dict +# A dictionary where keys are face indices and values are their respective contributions to the total length. +# float +# The total length of all intervals considering their overlaps. +# +# Example +# ------- +# >>> intervals_data = [ +# ... {'start': 0.0, 'end': 100.0, 'face_index': 0}, +# ... {'start': 50.0, 'end': 150.0, 'face_index': 1}, +# ... {'start': 140.0, 'end': 150.0, 'face_index': 2} +# ... ] +# >>> intervals_df = pd.DataFrame(intervals_data) +# >>> overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) +# >>> print(overlap_contributions) +# >>> print(total_length) +# """ +# events = [] +# for idx, row in intervals_df.iterrows(): +# events.append((row["start"], "start", row["face_index"])) +# events.append((row["end"], "end", row["face_index"])) +# +# events.sort() # Sort by position and then by start/end +# +# active_faces = set() +# last_position = None +# total_length = 0.0 +# overlap_contributions = {} +# +# for position, event_type, face_idx in events: +# if last_position is not None and active_faces: +# segment_length = position - last_position +# segment_weight = segment_length / len(active_faces) if active_faces else 0 +# for active_face in active_faces: +# overlap_contributions[active_face] = ( +# overlap_contributions.get(active_face, 0) + segment_weight +# ) +# total_length += segment_length +# +# if event_type == "start": +# active_faces.add(face_idx) +# elif event_type == "end": +# active_faces.remove(face_idx) +# +# last_position = position +# +# return overlap_contributions, total_length diff --git a/uxarray/grid/veccoordinates.py b/uxarray/grid/veccoordinates.py deleted file mode 100644 index b40d8d644..000000000 --- a/uxarray/grid/veccoordinates.py +++ /dev/null @@ -1,131 +0,0 @@ -import xarray as xr -import numpy as np - -from numba import config - -from uxarray.constants import ENABLE_JIT, ERROR_TOLERANCE -from uxarray.conventions import ugrid - -config.DISABLE_JIT = not ENABLE_JIT - - -# ====================================================================================================================== -# Conversion Functions -# ====================================================================================================================== -def _lonlat_rad_to_xyz( - lon: np.ndarray, lat: np.ndarray -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Docstring TODO.""" - 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: np.ndarray, y: np.ndarray, z: np.ndarray -) -> tuple[np.ndarray, np.ndarray]: - """Docstring TODO.""" - denom = np.abs(x * x + y * y + z * z) - - x /= denom - y /= denom - z /= denom - - lon = np.empty_like(x) - lat = np.empty_like(x) - - # is there a way to vectorize this loop? - for i in range(len(x)): - if np.absolute(z[i]) < (1.0 - ERROR_TOLERANCE): - lon[i] = np.arctan2(y[i], x[i]) - lat[i] = np.arcsin(z[i]) - - if lon[i] < 0.0: - lon[i] += 2.0 * np.pi - elif z[i] > 0.0: - lon[i] = 0.0 - lat[i] = 0.5 * np.pi - else: - lon[i] = 0.0 - lat[i] = -0.5 * np.pi - - return lon, lat - - -# don't use numba here unless we figure out ord=2 -def _normalize_xyz( - x: np.ndarray, y: np.ndarray, z: np.ndarray -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Docstring TODO.""" - x_norm = x / np.linalg.norm(x, ord=2) - y_norm = y / np.linalg.norm(y, ord=2) - z_norm = z / np.linalg.norm(z, ord=2) - return x_norm, y_norm, z_norm - - -# ====================================================================================================================== -# Population Functions -# ====================================================================================================================== -def _populate_node_latlon(grid) -> None: - """Docstring TODO.""" - x_norm, y_norm, z_norm = _normalize_xyz( - grid.node_x.values, grid.node_y.values, grid.node_z.values - ) - lon_rad, lat_rad = _xyz_to_lonlat_rad(x_norm, y_norm, z_norm) - - # Convert to degrees - lon = np.rad2deg(lon_rad) - lat = np.rad2deg(lat_rad) - - grid._ds["node_lon"] = xr.DataArray( - data=lon, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LON_ATTRS - ) - grid._ds["node_lat"] = xr.DataArray( - data=lat, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LAT_ATTRS - ) - - -def _populate_edge_latlon(grid) -> None: - pass - - -def _populate_face_latlon(grid) -> None: - pass - - -def _populate_node_xyz(grid) -> None: - """Docstring TODO.""" - - node_lon_rad = np.deg2rad(grid.node_lon.values) - node_lat_rad = np.deg2rad(grid.node_lat.values) - x, y, z = _lonlat_rad_to_xyz(node_lon_rad, node_lat_rad) - - grid._ds["node_x"] = xr.DataArray( - data=x, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_X_ATTRS - ) - grid._ds["node_y"] = xr.DataArray( - data=y, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_Y_ATTRS - ) - grid._ds["node_z"] = xr.DataArray( - data=z, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_Z_ATTRS - ) - - -def _populate_edge_xyz(grid) -> None: - pass - - -def _populate_face_xyz(grid) -> None: - pass - - -# ====================================================================================================================== -# Build Functions -# ====================================================================================================================== -def _build_face_centroids(): - pass - - -def _build_edge_centroids(): - pass diff --git a/uxarray/io/_exodus.py b/uxarray/io/_exodus.py index 205b68d36..8c91913fa 100644 --- a/uxarray/io/_exodus.py +++ b/uxarray/io/_exodus.py @@ -5,7 +5,8 @@ from uxarray.grid.connectivity import _replace_fill_values from uxarray.constants import INT_DTYPE, INT_FILL_VALUE -from uxarray.grid.coordinates import _get_lonlat_from_xyz, _get_xyz_from_lonlat + +from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_deg from uxarray.conventions import ugrid @@ -100,7 +101,7 @@ def _read_exodus(ext_ds): ) # populate lon/lat coordinates - lon, lat = _get_lonlat_from_xyz( + lon, lat = _xyz_to_lonlat_deg( ds["node_x"].values, ds["node_y"].values, ds["node_z"].values ) @@ -181,7 +182,7 @@ def _encode_exodus(ds, outfile=None): # Set the dim to 3 as we will always have x/y/z for cartesian grid # Note: Don't get orig dimension from Mesh2 attribute topology dimension if "node_x" not in ds: - x, y, z = _get_xyz_from_lonlat(ds["node_lon"].values, ds["node_lat"].values) + x, y, z = _lonlat_rad_to_xyz(ds["node_lon"].values, ds["node_lat"].values) c_data = xr.DataArray([x, y, z]) else: c_data = xr.DataArray( From 95ccb04b59bdcc77c1314cda374a8cfb1c6f4131 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 8 May 2024 13:50:49 -0500 Subject: [PATCH 03/14] update existing code and tests to use new functions --- test/test_arcs.py | 35 +- test/test_computing.py | 10 +- test/test_intersections.py | 76 ++-- uxarray/grid/arcs.py | 620 ++++++++++++++++--------------- uxarray/grid/area.py | 26 +- uxarray/grid/coordinates.py | 9 +- uxarray/grid/integrate.py | 714 ++++++++++++++++++------------------ 7 files changed, 741 insertions(+), 749 deletions(-) diff --git a/test/test_arcs.py b/test/test_arcs.py index c2e24f5aa..c9859e724 100644 --- a/test/test_arcs.py +++ b/test/test_arcs.py @@ -9,7 +9,7 @@ import uxarray as ux -from uxarray.grid.coordinates import node_lonlat_rad_to_xyz +from uxarray.grid.coordinates import _lonlat_rad_to_xyz from uxarray.grid.arcs import point_within_gca try: @@ -31,43 +31,40 @@ class TestIntersectionPoint(TestCase): def test_pt_within_gcr(self): # The GCR that's eexactly 180 degrees will have Value Error raised gcr_180degree_cart = [ - ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0]), - ux.grid.coordinates.node_lonlat_rad_to_xyz([np.pi, 0.0]) + _lonlat_rad_to_xyz(0.0, 0.0), + _lonlat_rad_to_xyz(np.pi, 0.0) ] - pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0]) + pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) with self.assertRaises(ValueError): point_within_gca(pt_same_lon_in, gcr_180degree_cart) gcr_180degree_cart = [ - ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, np.pi / 2.0]), - ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, -np.pi / 2.0]) + _lonlat_rad_to_xyz(0.0, np.pi / 2.0), + _lonlat_rad_to_xyz(0.0, -np.pi / 2.0) ] - pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0]) + pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) with self.assertRaises(ValueError): point_within_gca(pt_same_lon_in, gcr_180degree_cart) # Test when the point and the GCA all have the same longitude gcr_same_lon_cart = [ - ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 1.5]), - ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, -1.5]) + _lonlat_rad_to_xyz(0.0, 1.5), + _lonlat_rad_to_xyz(0.0, -1.5) ] - pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0]) + pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart)) - pt_same_lon_out = ux.grid.coordinates.node_lonlat_rad_to_xyz( - [0.0, 1.500000000000001]) + pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001) res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart) self.assertFalse(res) - pt_same_lon_out_2 = ux.grid.coordinates.node_lonlat_rad_to_xyz( - [0.1, 1.0]) + pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0) res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart) self.assertFalse(res) # And if we increase the digital place by one, it should be true again - pt_same_lon_out_add_one_place = ux.grid.coordinates.node_lonlat_rad_to_xyz( - [0.0, 1.5000000000000001]) + pt_same_lon_out_add_one_place = _lonlat_rad_to_xyz(0.0, 1.5000000000000001) res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart) self.assertTrue(res) @@ -117,10 +114,10 @@ def test_pt_within_gcr_antimeridian(self): # The first case should not work and the second should work v1_rad = [0.1, 0.0] v2_rad = [2 * np.pi - 0.1, 0.0] - v1_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz(v1_rad) - v2_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz(v2_rad) + v1_cart = _lonlat_rad_to_xyz(v1_rad[0], v1_rad[1]) + v2_cart = _lonlat_rad_to_xyz(v2_rad[0], v1_rad[1]) gcr_cart = np.array([v1_cart, v2_cart]) - pt_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.01, 0.0]) + pt_cart = _lonlat_rad_to_xyz(0.01, 0.0) with self.assertRaises(ValueError): point_within_gca(pt_cart, gcr_cart, is_directed=True) gcr_car_flipped = np.array([v2_cart, v1_cart]) diff --git a/test/test_computing.py b/test/test_computing.py index 4e9cde81b..a83d3b76c 100644 --- a/test/test_computing.py +++ b/test/test_computing.py @@ -1,7 +1,7 @@ from unittest import TestCase import numpy.testing as nt import numpy as np -from uxarray.grid.coordinates import normalize_in_place +from uxarray.grid.coordinates import _normalize_xyz import uxarray.utils.computing as ac_utils from uxarray.constants import ERROR_TOLERANCE @@ -12,8 +12,8 @@ class TestCrossProduct(TestCase): one.""" def test_cross_fma(self): - v1 = np.array(normalize_in_place([1.0, 2.0, 3.0])) - v2 = np.array(normalize_in_place([4.0, 5.0, 6.0])) + v1 = np.array(_normalize_xyz(*[1.0, 2.0, 3.0])) + v2 = np.array(_normalize_xyz(*[4.0, 5.0, 6.0])) np_cross = np.cross(v1, v2) fma_cross = ac_utils.cross_fma(v1, v2) @@ -26,8 +26,8 @@ class TestDotProduct(TestCase): one.""" def test_dot_fma(self): - v1 = np.array(normalize_in_place([1.0, 0.0, 0.0]), dtype=np.float64) - v2 = np.array(normalize_in_place([1.0, 0.0, 0.0]), dtype=np.float64) + v1 = np.array(_normalize_xyz(*[1.0, 0.0, 0.0]), dtype=np.float64) + v2 = np.array(_normalize_xyz(*[1.0, 0.0, 0.0]), dtype=np.float64) np_dot = np.dot(v1, v2) fma_dot = ac_utils.dot_fma(v1, v2) diff --git a/test/test_intersections.py b/test/test_intersections.py index 6291b93b0..b134d9029 100644 --- a/test/test_intersections.py +++ b/test/test_intersections.py @@ -17,16 +17,16 @@ def test_get_GCA_GCA_intersections_antimeridian(self): - GCA1 = node_lonlat_rad_to_xyz([np.deg2rad(170.0), np.deg2rad(89.99)]) + GCA1 = _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(89.99)) GCR1_cart = np.array([ - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(89.99)]), - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(10.0)]) + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(89.99)), + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(10.0)) ]) GCR2_cart = np.array([ - node_lonlat_rad_to_xyz([np.deg2rad(70.0), 0.0]), - node_lonlat_rad_to_xyz([np.deg2rad(179.0), 0.0]) + _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0), + _lonlat_rad_to_xyz(np.deg2rad(179.0), 0.0) ]) res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) @@ -34,14 +34,14 @@ def test_get_GCA_GCA_intersections_antimeridian(self): self.assertTrue(np.array_equal(res_cart, np.array([]))) GCR1_cart = np.array([ - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(89.0)]), - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(-10.0)]) + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(89.0)), + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(-10.0)) ]) GCR2_cart = np.array([ - node_lonlat_rad_to_xyz([np.deg2rad(70.0), 0.0]), - node_lonlat_rad_to_xyz([np.deg2rad(175.0), 0.0]) + _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0), + _lonlat_rad_to_xyz(np.deg2rad(175.0), 0.0) ]) res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) @@ -51,7 +51,7 @@ def test_get_GCA_GCA_intersections_antimeridian(self): np.allclose(np.linalg.norm(res_cart, axis=0), 1.0, atol=ERROR_TOLERANCE)) - res_lonlat_rad = node_xyz_to_lonlat_rad(res_cart.tolist()) + res_lonlat_rad = _xyz_to_lonlat_rad(res_cart[0], res_cart[1], res_cart[2]) # res_cart should be [170, 0] self.assertTrue( @@ -62,12 +62,12 @@ def test_get_GCA_GCA_intersections_antimeridian(self): def test_get_GCA_GCA_intersections_parallel(self): # Test the case where the two GCAs are parallel GCR1_cart = np.array([ - node_lonlat_rad_to_xyz([0.3 * np.pi, 0.0]), - node_lonlat_rad_to_xyz([0.5 * np.pi, 0.0]) + _lonlat_rad_to_xyz(0.3 * np.pi, 0.0), + _lonlat_rad_to_xyz(0.5 * np.pi, 0.0) ]) GCR2_cart = np.array([ - node_lonlat_rad_to_xyz([0.5 * np.pi, 0.0]), - node_lonlat_rad_to_xyz([-0.5 * np.pi - 0.01, 0.0]) + _lonlat_rad_to_xyz(0.5 * np.pi, 0.0), + _lonlat_rad_to_xyz(-0.5 * np.pi - 0.01, 0.0) ]) res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) self.assertTrue(np.array_equal(res_cart, np.array([]))) @@ -75,14 +75,14 @@ def test_get_GCA_GCA_intersections_parallel(self): def test_get_GCA_GCA_intersections_perpendicular(self): # Test the case where the two GCAs are perpendicular to each other GCR1_cart = np.array([ - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(0.0)]), - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(10.0)]) + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(0.0)), + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(10.0)) ]) GCR2_cart = np.array([ - node_lonlat_rad_to_xyz([0.5 * np.pi, 0.0]), - node_lonlat_rad_to_xyz([-0.5 * np.pi - 0.01, 0.0]) + _lonlat_rad_to_xyz(*[0.5 * np.pi, 0.0]), + _lonlat_rad_to_xyz(*[-0.5 * np.pi - 0.01, 0.0]) ]) res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) @@ -91,7 +91,7 @@ def test_get_GCA_GCA_intersections_perpendicular(self): np.allclose(np.linalg.norm(res_cart, axis=0), 1.0, atol=ERROR_TOLERANCE)) - res_lonlat_rad = node_xyz_to_lonlat_rad(res_cart.tolist()) + res_lonlat_rad = _xyz_to_lonlat_rad(*res_cart) self.assertTrue( np.allclose(res_lonlat_rad, np.array([np.deg2rad(170.0), @@ -102,14 +102,14 @@ class TestGCAconstLatIntersection(TestCase): def test_GCA_constLat_intersections_antimeridian(self): GCR1_cart = np.array([ - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(89.99)]), - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(10.0)]) + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(89.99)), + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(10.0)) ]) res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(60.0)), verbose=True) - res_lonlat_rad = node_xyz_to_lonlat_rad(res[0].tolist()) + res_lonlat_rad = _xyz_to_lonlat_rad(*(res[0].tolist())) self.assertTrue( np.allclose(res_lonlat_rad, np.array([np.deg2rad(170.0), @@ -117,10 +117,10 @@ def test_GCA_constLat_intersections_antimeridian(self): def test_GCA_constLat_intersections_empty(self): GCR1_cart = np.array([ - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(89.99)]), - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(10.0)]) + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(89.99)), + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(10.0)) ]) res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(-10.0)), verbose=False) @@ -128,10 +128,10 @@ def test_GCA_constLat_intersections_empty(self): def test_GCA_constLat_intersections_two_pts(self): GCR1_cart = np.array([ - node_lonlat_rad_to_xyz([np.deg2rad(10.0), - np.deg2rad(10)]), - node_lonlat_rad_to_xyz([np.deg2rad(170.0), - np.deg2rad(10.0)]) + _lonlat_rad_to_xyz(np.deg2rad(10.0), + np.deg2rad(10)), + _lonlat_rad_to_xyz(np.deg2rad(170.0), + np.deg2rad(10.0)) ]) max_lat = ux.grid.arcs.extreme_gca_latitude(GCR1_cart, 'max') diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index df4a9dc7f..82954ccab 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -1,312 +1,308 @@ -# import numpy as np -# -# # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place -# from uxarray.grid.coordinates import _xyz_to_lonlat_rad -# -# -# from uxarray.constants import ERROR_TOLERANCE -# -# -# def _to_list(obj): -# if not isinstance(obj, list): -# if isinstance(obj, np.ndarray): -# # Convert the NumPy array to a list using .tolist() -# obj = obj.tolist() -# else: -# # If not a list or NumPy array, return the object as-is -# obj = [obj] -# return obj -# -# -# def point_within_gca(pt, gca_cart, is_directed=False): -# """Check if a point lies on a given Great Circle Arc (GCA). The anti- -# meridian case is also considered. -# -# Parameters -# ---------- -# pt : numpy.ndarray (float) -# Cartesian coordinates of the point. -# gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr) -# Cartesian coordinates of the Great Circle Arc (GCR). -# is_directed : bool, optional, default = False -# If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected, -# and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False. -# For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1 -# -# Returns -# ------- -# bool -# True if the point lies between the two endpoints of the GCR, False otherwise. -# -# Raises -# ------ -# ValueError -# If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes. -# In such cases, consider breaking the GCR into two separate GCRs. -# -# ValueError -# If the input GCR spans more than 180 degrees (π radians). -# In such cases, consider breaking the GCR into two separate GCRs. -# -# Notes -# ----- -# The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and -# accounting for the anti-meridian case. -# -# The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude). -# In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation. -# -# The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations. -# -# Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. -# """ -# # Convert the cartesian coordinates to lonlat coordinates -# -# pt_lon, pt_lat = _xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) -# -# GCRv0_lon, GCRv0_lat = _xyz_to_lonlat_rad( -# gca_cart[0][0], gca_cart[0][1], gca_cart[0][2] -# ) -# GCRv1_lon, GCRv1_lat = _xyz_to_lonlat_rad( -# gca_cart[1][0], gca_cart[1][1], gca_cart[1][2] -# ) -# -# # Convert the list to np.float64 -# gca_cart[0] = np.array(gca_cart[0], dtype=np.float64) -# gca_cart[1] = np.array(gca_cart[1], dtype=np.float64) -# -# # First if the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes -# angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1]) -# if np.allclose(angle, np.pi, rtol=0, atol=ERROR_TOLERANCE): -# raise ValueError( -# "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. " -# "Consider breaking the Great Circle Arc" -# "into two Great Circle Arcs" -# ) -# -# if not np.allclose( -# np.dot(np.cross(gca_cart[0], gca_cart[1]), pt), 0, rtol=0, atol=ERROR_TOLERANCE -# ): -# return False -# -# if np.isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): -# # If the pt and the GCA are on the same longitude (the y coordinates are the same) -# if np.isclose(GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): -# # Now use the latitude to determine if the pt falls between the interval -# return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1]) -# else: -# # If the pt and the GCA are not on the same longitude when the GCA is a longnitude arc, then the pt is not on the GCA -# return False -# -# # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point -# if np.isclose( -# abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0, atol=ERROR_TOLERANCE -# ): -# # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0] -# if np.isclose(abs(pt_lonlat[1]), np.pi / 2, rtol=0, atol=ERROR_TOLERANCE): -# pt_lonlat[0] = GCRv0_lonlat[0] -# if not np.isclose( -# GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE -# ) and not np.isclose( -# GCRv1_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE -# ): -# return False -# else: -# # Determine the pole latitude and latitude extension -# if (GCRv0_lonlat[1] > 0 and GCRv1_lonlat[1] > 0) or ( -# GCRv0_lonlat[1] < 0 and GCRv1_lonlat[1] < 0 -# ): -# pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0 else -np.pi / 2 -# lat_extend = ( -# abs(np.pi / 2 - abs(GCRv0_lonlat[1])) -# + np.pi / 2 -# + abs(GCRv1_lonlat[1]) -# ) -# else: -# pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1]) -# lat_extend = ( -# abs(np.pi / 2 - abs(GCRv0_lonlat[1])) -# + np.pi / 2 -# + abs(GCRv1_lonlat[1]) -# ) -# -# if is_directed and lat_extend >= np.pi: -# raise ValueError( -# "The input GCA spans more than 180 degrees. Please divide it into two." -# ) -# -# return in_between(GCRv0_lonlat[1], pt_lonlat[1], pole_lat) or in_between( -# pole_lat, pt_lonlat[1], GCRv1_lonlat[1] -# ) -# -# if is_directed: -# # The anti-meridian case Sufficient condition: absolute difference between the longitudes of the two -# # vertices is greater than 180 degrees (π radians): abs(GCRv1_lon - GCRv0_lon) > π -# if abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]) > np.pi: -# # The necessary condition: the pt longitude is on the opposite side of the anti-meridian -# # Case 1: where 0 --> x0--> 180 -->x1 -->0 case is lager than the 180degrees (pi radians) -# if GCRv0_lonlat[0] <= np.pi <= GCRv1_lonlat[0]: -# raise ValueError( -# "The input Great Circle Arc span is larger than 180 degree, please break it into two." -# ) -# -# # The necessary condition: the pt longitude is on the opposite side of the anti-meridian -# # Case 2: The anti-meridian case where 180 -->x0 --> 0 lon --> x1 --> 180 lon -# elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0: -# return in_between( -# GCRv0_lonlat[0], pt_lonlat[0], 2 * np.pi -# ) or in_between(0, pt_lonlat[0], GCRv1_lonlat[0]) -# -# # The non-anti-meridian case. -# else: -# return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) -# else: -# # The undirected case -# # sort the longitude -# GCRv0_lonlat_min, GCRv1_lonlat_max = sorted([GCRv0_lonlat[0], GCRv1_lonlat[0]]) -# if np.pi > GCRv1_lonlat_max - GCRv0_lonlat_min >= 0.0: -# return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) -# else: -# return in_between(GCRv1_lonlat_max, pt_lonlat[0], 2 * np.pi) or in_between( -# 0.0, pt_lonlat[0], GCRv0_lonlat_min -# ) -# -# -# def in_between(p, q, r) -> bool: -# """Determines whether the number q is between p and r. -# -# Parameters -# ---------- -# p : float -# The lower bound. -# q : float -# The number to check. -# r : float -# The upper bound. -# -# Returns -# ------- -# bool -# True if q is between p and r, False otherwise. -# """ -# -# return p <= q <= r or r <= q <= p -# -# -# def _decide_pole_latitude(lat1, lat2): -# """Determine the pole latitude based on the latitudes of two points on a -# Great Circle Arc (GCA). -# -# This function calculates the combined latitude span from each point to its nearest pole -# and decides which pole (North or South) the smaller GCA will pass. This decision is crucial -# for handling GCAs that span exactly or more than 180 degrees in longitude, indicating -# the arc might pass through or close to one of the Earth's poles. -# -# Parameters -# ---------- -# lat1 : float -# Latitude of the first point in radians. Positive for the Northern Hemisphere, negative for the Southern. -# lat2 : float -# Latitude of the second point in radians. Positive for the Northern Hemisphere, negative for the Southern. -# -# Returns -# ------- -# float -# The latitude of the pole (np.pi/2 for the North Pole or -np.pi/2 for the South Pole) the GCA is closer to. -# -# Notes -# ----- -# The function assumes the input latitudes are valid (i.e., between -np.pi/2 and np.pi/2) and expressed in radians. -# The determination of which pole a GCA is closer to is based on the sum of the latitudinal spans from each point -# to its nearest pole, considering the shortest path on the sphere. -# """ -# # Calculate the total latitudinal span to the nearest poles -# lat_extend = abs(np.pi / 2 - abs(lat1)) + np.pi / 2 + abs(lat2) -# -# # Determine the closest pole based on the latitudinal span -# if lat_extend < np.pi: -# closest_pole = np.pi / 2 if lat1 > 0 else -np.pi / 2 -# else: -# closest_pole = -np.pi / 2 if lat1 > 0 else np.pi / 2 -# -# return closest_pole -# -# -# 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. -# -# Parameters -# ---------- -# gca_cart : numpy.ndarray -# An array containing two 3D vectors that define a great circle arc. -# -# extreme_type : str -# The type of extreme latitude to calculate. Must be either 'max' or 'min'. -# -# Returns -# ------- -# float -# The maximum or minimum latitude of the great circle arc in radians. -# -# Raises -# ------ -# ValueError -# If `extreme_type` is not 'max' or 'min'. -# """ -# extreme_type = extreme_type.lower() -# -# if extreme_type not in ("max", "min"): -# raise ValueError("extreme_type must be either 'max' or 'min'") -# -# n1, n2 = gca_cart -# dot_n1_n2 = np.dot(n1, n2) -# denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0) -# d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom -# -# d_a_max = ( -# np.clip(d_a_max, 0, 1) -# if np.isclose(d_a_max, [0, 1], atol=ERROR_TOLERANCE).any() -# else d_a_max -# ) -# _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2]) -# _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2]) -# -# if 0 < d_a_max < 1: -# node3 = (1 - d_a_max) * n1 + d_a_max * n2 -# # TODO: -# node3 = np.array(normalize_in_place(node3.tolist())) -# node3 = np.array(normalize) -# d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1)) -# -# return ( -# max(d_lat_rad, lat_n1, lat_n2) -# if extreme_type == "max" -# else min(d_lat_rad, lat_n1, lat_n2) -# ) -# else: -# return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2) +import numpy as np + +# 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.constants import ERROR_TOLERANCE + + +def _to_list(obj): + if not isinstance(obj, list): + if isinstance(obj, np.ndarray): + # Convert the NumPy array to a list using .tolist() + obj = obj.tolist() + else: + # If not a list or NumPy array, return the object as-is + obj = [obj] + return obj + + +def point_within_gca(pt, gca_cart, is_directed=False): + """Check if a point lies on a given Great Circle Arc (GCA). The anti- + meridian case is also considered. + + Parameters + ---------- + pt : numpy.ndarray (float) + Cartesian coordinates of the point. + gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr) + Cartesian coordinates of the Great Circle Arc (GCR). + is_directed : bool, optional, default = False + If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected, + and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False. + For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1 + + Returns + ------- + bool + True if the point lies between the two endpoints of the GCR, False otherwise. + + Raises + ------ + ValueError + If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes. + In such cases, consider breaking the GCR into two separate GCRs. + + ValueError + If the input GCR spans more than 180 degrees (π radians). + In such cases, consider breaking the GCR into two separate GCRs. + + Notes + ----- + The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and + accounting for the anti-meridian case. + + The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude). + In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation. + + The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations. + + Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. + """ + # Convert the cartesian coordinates to lonlat coordinates + pt_lonlat = _xyz_to_lonlat_rad(pt[0], pt[1], pt[2], scalar=True) + GCRv0_lonlat = _xyz_to_lonlat_rad( + gca_cart[0][0], gca_cart[0][1], gca_cart[0][2], scalar=True + ) + GCRv1_lonlat = _xyz_to_lonlat_rad( + gca_cart[1][0], gca_cart[1][1], gca_cart[1][2], scalar=True + ) + + # Convert the list to np.float64 + gca_cart[0] = np.array(gca_cart[0], dtype=np.float64) + gca_cart[1] = np.array(gca_cart[1], dtype=np.float64) + + # First if the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes + angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1]) + if np.allclose(angle, np.pi, rtol=0, atol=ERROR_TOLERANCE): + raise ValueError( + "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. " + "Consider breaking the Great Circle Arc" + "into two Great Circle Arcs" + ) + + if not np.allclose( + np.dot(np.cross(gca_cart[0], gca_cart[1]), pt), 0, rtol=0, atol=ERROR_TOLERANCE + ): + return False + + if np.isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): + # If the pt and the GCA are on the same longitude (the y coordinates are the same) + if np.isclose(GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE): + # Now use the latitude to determine if the pt falls between the interval + return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1]) + else: + # If the pt and the GCA are not on the same longitude when the GCA is a longnitude arc, then the pt is not on the GCA + return False + + # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point + if np.isclose( + abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0, atol=ERROR_TOLERANCE + ): + # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0] + if np.isclose(abs(pt_lonlat[1]), np.pi / 2, rtol=0, atol=ERROR_TOLERANCE): + pt_lonlat[0] = GCRv0_lonlat[0] + if not np.isclose( + GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE + ) and not np.isclose( + GCRv1_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE + ): + return False + else: + # Determine the pole latitude and latitude extension + if (GCRv0_lonlat[1] > 0 and GCRv1_lonlat[1] > 0) or ( + GCRv0_lonlat[1] < 0 and GCRv1_lonlat[1] < 0 + ): + pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0 else -np.pi / 2 + lat_extend = ( + abs(np.pi / 2 - abs(GCRv0_lonlat[1])) + + np.pi / 2 + + abs(GCRv1_lonlat[1]) + ) + else: + pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1]) + lat_extend = ( + abs(np.pi / 2 - abs(GCRv0_lonlat[1])) + + np.pi / 2 + + abs(GCRv1_lonlat[1]) + ) + + if is_directed and lat_extend >= np.pi: + raise ValueError( + "The input GCA spans more than 180 degrees. Please divide it into two." + ) + + return in_between(GCRv0_lonlat[1], pt_lonlat[1], pole_lat) or in_between( + pole_lat, pt_lonlat[1], GCRv1_lonlat[1] + ) + + if is_directed: + # The anti-meridian case Sufficient condition: absolute difference between the longitudes of the two + # vertices is greater than 180 degrees (π radians): abs(GCRv1_lon - GCRv0_lon) > π + if abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]) > np.pi: + # The necessary condition: the pt longitude is on the opposite side of the anti-meridian + # Case 1: where 0 --> x0--> 180 -->x1 -->0 case is lager than the 180degrees (pi radians) + if GCRv0_lonlat[0] <= np.pi <= GCRv1_lonlat[0]: + raise ValueError( + "The input Great Circle Arc span is larger than 180 degree, please break it into two." + ) + + # The necessary condition: the pt longitude is on the opposite side of the anti-meridian + # Case 2: The anti-meridian case where 180 -->x0 --> 0 lon --> x1 --> 180 lon + elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0: + return in_between( + GCRv0_lonlat[0], pt_lonlat[0], 2 * np.pi + ) or in_between(0, pt_lonlat[0], GCRv1_lonlat[0]) + + # The non-anti-meridian case. + else: + return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) + else: + # The undirected case + # sort the longitude + GCRv0_lonlat_min, GCRv1_lonlat_max = sorted([GCRv0_lonlat[0], GCRv1_lonlat[0]]) + if np.pi > GCRv1_lonlat_max - GCRv0_lonlat_min >= 0.0: + return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0]) + else: + return in_between(GCRv1_lonlat_max, pt_lonlat[0], 2 * np.pi) or in_between( + 0.0, pt_lonlat[0], GCRv0_lonlat_min + ) + + +def in_between(p, q, r) -> bool: + """Determines whether the number q is between p and r. + + Parameters + ---------- + p : float + The lower bound. + q : float + The number to check. + r : float + The upper bound. + + Returns + ------- + bool + True if q is between p and r, False otherwise. + """ + + return p <= q <= r or r <= q <= p + + +def _decide_pole_latitude(lat1, lat2): + """Determine the pole latitude based on the latitudes of two points on a + Great Circle Arc (GCA). + + This function calculates the combined latitude span from each point to its nearest pole + and decides which pole (North or South) the smaller GCA will pass. This decision is crucial + for handling GCAs that span exactly or more than 180 degrees in longitude, indicating + the arc might pass through or close to one of the Earth's poles. + + Parameters + ---------- + lat1 : float + Latitude of the first point in radians. Positive for the Northern Hemisphere, negative for the Southern. + lat2 : float + Latitude of the second point in radians. Positive for the Northern Hemisphere, negative for the Southern. + + Returns + ------- + float + The latitude of the pole (np.pi/2 for the North Pole or -np.pi/2 for the South Pole) the GCA is closer to. + + Notes + ----- + The function assumes the input latitudes are valid (i.e., between -np.pi/2 and np.pi/2) and expressed in radians. + The determination of which pole a GCA is closer to is based on the sum of the latitudinal spans from each point + to its nearest pole, considering the shortest path on the sphere. + """ + # Calculate the total latitudinal span to the nearest poles + lat_extend = abs(np.pi / 2 - abs(lat1)) + np.pi / 2 + abs(lat2) + + # Determine the closest pole based on the latitudinal span + if lat_extend < np.pi: + closest_pole = np.pi / 2 if lat1 > 0 else -np.pi / 2 + else: + closest_pole = -np.pi / 2 if lat1 > 0 else np.pi / 2 + + return closest_pole + + +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. + + Parameters + ---------- + gca_cart : numpy.ndarray + An array containing two 3D vectors that define a great circle arc. + + extreme_type : str + The type of extreme latitude to calculate. Must be either 'max' or 'min'. + + Returns + ------- + float + The maximum or minimum latitude of the great circle arc in radians. + + Raises + ------ + ValueError + If `extreme_type` is not 'max' or 'min'. + """ + extreme_type = extreme_type.lower() + + if extreme_type not in ("max", "min"): + raise ValueError("extreme_type must be either 'max' or 'min'") + + n1, n2 = gca_cart + dot_n1_n2 = np.dot(n1, n2) + denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0) + d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom + + d_a_max = ( + np.clip(d_a_max, 0, 1) + if np.isclose(d_a_max, [0, 1], atol=ERROR_TOLERANCE).any() + else d_a_max + ) + + _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2]) + _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2]) + + if 0 < d_a_max < 1: + node3 = (1 - d_a_max) * n1 + d_a_max * n2 + node3 = np.array(_normalize_xyz(node3[0], node3[1], node3[2])) + d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1)) + + return ( + max(d_lat_rad, lat_n1, lat_n2) + if extreme_type == "max" + else min(d_lat_rad, lat_n1, lat_n2) + ) + else: + return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2) diff --git a/uxarray/grid/area.py b/uxarray/grid/area.py index 5afa72e69..6e1692feb 100644 --- a/uxarray/grid/area.py +++ b/uxarray/grid/area.py @@ -1,14 +1,8 @@ import numpy as np -from numba import njit, config -from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT from uxarray.grid.coordinates import _lonlat_rad_to_xyz -config.DISABLE_JIT = not ENABLE_JIT - - -@njit(cache=ENABLE_JIT_CACHE) def calculate_face_area( x, y, z, quadrature_rule="gaussian", order=4, coords_type="spherical" ): @@ -67,13 +61,14 @@ def calculate_face_area( node3 = np.array([x[j + 2], y[j + 2], z[j + 2]], dtype=np.float64) if coords_type == "spherical": - node1 = np.array(_lonlat_rad_to_xyz(np.deg2rad(x[0]), np.deg2rad(y[0]))) - node2 = np.array( - _lonlat_rad_to_xyz(np.deg2rad(x[j + 1]), np.deg2rad(y[j + 1])) - ) - node3 = np.array( - _lonlat_rad_to_xyz(np.deg2rad(x[j + 2]), np.deg2rad(y[j + 2])) - ) + node1 = _lonlat_rad_to_xyz(np.deg2rad(x[0]), np.deg2rad(y[0])) + node1 = np.asarray(node1) + + node2 = _lonlat_rad_to_xyz(np.deg2rad(x[j + 1]), np.deg2rad(y[j + 1])) + node2 = np.asarray(node2) + + node3 = _lonlat_rad_to_xyz(np.deg2rad(x[j + 2]), np.deg2rad(y[j + 2])) + node3 = np.asarray(node3) for p in range(len(dW)): if quadrature_rule == "gaussian": @@ -97,7 +92,6 @@ def calculate_face_area( return area, jacobian -@njit(cache=ENABLE_JIT_CACHE) def get_all_face_area_from_coords( x, y, @@ -172,7 +166,6 @@ def get_all_face_area_from_coords( return area, jacobian -@njit(cache=ENABLE_JIT_CACHE) def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB): """Calculate Jacobian of a spherical triangle. This is a helper function for calculating face area. @@ -262,7 +255,6 @@ def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB): return dJacobian -@njit(cache=ENABLE_JIT_CACHE) def calculate_spherical_triangle_jacobian_barycentric(node1, node2, node3, dA, dB): """Calculate Jacobian of a spherical triangle. This is a helper function for calculating face area. @@ -341,7 +333,6 @@ def calculate_spherical_triangle_jacobian_barycentric(node1, node2, node3, dA, d return 0.5 * dJacobian -@njit(cache=ENABLE_JIT_CACHE) def get_gauss_quadratureDG(nCount): """Gauss Quadrature Points for integration. @@ -586,7 +577,6 @@ def get_gauss_quadratureDG(nCount): return dG, dW -@njit(cache=ENABLE_JIT_CACHE) def get_tri_quadratureDG(nOrder): """Triangular Quadrature Points for integration. diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 228c206f6..bec20f771 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -18,6 +18,7 @@ def _lonlat_rad_to_xyz( x = np.cos(lon) * np.cos(lat) y = np.sin(lon) * np.cos(lat) z = np.sin(lat) + return x, y, z @@ -29,6 +30,8 @@ def _xyz_to_lonlat_rad( ) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: """Converts Cartesian x, y, z coordinates in Spherical latitude and longitude coordinates.""" + + x, y, z = _normalize_xyz(x, y, z) denom = np.abs(x * x + y * y + z * z) x /= denom @@ -70,7 +73,6 @@ def _normalize_xyz( x: Union[np.ndarray, float], y: Union[np.ndarray, float], z: Union[np.ndarray, float], - scalar: bool = False, ) -> tuple[ Union[np.ndarray, float], Union[np.ndarray, float], Union[np.ndarray, float] ]: @@ -84,13 +86,12 @@ def _normalize_xyz( def _populate_node_latlon(grid) -> None: """Docstring TODO.""" - x_norm, y_norm, z_norm = _normalize_xyz( + lon_rad, lat_rad = _xyz_to_lonlat_rad( grid.node_x.values, grid.node_y.values, grid.node_z.values ) - lon_rad, lat_rad = _xyz_to_lonlat_rad(x_norm, y_norm, z_norm) # set longitude range to [-pi, pi] - lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi + # lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi # convert to degrees lon = np.rad2deg(lon_rad) diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index 0930fa6dc..058613514 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -1,353 +1,361 @@ -# import numpy as np -# from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE -# from uxarray.grid.intersections import gca_constLat_intersection -# -# # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad -# import pandas as pd -# -# DUMMY_EDGE_VALUE = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE] -# -# -# def _get_zonal_faces_weight_at_constLat( -# faces_edges_cart, -# latitude_cart, -# face_latlon_bound, -# is_directed=False, -# is_latlonface=False, -# is_face_GCA_list=None, -# ): -# """Utilize the sweep line algorithm to calculate the weight of each face at -# a constant latitude. -# -# Parameters -# ---------- -# face_edges_cart : np.ndarray -# A list of face polygon represented by edges in Cartesian coordinates. Shape: (n_faces, n_edges, 2, 3) -# -# latitude_cart : float -# The latitude in Cartesian coordinates (The normalized z coordinate) -# -# face_latlon_bound : np.ndarray -# The list of latitude and longitude bounds of faces. Shape: (n_faces,2, 2), -# [...,[lat_min, lat_max], [lon_min, lon_max],...] -# -# is_directed : bool, optional (default=False) -# If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected, -# and we will always assume the small circle (The one less than 180 degree) side is the GCA. -# -# is_latlonface : bool, optional, default=False -# A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means -# That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. -# Default is False. -# -# is_face_GCA_list : np.ndarray, optional (default=None) -# A list of boolean values that indicates if the edge in that face is a GCA. Shape: (n_faces,n_edges). -# True means edge face is a GCA. -# False mean this edge is a constant latitude. -# If None, all edges are considered as GCA. This attribute will overwrite the is_latlonface attribute. -# -# Returns -# ------- -# weights_df : pandas.DataFrame -# A DataFrame with the calculated weights of each face. The DataFrame has two columns: -# - 'face_index': The index of the face (integer). -# - 'weight': The calculated weight of the face in radian (float). -# The DataFrame is indexed by the face indices, providing a mapping from each face to its corresponding weight. -# """ -# intervals_list = [] -# -# # Iterate through all faces and their edges -# for face_index, face_edges in enumerate(faces_edges_cart): -# if is_face_GCA_list is not None: -# is_GCA_list = is_face_GCA_list[face_index] -# else: -# is_GCA_list = None -# face_interval_df = _get_zonal_face_interval( -# face_edges, -# latitude_cart, -# face_latlon_bound[face_index], -# is_directed=is_directed, -# is_latlonface=is_latlonface, -# is_GCA_list=is_GCA_list, -# ) -# for _, row in face_interval_df.iterrows(): -# intervals_list.append( -# {"start": row["start"], "end": row["end"], "face_index": face_index} -# ) -# -# intervals_df = pd.DataFrame(intervals_list) -# overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) -# -# # Calculate weights for each face -# weights = { -# face_index: overlap_contributions.get(face_index, 0.0) / total_length -# for face_index in range(len(faces_edges_cart)) -# } -# -# # Convert weights to DataFrame -# weights_df = pd.DataFrame(list(weights.items()), columns=["face_index", "weight"]) -# return weights_df -# -# -# def _is_edge_gca(is_GCA_list, is_latlonface, edges_z): -# """Determine if each edge is a Great Circle Arc (GCA) or a constant -# latitude line in a vectorized manner. -# -# Parameters: -# ---------- -# is_GCA_list : np.ndarray or None -# An array indicating whether each edge is a GCA (True) or a constant latitude line (False). -# Shape: (n_edges). If None, edge types are determined based on `is_latlonface` and the z-coordinates. -# is_latlonface : bool -# Flag indicating if all edges should be considered as lat-lon faces, which implies all edges -# are either constant latitude or longitude lines. -# edges_z : np.ndarray -# Array containing the z-coordinates for each vertex of the edges. This is used to determine -# whether edges are on the equator or if they are aligned in latitude when `is_GCA_list` is None. -# Shape should be (n_edges, 2). -# -# Returns: -# ------- -# np.ndarray -# A boolean array where each element indicates whether the corresponding edge is considered a GCA. -# True for GCA, False for constant latitude line. -# """ -# if is_GCA_list is not None: -# return is_GCA_list -# if is_latlonface: -# return ~np.isclose(edges_z[:, 0], edges_z[:, 1], atol=ERROR_TOLERANCE) -# return ~( -# np.isclose(edges_z[:, 0], 0, atol=ERROR_TOLERANCE) -# & np.isclose(edges_z[:, 1], 0, atol=ERROR_TOLERANCE) -# ) -# -# -# def _get_faces_constLat_intersection_info( -# face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed -# ): -# """Processes each edge of a face polygon in a vectorized manner to -# determine overlaps and calculate the intersections for a given latitude and -# the faces. -# -# Parameters: -# ---------- -# face_edges_cart : np.ndarray -# A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3). -# latitude_cart : float -# The latitude in Cartesian coordinates to which intersections or overlaps are calculated. -# is_GCA_list : np.ndarray or None -# An array indicating whether each edge is a GCA (True) or a constant latitude line (False). -# Shape: (n_edges). If None, the function will determine edge types based on `is_latlonface`. -# is_latlonface : bool -# Flag indicating if all faces are considered as lat-lon faces, meaning all edges are either -# constant latitude or longitude lines. This parameter overwrites the `is_GCA_list` if set to True. -# is_directed : bool -# Flag indicating if the GCA should be considered as directed (from v0 to v1). If False, -# the smaller circle (less than 180 degrees) side of the GCA is used. -# -# Returns: -# ------- -# 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. -# - overlap_flag (bool): A boolean indicating if any overlap with the latitude was detected. -# - overlap_edge (np.ndarray or None): The edge that overlaps with the latitude, if any; otherwise, None. -# """ -# valid_edges_mask = ~(np.any(face_edges_cart == DUMMY_EDGE_VALUE, axis=(1, 2))) -# -# # Apply mask to filter out dummy edges -# valid_edges = face_edges_cart[valid_edges_mask] -# -# # Extract Z coordinates for edge determination -# edges_z = valid_edges[:, :, 2] -# -# # Determine if each edge is GCA or constant latitude -# is_GCA = _is_edge_gca(is_GCA_list, is_latlonface, edges_z) -# -# # Check overlap with latitude -# overlaps_with_latitude = np.all( -# np.isclose(edges_z, latitude_cart, atol=ERROR_TOLERANCE), axis=1 -# ) -# overlap_flag = np.any(overlaps_with_latitude & ~is_GCA) -# -# # Identify overlap edges if needed -# intersections_pts_list_cart = [] -# if overlap_flag: -# overlap_index = np.where(overlaps_with_latitude & ~is_GCA)[0][0] -# intersections_pts_list_cart.extend(valid_edges[overlap_index]) -# else: -# # Calculate intersections (assuming a batch-capable intersection function) -# for idx, edge in enumerate(valid_edges): -# if is_GCA[idx]: -# intersections = gca_constLat_intersection( -# edge, latitude_cart, is_directed=is_directed -# ) -# if intersections.size == 0: -# continue -# elif intersections.shape[0] == 2: -# intersections_pts_list_cart.extend(intersections) -# else: -# intersections_pts_list_cart.append(intersections[0]) -# -# # Handle unique intersections and check for convex or concave cases -# unique_intersections = np.unique(intersections_pts_list_cart, axis=0) -# if len(unique_intersections) == 2: -# unique_intersection_lonlat = np.array( -# [node_xyz_to_lonlat_rad(pt.tolist()) for pt in unique_intersections] -# ) -# 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) != 0: -# raise ValueError( -# "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly" -# ) -# elif len(unique_intersections) == 0: -# raise ValueError( -# "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" -# ) -# -# -# def _get_zonal_face_interval( -# face_edges_cart, -# latitude_cart, -# face_latlon_bound, -# is_directed=False, -# is_latlonface=False, -# is_GCA_list=None, -# ): -# """Processes a face polygon represented by edges in Cartesian coordinates -# to find intervals where the face intersects with a given latitude. This -# function handles directed and undirected Great Circle Arcs (GCAs) and edges -# at constant latitude. -# -# Requires the face edges to be sorted in counter-clockwise order, and the span of the -# face in longitude should be less than pi. Also, all arcs/edges length should be within pi. -# -# Users can specify which edges are GCAs and which are constant latitude using `is_GCA_list`. -# However, edges on the equator are always treated as constant latitude edges regardless of -# `is_GCA_list`. -# -# Parameters -# ---------- -# face_edges_cart : np.ndarray -# A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3) -# latitude_cart : float -# The latitude in cartesian, the normalized Z coordinates. -# face_latlon_bound : np.ndarray -# The latitude and longitude bounds of the face. Shape: (2, 2), [[lat_min, lat_max], [lon_min, lon_max]] -# is_directed : bool, optional -# If True, the GCA is considered to be directed (from v0 to v1). If False, the GCA is undirected, -# and the smaller circle (less than 180 degrees) side of the GCA is used. Default is False. -# is_latlonface : bool, optional, default=False -# A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means -# That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. -# Default is False. This attribute will overwrite the is_latlonface attribute. -# is_GCA_list : np.ndarray, optional, default=False -# An array indicating if each edge is a GCA (True) or a constant latitude (False). Shape: (n_edges,). -# If None, all edges are considered as GCAs. Default is None. -# -# Returns -# ------- -# Intervals_df : pd.DataFrame -# A DataFrame containing the intervals stored as pandas Intervals for the face. -# The columns of the DataFrame are: ['start', 'end'] -# """ -# face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1] -# -# # Call the vectorized function to process all edges -# ( -# unique_intersections, -# pt_lon_min, -# pt_lon_max, -# ) = _get_faces_constLat_intersection_info( -# face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed -# ) -# -# # Convert intersection points to longitude-latitude -# longitudes = np.array( -# [node_xyz_to_lonlat_rad(pt.tolist())[0] for pt in unique_intersections] -# ) -# -# # Handle special wrap-around cases by checking the face bounds -# if face_lon_bound_left >= face_lon_bound_right: -# if not ( -# (pt_lon_max >= np.pi and pt_lon_min >= np.pi) -# or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi) -# ): -# # They're at different sides of the 0-lon, adding wrap-around points -# longitudes = np.append(longitudes, [0.0, 2 * np.pi]) -# -# # Ensure longitudes are sorted -# longitudes.sort() -# -# # Split the sorted longitudes into start and end points of intervals -# starts = longitudes[::2] # Start points -# ends = longitudes[1::2] # End points -# -# # Create the intervals DataFrame -# Intervals_df = pd.DataFrame({"start": starts, "end": ends}) -# # For consistency, sort the intervals by start -# interval_df_sorted = Intervals_df.sort_values(by="start").reset_index(drop=True) -# return interval_df_sorted -# -# -# def _process_overlapped_intervals(intervals_df): -# """Process the overlapped intervals using the sweep line algorithm, -# considering multiple intervals per face. -# -# Parameters -# ---------- -# intervals_df : pd.DataFrame -# The DataFrame that contains the intervals and the corresponding face index. -# The columns of the DataFrame are: ['start', 'end', 'face_index'] -# -# Returns -# ------- -# dict -# A dictionary where keys are face indices and values are their respective contributions to the total length. -# float -# The total length of all intervals considering their overlaps. -# -# Example -# ------- -# >>> intervals_data = [ -# ... {'start': 0.0, 'end': 100.0, 'face_index': 0}, -# ... {'start': 50.0, 'end': 150.0, 'face_index': 1}, -# ... {'start': 140.0, 'end': 150.0, 'face_index': 2} -# ... ] -# >>> intervals_df = pd.DataFrame(intervals_data) -# >>> overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) -# >>> print(overlap_contributions) -# >>> print(total_length) -# """ -# events = [] -# for idx, row in intervals_df.iterrows(): -# events.append((row["start"], "start", row["face_index"])) -# events.append((row["end"], "end", row["face_index"])) -# -# events.sort() # Sort by position and then by start/end -# -# active_faces = set() -# last_position = None -# total_length = 0.0 -# overlap_contributions = {} -# -# for position, event_type, face_idx in events: -# if last_position is not None and active_faces: -# segment_length = position - last_position -# segment_weight = segment_length / len(active_faces) if active_faces else 0 -# for active_face in active_faces: -# overlap_contributions[active_face] = ( -# overlap_contributions.get(active_face, 0) + segment_weight -# ) -# total_length += segment_length -# -# if event_type == "start": -# active_faces.add(face_idx) -# elif event_type == "end": -# active_faces.remove(face_idx) -# -# last_position = position -# -# return overlap_contributions, total_length +import numpy as np +from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE +from uxarray.grid.intersections import gca_constLat_intersection +from uxarray.grid.coordinates import _xyz_to_lonlat_rad +import pandas as pd + +DUMMY_EDGE_VALUE = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE] + + +def _get_zonal_faces_weight_at_constLat( + faces_edges_cart, + latitude_cart, + face_latlon_bound, + is_directed=False, + is_latlonface=False, + is_face_GCA_list=None, +): + """Utilize the sweep line algorithm to calculate the weight of each face at + a constant latitude. + + Parameters + ---------- + face_edges_cart : np.ndarray + A list of face polygon represented by edges in Cartesian coordinates. Shape: (n_faces, n_edges, 2, 3) + + latitude_cart : float + The latitude in Cartesian coordinates (The normalized z coordinate) + + face_latlon_bound : np.ndarray + The list of latitude and longitude bounds of faces. Shape: (n_faces,2, 2), + [...,[lat_min, lat_max], [lon_min, lon_max],...] + + is_directed : bool, optional (default=False) + If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected, + and we will always assume the small circle (The one less than 180 degree) side is the GCA. + + is_latlonface : bool, optional, default=False + A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means + That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. + Default is False. + + is_face_GCA_list : np.ndarray, optional (default=None) + A list of boolean values that indicates if the edge in that face is a GCA. Shape: (n_faces,n_edges). + True means edge face is a GCA. + False mean this edge is a constant latitude. + If None, all edges are considered as GCA. This attribute will overwrite the is_latlonface attribute. + + Returns + ------- + weights_df : pandas.DataFrame + A DataFrame with the calculated weights of each face. The DataFrame has two columns: + - 'face_index': The index of the face (integer). + - 'weight': The calculated weight of the face in radian (float). + The DataFrame is indexed by the face indices, providing a mapping from each face to its corresponding weight. + """ + intervals_list = [] + + # Iterate through all faces and their edges + for face_index, face_edges in enumerate(faces_edges_cart): + if is_face_GCA_list is not None: + is_GCA_list = is_face_GCA_list[face_index] + else: + is_GCA_list = None + face_interval_df = _get_zonal_face_interval( + face_edges, + latitude_cart, + face_latlon_bound[face_index], + is_directed=is_directed, + is_latlonface=is_latlonface, + is_GCA_list=is_GCA_list, + ) + for _, row in face_interval_df.iterrows(): + intervals_list.append( + {"start": row["start"], "end": row["end"], "face_index": face_index} + ) + + intervals_df = pd.DataFrame(intervals_list) + overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) + + # Calculate weights for each face + weights = { + face_index: overlap_contributions.get(face_index, 0.0) / total_length + for face_index in range(len(faces_edges_cart)) + } + + # Convert weights to DataFrame + weights_df = pd.DataFrame(list(weights.items()), columns=["face_index", "weight"]) + return weights_df + + +def _is_edge_gca(is_GCA_list, is_latlonface, edges_z): + """Determine if each edge is a Great Circle Arc (GCA) or a constant + latitude line in a vectorized manner. + + Parameters: + ---------- + is_GCA_list : np.ndarray or None + An array indicating whether each edge is a GCA (True) or a constant latitude line (False). + Shape: (n_edges). If None, edge types are determined based on `is_latlonface` and the z-coordinates. + is_latlonface : bool + Flag indicating if all edges should be considered as lat-lon faces, which implies all edges + are either constant latitude or longitude lines. + edges_z : np.ndarray + Array containing the z-coordinates for each vertex of the edges. This is used to determine + whether edges are on the equator or if they are aligned in latitude when `is_GCA_list` is None. + Shape should be (n_edges, 2). + + Returns: + ------- + np.ndarray + A boolean array where each element indicates whether the corresponding edge is considered a GCA. + True for GCA, False for constant latitude line. + """ + if is_GCA_list is not None: + return is_GCA_list + if is_latlonface: + return ~np.isclose(edges_z[:, 0], edges_z[:, 1], atol=ERROR_TOLERANCE) + return ~( + np.isclose(edges_z[:, 0], 0, atol=ERROR_TOLERANCE) + & np.isclose(edges_z[:, 1], 0, atol=ERROR_TOLERANCE) + ) + + +def _get_faces_constLat_intersection_info( + face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed +): + """Processes each edge of a face polygon in a vectorized manner to + determine overlaps and calculate the intersections for a given latitude and + the faces. + + Parameters: + ---------- + face_edges_cart : np.ndarray + A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3). + latitude_cart : float + The latitude in Cartesian coordinates to which intersections or overlaps are calculated. + is_GCA_list : np.ndarray or None + An array indicating whether each edge is a GCA (True) or a constant latitude line (False). + Shape: (n_edges). If None, the function will determine edge types based on `is_latlonface`. + is_latlonface : bool + Flag indicating if all faces are considered as lat-lon faces, meaning all edges are either + constant latitude or longitude lines. This parameter overwrites the `is_GCA_list` if set to True. + is_directed : bool + Flag indicating if the GCA should be considered as directed (from v0 to v1). If False, + the smaller circle (less than 180 degrees) side of the GCA is used. + + Returns: + ------- + 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. + - overlap_flag (bool): A boolean indicating if any overlap with the latitude was detected. + - overlap_edge (np.ndarray or None): The edge that overlaps with the latitude, if any; otherwise, None. + """ + valid_edges_mask = ~(np.any(face_edges_cart == DUMMY_EDGE_VALUE, axis=(1, 2))) + + # Apply mask to filter out dummy edges + valid_edges = face_edges_cart[valid_edges_mask] + + # Extract Z coordinates for edge determination + edges_z = valid_edges[:, :, 2] + + # Determine if each edge is GCA or constant latitude + is_GCA = _is_edge_gca(is_GCA_list, is_latlonface, edges_z) + + # Check overlap with latitude + overlaps_with_latitude = np.all( + np.isclose(edges_z, latitude_cart, atol=ERROR_TOLERANCE), axis=1 + ) + overlap_flag = np.any(overlaps_with_latitude & ~is_GCA) + + # Identify overlap edges if needed + intersections_pts_list_cart = [] + if overlap_flag: + overlap_index = np.where(overlaps_with_latitude & ~is_GCA)[0][0] + intersections_pts_list_cart.extend(valid_edges[overlap_index]) + else: + # Calculate intersections (assuming a batch-capable intersection function) + for idx, edge in enumerate(valid_edges): + if is_GCA[idx]: + intersections = gca_constLat_intersection( + edge, latitude_cart, is_directed=is_directed + ) + if intersections.size == 0: + continue + elif intersections.shape[0] == 2: + intersections_pts_list_cart.extend(intersections) + else: + intersections_pts_list_cart.append(intersections[0]) + + # Handle unique intersections and check for convex or concave cases + unique_intersections = np.unique(intersections_pts_list_cart, axis=0) + if len(unique_intersections) == 2: + # TODO: vectorize? + unique_intersection_lonlat = np.array( + [ + _xyz_to_lonlat_rad(pt[0], pt[1], pt[2], scalar=True) + for pt in unique_intersections + ] + ) + + 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) != 0: + raise ValueError( + "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly" + ) + elif len(unique_intersections) == 0: + raise ValueError( + "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" + ) + + +def _get_zonal_face_interval( + face_edges_cart, + latitude_cart, + face_latlon_bound, + is_directed=False, + is_latlonface=False, + is_GCA_list=None, +): + """Processes a face polygon represented by edges in Cartesian coordinates + to find intervals where the face intersects with a given latitude. This + function handles directed and undirected Great Circle Arcs (GCAs) and edges + at constant latitude. + + Requires the face edges to be sorted in counter-clockwise order, and the span of the + face in longitude should be less than pi. Also, all arcs/edges length should be within pi. + + Users can specify which edges are GCAs and which are constant latitude using `is_GCA_list`. + However, edges on the equator are always treated as constant latitude edges regardless of + `is_GCA_list`. + + Parameters + ---------- + face_edges_cart : np.ndarray + A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3) + latitude_cart : float + The latitude in cartesian, the normalized Z coordinates. + face_latlon_bound : np.ndarray + The latitude and longitude bounds of the face. Shape: (2, 2), [[lat_min, lat_max], [lon_min, lon_max]] + is_directed : bool, optional + If True, the GCA is considered to be directed (from v0 to v1). If False, the GCA is undirected, + and the smaller circle (less than 180 degrees) side of the GCA is used. Default is False. + is_latlonface : bool, optional, default=False + A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means + That all edge is either a longitude or constant latitude line. If False, then all edges are GCA. + Default is False. This attribute will overwrite the is_latlonface attribute. + is_GCA_list : np.ndarray, optional, default=False + An array indicating if each edge is a GCA (True) or a constant latitude (False). Shape: (n_edges,). + If None, all edges are considered as GCAs. Default is None. + + Returns + ------- + Intervals_df : pd.DataFrame + A DataFrame containing the intervals stored as pandas Intervals for the face. + The columns of the DataFrame are: ['start', 'end'] + """ + face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1] + + # Call the vectorized function to process all edges + ( + unique_intersections, + pt_lon_min, + pt_lon_max, + ) = _get_faces_constLat_intersection_info( + face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed + ) + + # Convert intersection points to longitude-latitude + # TODO: Vectorize + longitudes = np.array( + [ + _xyz_to_lonlat_rad(pt[0], pt[1], pt[2], scalar=True) + for pt in unique_intersections + ] + ) + + # Handle special wrap-around cases by checking the face bounds + if face_lon_bound_left >= face_lon_bound_right: + if not ( + (pt_lon_max >= np.pi and pt_lon_min >= np.pi) + or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi) + ): + # They're at different sides of the 0-lon, adding wrap-around points + longitudes = np.append(longitudes, [0.0, 2 * np.pi]) + + # Ensure longitudes are sorted + longitudes.sort() + + # Split the sorted longitudes into start and end points of intervals + starts = longitudes[::2] # Start points + ends = longitudes[1::2] # End points + + # Create the intervals DataFrame + Intervals_df = pd.DataFrame({"start": starts, "end": ends}) + # For consistency, sort the intervals by start + interval_df_sorted = Intervals_df.sort_values(by="start").reset_index(drop=True) + return interval_df_sorted + + +def _process_overlapped_intervals(intervals_df): + """Process the overlapped intervals using the sweep line algorithm, + considering multiple intervals per face. + + Parameters + ---------- + intervals_df : pd.DataFrame + The DataFrame that contains the intervals and the corresponding face index. + The columns of the DataFrame are: ['start', 'end', 'face_index'] + + Returns + ------- + dict + A dictionary where keys are face indices and values are their respective contributions to the total length. + float + The total length of all intervals considering their overlaps. + + Example + ------- + >>> intervals_data = [ + ... {'start': 0.0, 'end': 100.0, 'face_index': 0}, + ... {'start': 50.0, 'end': 150.0, 'face_index': 1}, + ... {'start': 140.0, 'end': 150.0, 'face_index': 2} + ... ] + >>> intervals_df = pd.DataFrame(intervals_data) + >>> overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) + >>> print(overlap_contributions) + >>> print(total_length) + """ + events = [] + for idx, row in intervals_df.iterrows(): + events.append((row["start"], "start", row["face_index"])) + events.append((row["end"], "end", row["face_index"])) + + events.sort() # Sort by position and then by start/end + + active_faces = set() + last_position = None + total_length = 0.0 + overlap_contributions = {} + + for position, event_type, face_idx in events: + if last_position is not None and active_faces: + segment_length = position - last_position + segment_weight = segment_length / len(active_faces) if active_faces else 0 + for active_face in active_faces: + overlap_contributions[active_face] = ( + overlap_contributions.get(active_face, 0) + segment_weight + ) + total_length += segment_length + + if event_type == "start": + active_faces.add(face_idx) + elif event_type == "end": + active_faces.remove(face_idx) + + last_position = position + + return overlap_contributions, total_length From 594c7f1516e1deed0f2958f422c8a858a44f9b10 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 8 May 2024 14:01:29 -0500 Subject: [PATCH 04/14] fix tests (1) --- test/test_helpers.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/test/test_helpers.py b/test/test_helpers.py index ec310c346..42cd72def 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 node_lonlat_rad_to_xyz +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 @@ -105,23 +105,22 @@ def test_grid_center(self): class TestCoordinatesConversion(TestCase): def test_normalize_in_place(self): - [x, y, z] = ux.grid.coordinates.normalize_in_place( - [random.random(), random.random(), - random.random()]) + x, y, z = _normalize_xyz( + random.random(), random.random(), + random.random()) self.assertLessEqual(np.absolute(np.sqrt(x * x + y * y + z * z) - 1), err_tolerance) def test_node_xyz_to_lonlat_rad(self): - [x, y, z] = ux.grid.coordinates.normalize_in_place([ + x, y, z = _normalize_xyz(*[ random.uniform(-1, 1), random.uniform(-1, 1), random.uniform(-1, 1) ]) - [lon, lat] = ux.grid.coordinates.node_xyz_to_lonlat_rad([x, y, z]) - [new_x, new_y, - new_z] = ux.grid.coordinates.node_lonlat_rad_to_xyz([lon, lat]) + lon, lat = _xyz_to_lonlat_rad(x, y, z) + new_x, new_y, new_z =_lonlat_rad_to_xyz(lon, lat) self.assertLessEqual(np.absolute(new_x - x), err_tolerance) self.assertLessEqual(np.absolute(new_y - y), err_tolerance) @@ -133,10 +132,9 @@ def test_node_latlon_rad_to_xyz(self): random.uniform(-0.5 * np.pi, 0.5 * np.pi) ] - [x, y, z] = ux.grid.coordinates.node_lonlat_rad_to_xyz([lon, lat]) + x, y, z = _lonlat_rad_to_xyz(lon, lat) - [new_lon, - new_lat] = ux.grid.coordinates.node_xyz_to_lonlat_rad([x, y, z]) + new_lon, new_lat = _xyz_to_lonlat_rad(x, y, z) self.assertLessEqual(np.absolute(new_lon - lon), err_tolerance) self.assertLessEqual(np.absolute(new_lat - lat), err_tolerance) @@ -342,7 +340,7 @@ def test_get_lonlat_face_edge_nodes_pipeline(self): face_edges_connectivity_cartesian = [] for i in range(len(face_edges_connectivity_lonlat)): edge = face_edges_connectivity_lonlat[i] - edge_cart = [node_lonlat_rad_to_xyz(node) for node in edge] + edge_cart = [_lonlat_rad_to_xyz(*node) for node in edge] face_edges_connectivity_cartesian.append(edge_cart) # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon From 7abfc9334a97ace91fc0da51084a2d1fc7a1b3a1 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 8 May 2024 15:10:15 -0500 Subject: [PATCH 05/14] fix tests (2) --- test/test_centroids.py | 13 ++--- test/test_geometry.py | 94 ++++++++++++++++++++----------------- uxarray/grid/arcs.py | 10 ++-- uxarray/grid/coordinates.py | 18 +++---- 4 files changed, 72 insertions(+), 63 deletions(-) diff --git a/test/test_centroids.py b/test/test_centroids.py index 4b5de9fd4..6eb69542f 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_in_place +from uxarray.grid.coordinates import _populate_face_centroids, _populate_edge_centroids, _normalize_xyz current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -21,8 +21,8 @@ def test_centroids_from_mean_verts_triangle(self): # Calculate the expected centroid expected_centroid = np.mean(test_triangle, axis=0) - [norm_x, norm_y, norm_z] = normalize_in_place( - [expected_centroid[0], expected_centroid[1], expected_centroid[2]]) + norm_x, norm_y, norm_z = _normalize_xyz( + expected_centroid[0], expected_centroid[1], expected_centroid[2]) # Open the dataset and find the centroids grid = ux.open_grid(test_triangle) @@ -42,8 +42,8 @@ def test_centroids_from_mean_verts_pentagon(self): # Calculate the expected centroid expected_centroid = np.mean(test_triangle, axis=0) - [norm_x, norm_y, norm_z] = normalize_in_place( - [expected_centroid[0], expected_centroid[1], expected_centroid[2]]) + norm_x, norm_y, norm_z = _normalize_xyz( + expected_centroid[0], expected_centroid[1], expected_centroid[2]) # Open the dataset and find the centroids grid = ux.open_grid(test_triangle) @@ -65,7 +65,8 @@ def test_centroids_from_mean_verts_scrip(self): _populate_face_centroids(uxgrid, repopulate=True) - computed_face_x = (uxgrid.face_lon.values + 180) % 360 - 180 + # computed_face_x = (uxgrid.face_lon.values + 180) % 360 - 180 + computed_face_x = uxgrid.face_lon.values computed_face_y = uxgrid.face_lat.values nt.assert_array_almost_equal(expected_face_x, computed_face_x) diff --git a/test/test_geometry.py b/test/test_geometry.py index beaa105b7..8b03d2430 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -9,7 +9,7 @@ import uxarray as ux from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE import uxarray.utils.computing as ac_utils -from uxarray.grid.coordinates import _populate_lonlat_coord, node_lonlat_rad_to_xyz +from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad from uxarray.grid.arcs import extreme_gca_latitude from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds @@ -67,7 +67,7 @@ def test_pole_point_inside_polygon_from_vertice_north(self): # Normalize the vertices to ensure they lie on the unit sphere for i, vertex in enumerate(vertices): float_vertex = [float(coord) for coord in vertex] - vertices[i] = ux.grid.coordinates.normalize_in_place(float_vertex) + vertices[i] = _normalize_xyz(*float_vertex) # Create face_edge_cart from the vertices face_edge_cart = np.array([[vertices[0], vertices[1]], @@ -93,7 +93,7 @@ def test_pole_point_inside_polygon_from_vertice_south(self): # Normalize the vertices to ensure they lie on the unit sphere for i, vertex in enumerate(vertices): float_vertex = [float(coord) for coord in vertex] - vertices[i] = ux.grid.coordinates.normalize_in_place(float_vertex) + vertices[i] = _normalize_xyz(*float_vertex) # Create face_edge_cart from the vertices, since we are using the south pole, and want retrive the smaller face # we need to reverse the order of the vertices @@ -121,7 +121,7 @@ def test_pole_point_inside_polygon_from_vertice_pole(self): # Normalize the vertices to ensure they lie on the unit sphere for i, vertex in enumerate(vertices): float_vertex = [float(coord) for coord in vertex] - vertices[i] = ux.grid.coordinates.normalize_in_place(float_vertex) + vertices[i] = _normalize_xyz(*float_vertex) # Create face_edge_cart from the vertices face_edge_cart = np.array([[vertices[0], vertices[1]], @@ -147,7 +147,7 @@ def test_pole_point_inside_polygon_from_vertice_cross(self): # Normalize the vertices to ensure they lie on the unit sphere for i, vertex in enumerate(vertices): float_vertex = [float(coord) for coord in vertex] - vertices[i] = ux.grid.coordinates.normalize_in_place(float_vertex) + vertices[i] = _normalize_xyz(*float_vertex) # Create face_edge_cart from the vertices face_edge_cart = np.array([[vertices[0], vertices[1]], @@ -191,13 +191,13 @@ def _max_latitude_rad_iterative(self, gca_cart): # Convert input vectors to radians and Cartesian coordinates v1_cart, v2_cart = gca_cart - b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v1_cart.tolist()) - c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v2_cart.tolist()) + b_lonlat = _xyz_to_lonlat_rad(*v1_cart.tolist()) + c_lonlat = _xyz_to_lonlat_rad(*v2_cart.tolist()) # Initialize variables for the iterative process v_temp = ac_utils.cross_fma(v1_cart, v2_cart) v0 = ac_utils.cross_fma(v_temp, v1_cart) - v0 = ux.grid.coordinates.normalize_in_place(v0.tolist()) + v0 = _normalize_xyz(*v0.tolist()) max_section = [v1_cart, v2_cart] # Iteratively find the maximum latitude @@ -207,7 +207,7 @@ def _max_latitude_rad_iterative(self, gca_cart): v_b, v_c = max_section angle_v1_v2_rad = ux.grid.arcs._angle_of_2_vectors(v_b, v_c) v0 = ac_utils.cross_fma(v_temp, v_b) - v0 = ux.grid.coordinates.normalize_in_place(v0.tolist()) + v0 = _normalize_xyz(*v0.tolist()) avg_angle_rad = angle_v1_v2_rad / 10.0 for i in range(10): @@ -217,10 +217,13 @@ def _max_latitude_rad_iterative(self, gca_cart): angle_rad_prev) * np.array(v0) w2_new = np.cos(angle_rad_next) * v_b + np.sin( angle_rad_next) * np.array(v0) - w1_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad( - w1_new.tolist()) - w2_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad( - w2_new.tolist()) + w1_lonlat = _xyz_to_lonlat_rad( + *w1_new.tolist()) + w2_lonlat = _xyz_to_lonlat_rad( + *w2_new.tolist()) + + w1_lonlat = np.asarray(w1_lonlat) + w2_lonlat = np.asarray(w2_lonlat) # Adjust latitude boundaries to avoid error accumulation if i == 0: @@ -241,10 +244,10 @@ def _max_latitude_rad_iterative(self, gca_cart): max_section = [w1_new, w2_new] if i != 9 else [w1_new, v_c] # Update longitude and latitude for the next iteration - b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad( - max_section[0].tolist()) - c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad( - max_section[1].tolist()) + b_lonlat = _xyz_to_lonlat_rad( + *max_section[0].tolist()) + c_lonlat = _xyz_to_lonlat_rad( + *max_section[1].tolist()) return np.average([b_lonlat[1], c_lonlat[1]]) @@ -275,13 +278,13 @@ def _min_latitude_rad_iterative(self, gca_cart): # Convert input vectors to radians and Cartesian coordinates v1_cart, v2_cart = gca_cart - b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v1_cart.tolist()) - c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v2_cart.tolist()) + b_lonlat = _xyz_to_lonlat_rad(*v1_cart.tolist()) + c_lonlat = _xyz_to_lonlat_rad(*v2_cart.tolist()) # Initialize variables for the iterative process v_temp = ac_utils.cross_fma(v1_cart, v2_cart) v0 = ac_utils.cross_fma(v_temp, v1_cart) - v0 = np.array(ux.grid.coordinates.normalize_in_place(v0.tolist())) + v0 = np.array(_normalize_xyz(*v0.tolist())) min_section = [v1_cart, v2_cart] # Iteratively find the minimum latitude @@ -291,7 +294,7 @@ def _min_latitude_rad_iterative(self, gca_cart): v_b, v_c = min_section angle_v1_v2_rad = ux.grid.arcs._angle_of_2_vectors(v_b, v_c) v0 = ac_utils.cross_fma(v_temp, v_b) - v0 = np.array(ux.grid.coordinates.normalize_in_place(v0.tolist())) + v0 = np.array(_normalize_xyz(*v0.tolist())) avg_angle_rad = angle_v1_v2_rad / 10.0 for i in range(10): @@ -301,10 +304,13 @@ def _min_latitude_rad_iterative(self, gca_cart): angle_rad_prev) * v0 w2_new = np.cos(angle_rad_next) * v_b + np.sin( angle_rad_next) * v0 - w1_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad( - w1_new.tolist()) - w2_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad( - w2_new.tolist()) + w1_lonlat = _xyz_to_lonlat_rad( + *w1_new.tolist()) + w2_lonlat = _xyz_to_lonlat_rad( + *w2_new.tolist()) + + w1_lonlat = np.asarray(w1_lonlat) + w2_lonlat = np.asarray(w2_lonlat) # Adjust latitude boundaries to avoid error accumulation if i == 0: @@ -325,18 +331,18 @@ def _min_latitude_rad_iterative(self, gca_cart): min_section = [w1_new, w2_new] if i != 9 else [w1_new, v_c] # Update longitude and latitude for the next iteration - b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad( - min_section[0].tolist()) - c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad( - min_section[1].tolist()) + b_lonlat = _xyz_to_lonlat_rad( + *min_section[0].tolist()) + c_lonlat = _xyz_to_lonlat_rad( + *min_section[1].tolist()) return np.average([b_lonlat[1], c_lonlat[1]]) def test_extreme_gca_latitude_max(self): # Define a great circle arc that is symmetrical around 0 degrees longitude gca_cart = np.array([ - ux.grid.coordinates.normalize_in_place([0.5, 0.5, 0.5]), - ux.grid.coordinates.normalize_in_place([-0.5, 0.5, 0.5]) + _normalize_xyz(*[0.5, 0.5, 0.5]), + _normalize_xyz(*[-0.5, 0.5, 0.5]) ]) # Calculate the maximum latitude @@ -363,8 +369,8 @@ def test_extreme_gca_latitude_max(self): def test_extreme_gca_latitude_min(self): # Define a great circle arc that is symmetrical around 0 degrees longitude gca_cart = np.array([ - ux.grid.coordinates.normalize_in_place([0.5, 0.5, -0.5]), - ux.grid.coordinates.normalize_in_place([-0.5, 0.5, -0.5]) + _normalize_xyz(*[0.5, 0.5, -0.5]), + _normalize_xyz(*[-0.5, 0.5, -0.5]) ]) # Calculate the minimum latitude @@ -441,7 +447,7 @@ def test_populate_bounds_normal(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = max(np.deg2rad(60.0), extreme_gca_latitude(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")) lat_min = min(np.deg2rad(10.0), @@ -470,7 +476,7 @@ def test_populate_bounds_antimeridian(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = max(np.deg2rad(60.0), extreme_gca_latitude(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")) lat_min = min(np.deg2rad(10.0), @@ -499,7 +505,7 @@ def test_populate_bounds_node_on_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.pi / 2 lat_min = min(np.deg2rad(10.0), extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) @@ -527,7 +533,7 @@ def test_populate_bounds_edge_over_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.pi / 2 lat_min = min(np.deg2rad(60.0), extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) @@ -555,7 +561,7 @@ def test_populate_bounds_pole_inside(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.pi / 2 lat_min = min(np.deg2rad(60.0), extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) @@ -666,7 +672,7 @@ def test_populate_bounds_edge_over_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.pi / 2 lat_min = np.deg2rad(60.0) lon_min = np.deg2rad(210.0) @@ -694,7 +700,7 @@ def test_populate_bounds_pole_inside(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.pi / 2 lat_min = np.deg2rad(60.0) lon_min = 0 @@ -752,7 +758,7 @@ def test_populate_bounds_antimeridian(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.deg2rad(60.0) lat_min = np.deg2rad(10.0) lon_min = np.deg2rad(350.0) @@ -780,7 +786,7 @@ def test_populate_bounds_node_on_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.pi / 2 lat_min = np.deg2rad(10.0) lon_min = np.deg2rad(10.0) @@ -808,7 +814,7 @@ def test_populate_bounds_edge_over_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.pi / 2 lat_min = np.deg2rad(60.0) lon_min = np.deg2rad(210.0) @@ -836,7 +842,7 @@ def test_populate_bounds_pole_inside(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad] + vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] lat_max = np.pi / 2 lat_min = np.deg2rad(60.0) lon_min = 0 diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 82954ccab..ec18195b4 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -60,12 +60,12 @@ def point_within_gca(pt, gca_cart, is_directed=False): Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. """ # Convert the cartesian coordinates to lonlat coordinates - pt_lonlat = _xyz_to_lonlat_rad(pt[0], pt[1], pt[2], scalar=True) - GCRv0_lonlat = _xyz_to_lonlat_rad( - gca_cart[0][0], gca_cart[0][1], gca_cart[0][2], scalar=True + pt_lonlat = np.array(_xyz_to_lonlat_rad(pt[0], pt[1], pt[2], scalar=True)) + GCRv0_lonlat = np.array( + _xyz_to_lonlat_rad(gca_cart[0][0], gca_cart[0][1], gca_cart[0][2], scalar=True) ) - GCRv1_lonlat = _xyz_to_lonlat_rad( - gca_cart[1][0], gca_cart[1][1], gca_cart[1][2], scalar=True + GCRv1_lonlat = np.array( + _xyz_to_lonlat_rad(gca_cart[1][0], gca_cart[1][1], gca_cart[1][2], scalar=True) ) # Convert the list to np.float64 diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index bec20f771..b2b027d1e 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -61,11 +61,13 @@ def _xyz_to_lonlat_deg( ) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z) - # set lon range to [-pi, pi] - lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi + # # set lon range to [-pi, pi] + # lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi lon = np.rad2deg(lon_rad) lat = np.rad2deg(lat_rad) + + lon = (lon + 180) % 360 - 180 return lon, lat @@ -76,7 +78,7 @@ def _normalize_xyz( ) -> tuple[ Union[np.ndarray, float], Union[np.ndarray, float], Union[np.ndarray, float] ]: - denom = np.linalg.norm([x, y, z], ord=2) + denom = np.linalg.norm(np.asarray([x, y, z]), ord=2) x_norm = x / denom y_norm = y / denom @@ -146,9 +148,9 @@ def _build_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face): def _build_edge_centroids(node_x, node_y, node_z, edge_nodes): - centroid_x = np.mean(node_x[edge_nodes]) - centroid_y = np.mean(node_y[edge_nodes]) - centroid_z = np.mean(node_z[edge_nodes]) + centroid_x = np.mean(node_x[edge_nodes], axis=1) + centroid_y = np.mean(node_y[edge_nodes], axis=1) + centroid_z = np.mean(node_z[edge_nodes], axis=1) return _normalize_xyz(centroid_x, centroid_y, centroid_z) @@ -172,7 +174,7 @@ def _populate_face_centroids(grid, repopulate=False): centroid_x, centroid_y, centroid_z = grid.face_x, grid.face_y, grid.face_z # Convert from xyz to latlon - centroid_lon, centroid_lat = _xyz_to_lonlat_rad( + centroid_lon, centroid_lat = _xyz_to_lonlat_deg( centroid_x, centroid_y, centroid_z ) else: @@ -233,7 +235,7 @@ def _populate_edge_centroids(grid, repopulate=False): centroid_x, centroid_y, centroid_z = grid.edge_x, grid.edge_y, grid.edge_z # Convert from xyz to latlon - centroid_lon, centroid_lat = _xyz_to_lonlat_rad( + centroid_lon, centroid_lat = _xyz_to_lonlat_deg( centroid_x, centroid_y, centroid_z ) else: From 4965d935f225c19de5f34893688a7b1a2e7e166b Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 8 May 2024 15:54:17 -0500 Subject: [PATCH 06/14] fix tests (3) --- test/test_integrate.py | 30 +++++++++++++++--------------- uxarray/grid/integrate.py | 6 +----- 2 files changed, 16 insertions(+), 20 deletions(-) diff --git a/test/test_integrate.py b/test/test_integrate.py index af26941c6..078f355cf 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -8,7 +8,7 @@ import numpy.testing as nt import uxarray as ux -from uxarray.grid.coordinates import node_lonlat_rad_to_xyz +from uxarray.grid.coordinates import _lonlat_rad_to_xyz from uxarray.grid.integrate import _get_zonal_face_interval, _process_overlapped_intervals, _get_zonal_faces_weight_at_constLat current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -67,7 +67,7 @@ def test_get_zonal_face_interval(self): [1.6 * np.pi, -0.25 * np.pi], [0.4 * np.pi, -0.25 * np.pi], [0.4 * np.pi, 0.25 * np.pi]] - vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat] + vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat] face_edge_nodes = np.array([[vertices[0], vertices[1]], [vertices[1], vertices[2]], @@ -101,7 +101,7 @@ def test_get_zonal_face_interval_GCA_constLat(self): [0.4 * np.pi, -0.25 * np.pi], [0.4 * np.pi, 0.25 * np.pi]] - vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat] + vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat] face_edge_nodes = np.array([[vertices[0], vertices[1]], [vertices[1], vertices[2]], @@ -132,7 +132,7 @@ def test_get_zonal_face_interval_equator(self): vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi], [-0.4 * np.pi, 0.0], [0.4 * np.pi, 0.0], [0.4 * np.pi, 0.25 * np.pi]] - vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat] + vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat] face_edge_nodes = np.array([[vertices[0], vertices[1]], [vertices[1], vertices[2]], @@ -284,10 +284,10 @@ def test_get_zonal_faces_weight_at_constLat_equator(self): [np.pi, -0.25 * np.pi], [np.pi, 0.0]] # Convert the face vertices to xyz coordinates - face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0] - face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1] - face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2] - face_3 = [node_lonlat_rad_to_xyz(v) for v in face_3] + face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0] + face_1 = [_lonlat_rad_to_xyz(*v) for v in face_1] + face_2 = [_lonlat_rad_to_xyz(*v) for v in face_2] + face_3 = [_lonlat_rad_to_xyz(*v) for v in face_3] face_0_edge_nodes = np.array([[face_0[0], face_0[1]], [face_0[1], face_0[2]], @@ -347,10 +347,10 @@ def test_get_zonal_faces_weight_at_constLat_regular(self): [1.6 * np.pi, -0.01 * np.pi], [1.6 * np.pi, 0.25 * np.pi]] # Convert the face vertices to xyz coordinates - face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0] - face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1] - face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2] - face_3 = [node_lonlat_rad_to_xyz(v) for v in face_3] + face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0] + face_1 = [_lonlat_rad_to_xyz(*v) for v in face_1] + face_2 = [_lonlat_rad_to_xyz(*v) for v in face_2] + face_3 = [_lonlat_rad_to_xyz(*v) for v in face_3] face_0_edge_nodes = np.array([[face_0[0], face_0[1]], [face_0[1], face_0[2]], @@ -410,9 +410,9 @@ def test_get_zonal_faces_weight_at_constLat_latlonface(self): [np.deg2rad(40), np.deg2rad(20)], [np.deg2rad(40), np.deg2rad(40)]] # Convert the face vertices to xyz coordinates - face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0] - face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1] - face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2] + face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0] + face_1 = [_lonlat_rad_to_xyz(*v) for v in face_1] + face_2 = [_lonlat_rad_to_xyz(*v) for v in face_2] face_0_edge_nodes = np.array([[face_0[0], face_0[1]], diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index 058613514..b18a25427 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -269,12 +269,8 @@ def _get_zonal_face_interval( ) # Convert intersection points to longitude-latitude - # TODO: Vectorize longitudes = np.array( - [ - _xyz_to_lonlat_rad(pt[0], pt[1], pt[2], scalar=True) - for pt in unique_intersections - ] + [_xyz_to_lonlat_rad(*pt.tolist())[0] for pt in unique_intersections] ) # Handle special wrap-around cases by checking the face bounds From cfab0ddd25384221f9007fd2d4240ec6e2f94faf Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 8 May 2024 16:40:21 -0500 Subject: [PATCH 07/14] fix tests (4) --- test/test_scrip.py | 39 ++++++++--------- test/test_ugrid.py | 38 ++++++++--------- uxarray/grid/coordinates.py | 85 +++++++++++++++++++++++++++++++------ 3 files changed, 112 insertions(+), 50 deletions(-) diff --git a/test/test_scrip.py b/test/test_scrip.py index b0c4426fb..9f234f6f7 100644 --- a/test/test_scrip.py +++ b/test/test_scrip.py @@ -38,25 +38,26 @@ def test_read_ugrid(self): nt.assert_equal(uxgrid_RLL10deg_ne4.node_lon.size, constants.NNODES_ov_RLL10deg_CSne4) - def test_read_ugrid_opendap(self): - """Read an ugrid model from an OPeNDAP URL.""" - - try: - # make sure we can read the ugrid file from the OPeNDAP URL - url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc" - uxgrid_url = ux.open_grid(url, drop_variables="siglay") - - except OSError: - # print warning and pass if we can't connect to the OPeNDAP server - warnings.warn(f'Could not connect to OPeNDAP server: {url}') - pass - - else: - - assert isinstance(getattr(uxgrid_url, "node_lon"), xr.DataArray) - assert isinstance(getattr(uxgrid_url, "node_lat"), xr.DataArray) - assert isinstance(getattr(uxgrid_url, "face_node_connectivity"), - xr.DataArray) + # TODO: UNCOMMENT + # def test_read_ugrid_opendap(self): + # """Read an ugrid model from an OPeNDAP URL.""" + # + # try: + # # make sure we can read the ugrid file from the OPeNDAP URL + # url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc" + # uxgrid_url = ux.open_grid(url, drop_variables="siglay") + # + # except OSError: + # # print warning and pass if we can't connect to the OPeNDAP server + # warnings.warn(f'Could not connect to OPeNDAP server: {url}') + # pass + # + # else: + # + # assert isinstance(getattr(uxgrid_url, "node_lon"), xr.DataArray) + # assert isinstance(getattr(uxgrid_url, "node_lat"), xr.DataArray) + # assert isinstance(getattr(uxgrid_url, "face_node_connectivity"), + # xr.DataArray) def test_encode_ugrid(self): """Read an Exodus dataset and encode that as a UGRID format.""" diff --git a/test/test_ugrid.py b/test/test_ugrid.py index b0c4426fb..9f7576ce2 100644 --- a/test/test_ugrid.py +++ b/test/test_ugrid.py @@ -38,25 +38,25 @@ def test_read_ugrid(self): nt.assert_equal(uxgrid_RLL10deg_ne4.node_lon.size, constants.NNODES_ov_RLL10deg_CSne4) - def test_read_ugrid_opendap(self): - """Read an ugrid model from an OPeNDAP URL.""" - - try: - # make sure we can read the ugrid file from the OPeNDAP URL - url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc" - uxgrid_url = ux.open_grid(url, drop_variables="siglay") - - except OSError: - # print warning and pass if we can't connect to the OPeNDAP server - warnings.warn(f'Could not connect to OPeNDAP server: {url}') - pass - - else: - - assert isinstance(getattr(uxgrid_url, "node_lon"), xr.DataArray) - assert isinstance(getattr(uxgrid_url, "node_lat"), xr.DataArray) - assert isinstance(getattr(uxgrid_url, "face_node_connectivity"), - xr.DataArray) + # def test_read_ugrid_opendap(self): + # """Read an ugrid model from an OPeNDAP URL.""" + # + # try: + # # make sure we can read the ugrid file from the OPeNDAP URL + # url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc" + # uxgrid_url = ux.open_grid(url, drop_variables="siglay") + # + # except OSError: + # # print warning and pass if we can't connect to the OPeNDAP server + # warnings.warn(f'Could not connect to OPeNDAP server: {url}') + # pass + # + # else: + # + # assert isinstance(getattr(uxgrid_url, "node_lon"), xr.DataArray) + # assert isinstance(getattr(uxgrid_url, "node_lat"), xr.DataArray) + # assert isinstance(getattr(uxgrid_url, "face_node_connectivity"), + # xr.DataArray) def test_encode_ugrid(self): """Read an Exodus dataset and encode that as a UGRID format.""" diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index b2b027d1e..6631ed284 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -1,6 +1,8 @@ import xarray as xr import numpy as np +import warnings + from uxarray.constants import ERROR_TOLERANCE from uxarray.conventions import ugrid @@ -27,16 +29,17 @@ def _xyz_to_lonlat_rad( y: Union[np.ndarray, float], z: Union[np.ndarray, float], scalar: bool = False, + normalize: bool = True, ) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: """Converts Cartesian x, y, z coordinates in Spherical latitude and longitude coordinates.""" - x, y, z = _normalize_xyz(x, y, z) - denom = np.abs(x * x + y * y + z * z) - - x /= denom - y /= denom - z /= denom + 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) lat = np.arcsin(z) @@ -58,8 +61,9 @@ def _xyz_to_lonlat_deg( x: Union[np.ndarray, float], y: Union[np.ndarray, float], z: Union[np.ndarray, float], + normalize: bool = True, ) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: - lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z) + lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z, normalize=normalize) # # set lon range to [-pi, pi] # lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi @@ -78,6 +82,7 @@ def _normalize_xyz( ) -> tuple[ Union[np.ndarray, float], Union[np.ndarray, float], Union[np.ndarray, float] ]: + # TODO denom = np.linalg.norm(np.asarray([x, y, z]), ord=2) x_norm = x / denom @@ -156,6 +161,18 @@ def _build_edge_centroids(node_x, node_y, node_z, edge_nodes): def _populate_face_centroids(grid, repopulate=False): + """Finds the centroids of faces using cartesian averaging based off the + vertices. The centroid is defined as the average of the x, y, z + coordinates, normalized. This cannot be guaranteed to work on concave + polygons. + + Parameters + ---------- + repopulate : bool, optional + Bool used to turn on/off repopulating the face coordinates of the centroids + """ + warnings.warn("This cannot be guaranteed to work correctly on concave polygons") + node_x = grid.node_x.values node_y = grid.node_y.values node_z = grid.node_z.values @@ -165,7 +182,7 @@ def _populate_face_centroids(grid, repopulate=False): if "face_lon" not in grid._ds or repopulate: # Construct the centroids if there are none stored if "face_x" not in grid._ds: - centroid_x, centroid_y, centroid_z = _build_face_centroids( + centroid_x, centroid_y, centroid_z = _construct_face_centroids( node_x, node_y, node_z, face_nodes, n_nodes_per_face ) @@ -173,9 +190,9 @@ def _populate_face_centroids(grid, repopulate=False): # If there are cartesian centroids already use those instead centroid_x, centroid_y, centroid_z = grid.face_x, grid.face_y, grid.face_z - # Convert from xyz to latlon + # Convert from xyz to latlon TODO centroid_lon, centroid_lat = _xyz_to_lonlat_deg( - centroid_x, centroid_y, centroid_z + centroid_x, centroid_y, centroid_z, normalize=False ) else: # Convert to xyz if there are latlon centroids already stored @@ -207,6 +224,28 @@ 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.""" + centroids = np.zeros((3, face_nodes.shape[0]), dtype=np.float64) + + for face_idx, n_max_nodes in enumerate(n_nodes_per_face): + # compute cartesian average + centroid_x = np.mean(node_x[face_nodes[face_idx, 0:n_max_nodes]]) + centroid_y = np.mean(node_y[face_nodes[face_idx, 0:n_max_nodes]]) + centroid_z = np.mean(node_z[face_nodes[face_idx, 0:n_max_nodes]]) + + # normalize coordinates + x_norm, y_norm, z_norm = _normalize_xyz(centroid_x, centroid_y, centroid_z) + + # store xyz + centroids[0, face_idx] = x_norm + centroids[1, face_idx] = y_norm + centroids[2, face_idx] = z_norm + + return centroids[0, :], centroids[1, :], centroids[2, :] + + 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 @@ -226,7 +265,7 @@ def _populate_edge_centroids(grid, repopulate=False): if "edge_lon" not in grid._ds or repopulate: # Construct the centroids if there are none stored if "edge_x" not in grid._ds: - centroid_x, centroid_y, centroid_z = _build_edge_centroids( + centroid_x, centroid_y, centroid_z = _construct_edge_centroids( node_x, node_y, node_z, edge_nodes_con ) @@ -236,7 +275,7 @@ def _populate_edge_centroids(grid, repopulate=False): # Convert from xyz to latlon centroid_lon, centroid_lat = _xyz_to_lonlat_deg( - centroid_x, centroid_y, centroid_z + centroid_x, centroid_y, centroid_z, normalize=False ) else: # Convert to xyz if there are latlon centroids already stored @@ -276,6 +315,28 @@ def _populate_edge_centroids(grid, repopulate=False): ) +def _construct_edge_centroids(node_x, node_y, node_z, edge_nodes_con): + """Constructs the xyz centroid coordinate for each edge using Cartesian + Averaging.""" + centroids = np.zeros((3, edge_nodes_con.shape[0]), dtype=np.float64) + + for edge_idx, connectivity in enumerate(edge_nodes_con): + # compute cartesian average + centroid_x = np.mean(node_x[connectivity[0:]]) + centroid_y = np.mean(node_y[connectivity[0:]]) + centroid_z = np.mean(node_z[connectivity[0:]]) + + # normalize coordinates + x_norm, y_nom, z_norm = _normalize_xyz(centroid_x, centroid_y, centroid_z) + + # store xyz + centroids[0, edge_idx] = x_norm + centroids[1, edge_idx] = y_nom + centroids[2, edge_idx] = z_norm + + return centroids[0, :], centroids[1, :], centroids[2, :] + + def _set_desired_longitude_range(ds): """Sets the longitude range to [-180, 180] for all longitude variables.""" From c1cd631551a9ca8bd544df37643978a9eee893cb Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 8 May 2024 20:30:52 -0500 Subject: [PATCH 08/14] fix tests (5) --- uxarray/grid/coordinates.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 6631ed284..43c1e7b9b 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -83,7 +83,7 @@ def _normalize_xyz( Union[np.ndarray, float], Union[np.ndarray, float], Union[np.ndarray, float] ]: # TODO - denom = np.linalg.norm(np.asarray([x, y, z]), ord=2) + denom = np.linalg.norm(np.asarray(np.array([x, y, z])), ord=2, axis=0) x_norm = x / denom y_norm = y / denom From edb1b794eab0003328e2cf98ac0d60eed87210de Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 8 May 2024 20:41:50 -0500 Subject: [PATCH 09/14] fix tests (6) --- test/test_coords.py | 31 ------------------------------- 1 file changed, 31 deletions(-) delete mode 100644 test/test_coords.py diff --git a/test/test_coords.py b/test/test_coords.py deleted file mode 100644 index 8e3fcc5ec..000000000 --- a/test/test_coords.py +++ /dev/null @@ -1,31 +0,0 @@ -import uxarray as ux - -from uxarray.constants import INT_FILL_VALUE -import numpy.testing as nt -import os - -import pytest - -from pathlib import Path - -current_path = Path(os.path.dirname(os.path.realpath(__file__))) - -GRID_PATHS = [ - current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc", - current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" -] - - - - -def test_lonlat_to_xyz(): - from uxarray.grid.veccoordinates import _populate_node_xyz - - for grid_path in GRID_PATHS: - uxgrid = ux.open_grid(grid_path) - - _populate_node_xyz(uxgrid) - - -def test_xyz_to_lonlat(): - pass From f133d0a9de9134c1e81a5469a220914a6fcc0a2a Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Thu, 9 May 2024 01:45:14 -0500 Subject: [PATCH 10/14] Update API reference --- docs/internal_api/index.rst | 14 ++-- docs/user_api/index.rst | 10 --- uxarray/grid/coordinates.py | 141 +++++++++++++++++------------------- 3 files changed, 75 insertions(+), 90 deletions(-) diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst index fc90d681d..b233e2413 100644 --- a/docs/internal_api/index.rst +++ b/docs/internal_api/index.rst @@ -147,15 +147,17 @@ Coordinates .. autosummary:: :toctree: generated/ - grid.coordinates._get_lonlat_from_xyz - grid.coordinates._get_xyz_from_lonlat - grid.coordinates._populate_cartesian_xyz_coord - grid.coordinates._populate_lonlat_coord + grid.coordinates._lonlat_rad_to_xyz + grid.coordinates._xyz_to_lonlat_rad + grid.coordinates._xyz_to_lonlat_deg + grid.coordinates._normalize_xyz + grid.coordinates._populate_node_latlon + grid.coordinates._populate_node_xyz grid.coordinates._populate_face_centroids - grid.coordinates._construct_face_centroids - grid.coordinates._set_desired_longitude_range grid.coordinates._populate_edge_centroids + grid.coordinates._construct_face_centroids grid.coordinates._construct_edge_centroids + grid.coordinates._set_desired_longitude_range Arcs diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst index 7f67456ff..ad51f3315 100644 --- a/docs/user_api/index.rst +++ b/docs/user_api/index.rst @@ -369,16 +369,6 @@ Connectivity grid.connectivity.close_face_nodes -Coordinates ------------ -.. autosummary:: - :toctree: generated/ - - grid.coordinates.node_lonlat_rad_to_xyz - grid.coordinates.node_xyz_to_lonlat_rad - grid.coordinates.normalize_in_place - - Arcs ---- .. autosummary:: diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 43c1e7b9b..9d3f34620 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -32,7 +32,28 @@ def _xyz_to_lonlat_rad( normalize: bool = True, ) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: """Converts Cartesian x, y, z coordinates in Spherical latitude and - longitude coordinates.""" + 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 + scalar: bool + Flag to compute a scalar conversion (i.e. x, y, z are each floats) + 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) @@ -41,8 +62,8 @@ def _xyz_to_lonlat_rad( y /= denom z /= denom - lon = np.arctan2(y, x) - lat = np.arcsin(z) + 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) @@ -63,10 +84,28 @@ def _xyz_to_lonlat_deg( z: Union[np.ndarray, float], normalize: bool = True, ) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: - lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z, normalize=normalize) + """Converts Cartesian x, y, z coordinates in Spherical latitude and + longitude coordinates in degrees. - # # set lon range to [-pi, pi] - # lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi + 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) @@ -82,8 +121,10 @@ def _normalize_xyz( ) -> tuple[ Union[np.ndarray, float], Union[np.ndarray, float], Union[np.ndarray, float] ]: - # TODO - denom = np.linalg.norm(np.asarray(np.array([x, y, z])), ord=2, axis=0) + """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 @@ -92,15 +133,12 @@ def _normalize_xyz( def _populate_node_latlon(grid) -> None: - """Docstring TODO.""" + """Populates the latitude and longitude 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 ) - # set longitude range to [-pi, pi] - # lon_rad = (lon_rad + np.pi) % (2 * np.pi) - np.pi - - # convert to degrees lon = np.rad2deg(lon_rad) lat = np.rad2deg(lat_rad) @@ -112,16 +150,9 @@ def _populate_node_latlon(grid) -> None: ) -def _populate_edge_latlon(grid) -> None: - pass - - -def _populate_face_latlon(grid) -> None: - pass - - def _populate_node_xyz(grid) -> None: - """Docstring TODO.""" + """Populates the Cartesiain node coordinates of a Grid (`node_x`, `node_y` + and `node_z`)""" node_lon_rad = np.deg2rad(grid.node_lon.values) node_lat_rad = np.deg2rad(grid.node_lat.values) @@ -138,28 +169,6 @@ def _populate_node_xyz(grid) -> None: ) -def _build_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face): - 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 _build_edge_centroids(node_x, node_y, node_z, edge_nodes): - centroid_x = np.mean(node_x[edge_nodes], axis=1) - centroid_y = np.mean(node_y[edge_nodes], axis=1) - centroid_z = np.mean(node_z[edge_nodes], axis=1) - - return _normalize_xyz(centroid_x, centroid_y, centroid_z) - - def _populate_face_centroids(grid, repopulate=False): """Finds the centroids of faces using cartesian averaging based off the vertices. The centroid is defined as the average of the x, y, z @@ -227,23 +236,18 @@ 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.""" - centroids = np.zeros((3, face_nodes.shape[0]), dtype=np.float64) - for face_idx, n_max_nodes in enumerate(n_nodes_per_face): - # compute cartesian average - centroid_x = np.mean(node_x[face_nodes[face_idx, 0:n_max_nodes]]) - centroid_y = np.mean(node_y[face_nodes[face_idx, 0:n_max_nodes]]) - centroid_z = np.mean(node_z[face_nodes[face_idx, 0:n_max_nodes]]) - - # normalize coordinates - x_norm, y_norm, z_norm = _normalize_xyz(centroid_x, centroid_y, centroid_z) + 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) - # store xyz - centroids[0, face_idx] = x_norm - centroids[1, face_idx] = y_norm - centroids[2, face_idx] = z_norm + 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 centroids[0, :], centroids[1, :], centroids[2, :] + return _normalize_xyz(centroid_x, centroid_y, centroid_z) def _populate_edge_centroids(grid, repopulate=False): @@ -315,26 +319,15 @@ def _populate_edge_centroids(grid, repopulate=False): ) -def _construct_edge_centroids(node_x, node_y, node_z, edge_nodes_con): +def _construct_edge_centroids(node_x, node_y, node_z, edge_node_conn): """Constructs the xyz centroid coordinate for each edge using Cartesian Averaging.""" - centroids = np.zeros((3, edge_nodes_con.shape[0]), dtype=np.float64) - - for edge_idx, connectivity in enumerate(edge_nodes_con): - # compute cartesian average - centroid_x = np.mean(node_x[connectivity[0:]]) - centroid_y = np.mean(node_y[connectivity[0:]]) - centroid_z = np.mean(node_z[connectivity[0:]]) - # normalize coordinates - x_norm, y_nom, z_norm = _normalize_xyz(centroid_x, centroid_y, centroid_z) + centroid_x = np.mean(node_x[edge_node_conn], axis=1) + centroid_y = np.mean(node_y[edge_node_conn], axis=1) + centroid_z = np.mean(node_z[edge_node_conn], axis=1) - # store xyz - centroids[0, edge_idx] = x_norm - centroids[1, edge_idx] = y_nom - centroids[2, edge_idx] = z_norm - - return centroids[0, :], centroids[1, :], centroids[2, :] + return _normalize_xyz(centroid_x, centroid_y, centroid_z) def _set_desired_longitude_range(ds): From 2bef0e94f25b603b47ccdc5e63a464c7930b417e Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Thu, 9 May 2024 13:45:14 -0500 Subject: [PATCH 11/14] generalize implementation --- uxarray/grid/arcs.py | 6 +++--- uxarray/grid/coordinates.py | 17 ++++++++++------- uxarray/grid/integrate.py | 5 +---- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index ec18195b4..381d02e41 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -60,12 +60,12 @@ def point_within_gca(pt, gca_cart, is_directed=False): Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. """ # Convert the cartesian coordinates to lonlat coordinates - pt_lonlat = np.array(_xyz_to_lonlat_rad(pt[0], pt[1], pt[2], scalar=True)) + pt_lonlat = np.array(_xyz_to_lonlat_rad(pt[0], pt[1], pt[2])) GCRv0_lonlat = np.array( - _xyz_to_lonlat_rad(gca_cart[0][0], gca_cart[0][1], gca_cart[0][2], scalar=True) + _xyz_to_lonlat_rad(gca_cart[0][0], gca_cart[0][1], gca_cart[0][2]) ) GCRv1_lonlat = np.array( - _xyz_to_lonlat_rad(gca_cart[1][0], gca_cart[1][1], gca_cart[1][2], scalar=True) + _xyz_to_lonlat_rad(gca_cart[1][0], gca_cart[1][1], gca_cart[1][2]) ) # Convert the list to np.float64 diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 9d3f34620..e1d77e2a3 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -28,7 +28,6 @@ def _xyz_to_lonlat_rad( x: Union[np.ndarray, float], y: Union[np.ndarray, float], z: Union[np.ndarray, float], - scalar: bool = False, normalize: bool = True, ) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: """Converts Cartesian x, y, z coordinates in Spherical latitude and @@ -68,13 +67,17 @@ def _xyz_to_lonlat_rad( # set longitude range to [0, pi] lon = np.mod(lon, 2 * np.pi) - if scalar: - if np.abs(z) > 1.0 - ERROR_TOLERANCE: - lat = np.sign(z) * np.pi / 2 - else: - # adjust z values near +- 1 - lat = np.where(np.abs(z) > 1.0 - ERROR_TOLERANCE, np.sign(z) * np.pi / 2, lat) + lat = np.where(np.abs(z) > 1.0 - ERROR_TOLERANCE, np.sign(z) * np.pi / 2, lat) + + # if lon.ndim == 0: + # return + # if scalar: + # if np.abs(z) > 1.0 - ERROR_TOLERANCE: + # lat = np.sign(z) * np.pi / 2 + # else: + # # adjust z values near +- 1 + # lat = np.where(np.abs(z) > 1.0 - ERROR_TOLERANCE, np.sign(z) * np.pi / 2, lat) return lon, lat diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index b18a25427..16cc13194 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -193,10 +193,7 @@ def _get_faces_constLat_intersection_info( if len(unique_intersections) == 2: # TODO: vectorize? unique_intersection_lonlat = np.array( - [ - _xyz_to_lonlat_rad(pt[0], pt[1], pt[2], scalar=True) - for pt in unique_intersections - ] + [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections] ) sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0) From 97a534ba55102b9f56f32691c96cfa3273581262 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Thu, 9 May 2024 13:52:28 -0500 Subject: [PATCH 12/14] update _lonlat_rad_to_xyz --- uxarray/grid/coordinates.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index e1d77e2a3..17d91b5f5 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -67,7 +67,10 @@ def _xyz_to_lonlat_rad( # set longitude range to [0, pi] lon = np.mod(lon, 2 * np.pi) - lat = np.where(np.abs(z) > 1.0 - ERROR_TOLERANCE, np.sign(z) * np.pi / 2, lat) + 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) # if lon.ndim == 0: # return From f94d67f7d07c6f1297d7a452fbe3dbee09711314 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Thu, 9 May 2024 14:41:20 -0500 Subject: [PATCH 13/14] update function signature --- uxarray/grid/coordinates.py | 23 ++++------------------- 1 file changed, 4 insertions(+), 19 deletions(-) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 17d91b5f5..7032d8d22 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -12,9 +12,7 @@ def _lonlat_rad_to_xyz( lon: Union[np.ndarray, float], lat: Union[np.ndarray, float], -) -> tuple[ - Union[np.ndarray, float], Union[np.ndarray, float], Union[np.ndarray, float] -]: +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Converts Spherical latitude and longitude coordinates into Cartesian x, y, z coordinates.""" x = np.cos(lon) * np.cos(lat) @@ -29,7 +27,7 @@ def _xyz_to_lonlat_rad( y: Union[np.ndarray, float], z: Union[np.ndarray, float], normalize: bool = True, -) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Converts Cartesian x, y, z coordinates in Spherical latitude and longitude coordinates in degrees. @@ -41,8 +39,6 @@ def _xyz_to_lonlat_rad( Cartesiain y coordinates z: Union[np.ndarray, float] Cartesian z coordinates - scalar: bool - Flag to compute a scalar conversion (i.e. x, y, z are each floats) normalize: bool Flag to select whether to normalize the coordinates @@ -72,15 +68,6 @@ def _xyz_to_lonlat_rad( lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) lon = np.where(z_mask, 0.0, lon) - # if lon.ndim == 0: - # return - - # if scalar: - # if np.abs(z) > 1.0 - ERROR_TOLERANCE: - # lat = np.sign(z) * np.pi / 2 - # else: - # # adjust z values near +- 1 - # lat = np.where(np.abs(z) > 1.0 - ERROR_TOLERANCE, np.sign(z) * np.pi / 2, lat) return lon, lat @@ -89,7 +76,7 @@ def _xyz_to_lonlat_deg( y: Union[np.ndarray, float], z: Union[np.ndarray, float], normalize: bool = True, -) -> tuple[Union[np.ndarray, float], Union[np.ndarray, float]]: +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Converts Cartesian x, y, z coordinates in Spherical latitude and longitude coordinates in degrees. @@ -124,9 +111,7 @@ def _normalize_xyz( x: Union[np.ndarray, float], y: Union[np.ndarray, float], z: Union[np.ndarray, float], -) -> tuple[ - Union[np.ndarray, float], Union[np.ndarray, float], 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 From 45d5841ee5f69fd8f9536267fbf6d6c93e9fc733 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 15 May 2024 16:17:50 -0500 Subject: [PATCH 14/14] vectorize conversion in tests, when possible --- test/test_geometry.py | 63 ++++++++++++++++++++++--------------------- 1 file changed, 33 insertions(+), 30 deletions(-) diff --git a/test/test_geometry.py b/test/test_geometry.py index 8b03d2430..b73e2d8b9 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -235,7 +235,7 @@ def _max_latitude_rad_iterative(self, gca_cart): max_lat = max(max_lat, w1_lonlat[1], w2_lonlat[1]) if np.abs(w2_lonlat[1] - w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[ - 1] == max_lat == w2_lonlat[1]: + 1] == max_lat == w2_lonlat[1]: max_section = [w1_new, w2_new] break if np.abs(max_lat - w1_lonlat[1]) <= ERROR_TOLERANCE: @@ -322,7 +322,7 @@ def _min_latitude_rad_iterative(self, gca_cart): min_lat = min(min_lat, w1_lonlat[1], w2_lonlat[1]) if np.abs(w2_lonlat[1] - w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[ - 1] == min_lat == w2_lonlat[1]: + 1] == min_lat == w2_lonlat[1]: min_section = [w1_new, w2_new] break if np.abs(min_lat - w1_lonlat[1]) <= ERROR_TOLERANCE: @@ -438,7 +438,6 @@ def test_insert_pt_in_empty_state(self): np.testing.assert_array_equal(result, expected) - class TestLatlonBoundsGCA(TestCase): def test_populate_bounds_normal(self): # Generate a normal face that is not crossing the antimeridian or the poles @@ -447,7 +446,9 @@ def test_populate_bounds_normal(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T + # vertices_cart = [_lonlat_rad_to_xyz(vertices_rad[0], vertices_rad[1])] lat_max = max(np.deg2rad(60.0), extreme_gca_latitude(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")) lat_min = min(np.deg2rad(10.0), @@ -476,7 +477,7 @@ def test_populate_bounds_antimeridian(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = max(np.deg2rad(60.0), extreme_gca_latitude(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")) lat_min = min(np.deg2rad(10.0), @@ -505,7 +506,7 @@ def test_populate_bounds_node_on_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = min(np.deg2rad(10.0), extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) @@ -533,7 +534,7 @@ def test_populate_bounds_edge_over_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = min(np.deg2rad(60.0), extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) @@ -561,7 +562,7 @@ def test_populate_bounds_pole_inside(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = min(np.deg2rad(60.0), extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")) @@ -582,6 +583,7 @@ def test_populate_bounds_pole_inside(self): bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat) nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + class TestLatlonBoundsLatLonFace(TestCase): def test_populate_bounds_normal(self): @@ -672,7 +674,7 @@ def test_populate_bounds_edge_over_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = np.deg2rad(60.0) lon_min = np.deg2rad(210.0) @@ -700,7 +702,7 @@ def test_populate_bounds_pole_inside(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = np.deg2rad(60.0) lon_min = 0 @@ -758,7 +760,7 @@ def test_populate_bounds_antimeridian(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.deg2rad(60.0) lat_min = np.deg2rad(10.0) lon_min = np.deg2rad(350.0) @@ -786,7 +788,7 @@ def test_populate_bounds_node_on_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = np.deg2rad(10.0) lon_min = np.deg2rad(10.0) @@ -814,7 +816,7 @@ def test_populate_bounds_edge_over_pole(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = np.deg2rad(60.0) lon_min = np.deg2rad(210.0) @@ -842,7 +844,7 @@ def test_populate_bounds_pole_inside(self): # Convert everything into radians vertices_rad = np.radians(vertices_lonlat) - vertices_cart = [_lonlat_rad_to_xyz(*v) for v in vertices_rad] + vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T lat_max = np.pi / 2 lat_min = np.deg2rad(60.0) lon_min = 0 @@ -863,6 +865,7 @@ def test_populate_bounds_pole_inside(self): is_GCA_list=[True, False, True, False]) nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + class TestLatlonBoundsMix(TestCase): def test_populate_bounds_GCA_mix(self): face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] @@ -873,11 +876,10 @@ def test_populate_bounds_GCA_mix(self): faces = [face_1, face_2, face_3, face_4] # Hand calculated bounds for the above faces in radians - expected_bounds = [[[0.17453293, 1.07370494],[0.17453293, 0.87266463]], - [[0.17453293, 1.10714872],[6.10865238, 0.87266463]], - [[1.04719755, 1.57079633],[3.66519143, 0.52359878]], - [[1.04719755,1.57079633],[0., 6.28318531]]] - + expected_bounds = [[[0.17453293, 1.07370494], [0.17453293, 0.87266463]], + [[0.17453293, 1.10714872], [6.10865238, 0.87266463]], + [[1.04719755, 1.57079633], [3.66519143, 0.52359878]], + [[1.04719755, 1.57079633], [0., 6.28318531]]] grid = ux.Grid.from_face_vertices(faces, latlon=True) face_bounds = grid.bounds.values @@ -892,14 +894,13 @@ def test_populate_bounds_LatlonFace_mix(self): faces = [face_1, face_2, face_3, face_4] - expected_bounds = [[[np.deg2rad(10.0), np.deg2rad(60.0)],[np.deg2rad(10.0), np.deg2rad(50.0)]], - [[np.deg2rad(10.0), np.deg2rad(60.0)],[np.deg2rad(350.0), np.deg2rad(50.0)]], - [[np.deg2rad(60.0), np.pi/2],[np.deg2rad(210.0), np.deg2rad(30.0)]], - [[np.deg2rad(60.0),np.pi/2],[0., 2*np.pi]]] - + expected_bounds = [[[np.deg2rad(10.0), np.deg2rad(60.0)], [np.deg2rad(10.0), np.deg2rad(50.0)]], + [[np.deg2rad(10.0), np.deg2rad(60.0)], [np.deg2rad(350.0), np.deg2rad(50.0)]], + [[np.deg2rad(60.0), np.pi / 2], [np.deg2rad(210.0), np.deg2rad(30.0)]], + [[np.deg2rad(60.0), np.pi / 2], [0., 2 * np.pi]]] grid = ux.Grid.from_face_vertices(faces, latlon=True) - bounds_xarray = _populate_bounds(grid,is_latlonface=True, return_array=True) + bounds_xarray = _populate_bounds(grid, is_latlonface=True, return_array=True) face_bounds = bounds_xarray.values for i in range(len(faces)): nt.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE) @@ -912,16 +913,18 @@ def test_populate_bounds_GCAList_mix(self): faces = [face_1, face_2, face_3, face_4] - expected_bounds = [[[np.deg2rad(10.0), np.deg2rad(60.0)],[np.deg2rad(10.0), np.deg2rad(50.0)]], - [[np.deg2rad(10.0), np.deg2rad(60.0)],[np.deg2rad(350.0), np.deg2rad(50.0)]], - [[np.deg2rad(60.0), np.pi/2],[np.deg2rad(210.0), np.deg2rad(30.0)]], - [[np.deg2rad(60.0),np.pi/2],[0., 2*np.pi]]] + expected_bounds = [[[np.deg2rad(10.0), np.deg2rad(60.0)], [np.deg2rad(10.0), np.deg2rad(50.0)]], + [[np.deg2rad(10.0), np.deg2rad(60.0)], [np.deg2rad(350.0), np.deg2rad(50.0)]], + [[np.deg2rad(60.0), np.pi / 2], [np.deg2rad(210.0), np.deg2rad(30.0)]], + [[np.deg2rad(60.0), np.pi / 2], [0., 2 * np.pi]]] grid = ux.Grid.from_face_vertices(faces, latlon=True) - bounds_xarray = _populate_bounds(grid,is_face_GCA_list=[[True, False, True, False]]*4, return_array=True) + bounds_xarray = _populate_bounds(grid, is_face_GCA_list=[[True, False, True, False]] * 4, return_array=True) face_bounds = bounds_xarray.values for i in range(len(faces)): nt.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE) + + class TestGeoDataFrame(TestCase): def test_to_gdf(self):