Skip to content

Commit

Permalink
Optimize the get_cartesian_face_edge_nodes and `get_lonlat_rad_face…
Browse files Browse the repository at this point in the history 
…_edge_nodes` (#799)

* Initial Draft Algorithm

* Revert "Initial Draft Algorithm"

This reverts commit 99573b8.

* Initial commit

* Revert "Initial commit"

This reverts commit 1ba6880.

* initial commit

* Finished implementation

* Fixed faces with holes

* Fix precommit

* test

* add examples to doc strings

* initial commit for vectorizing _get_cartesian_face_edge_nodes and _get_lonlat_rad_face_edge_nodes

* add examples, modify test_helpers to test the new implementation(cur fail cartesian with fill)

* Fix precommit

* continue testcase

* fixed input with INT_FILL_VALUE in face_node_conn

* fix test case

* more tests

* rename helper

* rename

* update usage of _get_cartesian_face_edge_nodes and _get_lonlat_rad_face_edge_nodes

* update helper index

---------

Co-authored-by: Philip Chmielowiec <[email protected]>
Co-authored-by: amberchen122 <[email protected]>
Co-authored-by: amberchen122 <amber@DESKTOP-SS0HRJ9>
  • Loading branch information
@hongyuchen1030 @philipc2 @amberchen122
4 people authored Jun 30, 2024
1 parent 6719a54 commit 3b27d3f
Show file tree
Hide file tree
Showing 6 changed files with 872 additions and 183 deletions.
1 change: 1 addition & 0 deletions docs/internal_api/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,6 +176,7 @@ Utils

grid.utils._newton_raphson_solver_for_gca_constLat
grid.utils._inv_jacobian
grid.utils._swap_first_fill_value_with_last
grid.utils._get_cartesiain_face_edge_nodes
grid.utils._get_lonlat_rad_face_edge_nodes

Expand Down
423 changes: 393 additions & 30 deletions test/test_geometry.py
View file

Large diffs are not rendered by default.

7 changes: 7 additions & 0 deletions test/test_grid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -913,6 +913,7 @@ def test_from_face_vertices(self):


class TestLatlonBounds(TestCase):
gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
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]]
face_2 = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
Expand All @@ -932,3 +933,9 @@ def test_populate_bounds_GCA_mix(self):
bounds_xarray = grid.bounds
face_bounds = bounds_xarray.values
nt.assert_allclose(grid.bounds.values, expected_bounds, atol=ERROR_TOLERANCE)

def test_populate_bounds_MPAS(self):
xrds = xr.open_dataset(self.gridfile_mpas)
uxgrid = ux.Grid.from_dataset(xrds, use_dual=True)
bounds_xarray = uxgrid.bounds
pass
304 changes: 248 additions & 56 deletions test/test_helpers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,111 +241,303 @@ def test_angle_of_2_vectors(self):

class TestFaceEdgeConnectivityHelper(TestCase):

def test_get_cartesian_face_edge_nodes(self):
# Load the dataset

uxds = ux.open_dataset(gridfile_geoflowsmall_grid,
gridfile_geoflowsmall_var)
def test_get_cartesian_face_edge_nodes_pipeline(self):
# Create the vertices for the grid, based around the North Pole
vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]

# Initialize array
# Normalize the vertices
vertices = [x / np.linalg.norm(x) for x in vertices]

face_edges_connectivity_cartesian = []
# Construct the grid from the vertices
grid = ux.Grid.from_face_vertices(vertices, latlon=False)

# Get the connectivity
# Extract the necessary grid data
face_node_conn = grid.face_node_connectivity.values
n_nodes_per_face = np.array([len(face) for face in face_node_conn])
n_face = len(face_node_conn)

for i in range(len(uxds.uxgrid.face_node_connectivity)):
face_edges_connectivity_cartesian.append(
_get_cartesian_face_edge_nodes(
uxds.uxgrid.face_node_connectivity.values[i],
uxds.uxgrid.face_edge_connectivity.values[i],
uxds.uxgrid.edge_node_connectivity.values,
uxds.uxgrid.node_x.values, uxds.uxgrid.node_y.values,
uxds.uxgrid.node_z.values))
n_max_face_edges = max(n_nodes_per_face)
node_x = grid.node_x.values
node_y = grid.node_y.values
node_z = grid.node_z.values

# Stack the arrays to get the desired (3,3) array
# Call the function to test
face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
)

face_edges_connectivity_cartesian = np.vstack(
face_edges_connectivity_cartesian)
# Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
result = ux.grid.geometry._pole_point_inside_polygon(
'North', face_edges_connectivity_cartesian[0]
)

assert (face_edges_connectivity_cartesian.ndim == 3)
# Assert that the result is True
self.assertTrue(result)

def test_get_cartesian_face_edge_nodes_pipeline(self):
def test_get_cartesian_face_edge_nodes_filled_value(self):
# Create the vertices for the grid, based around the North Pole
vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]

vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5],
[0.5, -0.5, 0.5]]

#Normalize the vertices
# Normalize the vertices
vertices = [x / np.linalg.norm(x) for x in vertices]
vertices.append([INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE])

# Construct the grid from the vertices
grid = ux.Grid.from_face_vertices(vertices, latlon=False)

# Extract the necessary grid data
face_node_conn = grid.face_node_connectivity.values
n_nodes_per_face = np.array([len(face) for face in face_node_conn])
n_face = len(face_node_conn)
n_max_face_edges = max(n_nodes_per_face)
node_x = grid.node_x.values
node_y = grid.node_y.values
node_z = grid.node_z.values

# Call the function to test
face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
grid.face_node_connectivity.values[0],
grid.face_edge_connectivity.values[0],
grid.edge_node_connectivity.values, grid.node_x.values,
grid.node_y.values, grid.node_z.values)
face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
)

# Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
result = ux.grid.geometry._pole_point_inside_polygon(
'North', face_edges_connectivity_cartesian)
'North', face_edges_connectivity_cartesian[0]
)

# Assert that the result is True

self.assertTrue(result)

def test_get_cartesian_face_edge_nodes_filled_value(self):
# Create the vertices for the grid, based around the North Pole
def test_get_cartesian_face_edge_nodes_filled_value2(self):
# The face vertices order in counter-clockwise
# face_conn = [[0,1,2],[1,3,4,2]]

#Each vertex is a 2D vector represent the longitude and latitude in degree. Call the node_lonlat_to_xyz to convert it to 3D vector
v0_deg = [10,10]
v1_deg = [15,15]
v2_deg = [5,15]
v3_deg = [15,45]
v4_deg = [5,45]

# First convert them into radians
v0_rad = np.deg2rad(v0_deg)
v1_rad = np.deg2rad(v1_deg)
v2_rad = np.deg2rad(v2_deg)
v3_rad = np.deg2rad(v3_deg)
v4_rad = np.deg2rad(v4_deg)

# It should look like following when passing in the _get_cartesian_face_edge_nodes
# [[v0_cart,v1_cart,v2_cart, [Fill_Value,Fill_Value,Fill_Value]],[v1_cart,v3_cart,v4_cart,v2_cart]]
v0_cart = _lonlat_rad_to_xyz(v0_rad[0],v0_rad[1])
v1_cart = _lonlat_rad_to_xyz(v1_rad[0],v1_rad[1])
v2_cart = _lonlat_rad_to_xyz(v2_rad[0],v2_rad[1])
v3_cart = _lonlat_rad_to_xyz(v3_rad[0],v3_rad[1])
v4_cart = _lonlat_rad_to_xyz(v4_rad[0],v4_rad[1])

face_node_conn = np.array([[0, 1, 2, INT_FILL_VALUE],[1, 3, 4, 2]])
n_face = 2
n_max_face_edges = 4
n_nodes_per_face = np.array([len(face) for face in face_node_conn])
node_x = np.array([v0_cart[0],v1_cart[0],v2_cart[0],v3_cart[0],v4_cart[0]])
node_y = np.array([v0_cart[1],v1_cart[1],v2_cart[1],v3_cart[1],v4_cart[1]])
node_z = np.array([v0_cart[2],v1_cart[2],v2_cart[2],v3_cart[2],v4_cart[2]])

# call the function to test
face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
)

# Define correct result
correct_result = np.array([
[
[
[v0_cart[0], v0_cart[1], v0_cart[2]],
[v1_cart[0], v1_cart[1], v1_cart[2]]
],
[
[v1_cart[0], v1_cart[1], v1_cart[2]],
[v2_cart[0], v2_cart[1], v2_cart[2]]
],
[
[v2_cart[0], v2_cart[1], v2_cart[2]],
[v0_cart[0], v0_cart[1], v0_cart[2]]
],
[
[INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE],
[INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE]
]
],
[
[
[v1_cart[0], v1_cart[1], v1_cart[2]],
[v3_cart[0], v3_cart[1], v3_cart[2]]
],
[
[v3_cart[0], v3_cart[1], v3_cart[2]],
[v4_cart[0], v4_cart[1], v4_cart[2]]
],
[
[v4_cart[0], v4_cart[1], v4_cart[2]],
[v2_cart[0], v2_cart[1], v2_cart[2]]
],
[
[v2_cart[0], v2_cart[1], v2_cart[2]],
[v1_cart[0], v1_cart[1], v1_cart[2]]
]
]
])


vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5],
[0.5, -0.5, 0.5]]
# Assert that the result is correct
self.assertEqual(face_edges_connectivity_cartesian.shape, correct_result.shape)

#Normalize the vertices

def test_get_lonlat_face_edge_nodes_pipeline(self):
# Create the vertices for the grid, based around the North Pole
vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]

# Normalize the vertices
vertices = [x / np.linalg.norm(x) for x in vertices]
vertices.append([INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE])

# Construct the grid from the vertices
grid = ux.Grid.from_face_vertices(vertices, latlon=False)
face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
grid.face_node_connectivity.values[0],
grid.face_edge_connectivity.values[0],
grid.edge_node_connectivity.values, grid.node_x.values,
grid.node_y.values, grid.node_z.values)

# Extract the necessary grid data
face_node_conn = grid.face_node_connectivity.values
n_nodes_per_face = np.array([len(face) for face in face_node_conn])
n_face = len(face_node_conn)
n_max_face_edges = max(n_nodes_per_face)
node_lon = grid.node_lon.values
node_lat = grid.node_lat.values

# Call the function to test
face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
)

# Convert the first face's edges to Cartesian coordinates
face_edges_connectivity_lonlat = face_edges_connectivity_lonlat[0]
face_edges_connectivity_cartesian = []
for edge in face_edges_connectivity_lonlat:
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
result = ux.grid.geometry._pole_point_inside_polygon(
'North', face_edges_connectivity_cartesian)
'North', np.array(face_edges_connectivity_cartesian)
)

# Assert that the result is True
self.assertTrue(result)

def test_get_lonlat_face_edge_nodes_pipeline(self):
def test_get_lonlat_face_edge_nodes_filled_value(self):
# Create the vertices for the grid, based around the North Pole
vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]

vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5],
[0.5, -0.5, 0.5]]

#Normalize the vertices
# Normalize the vertices
vertices = [x / np.linalg.norm(x) for x in vertices]
vertices.append([INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE])

# Construct the grid from the vertices
grid = ux.Grid.from_face_vertices(vertices, latlon=False)

# Extract the necessary grid data
face_node_conn = grid.face_node_connectivity.values
n_nodes_per_face = np.array([len(face) for face in face_node_conn])
n_face = len(face_node_conn)
n_max_face_edges = max(n_nodes_per_face)
node_lon = grid.node_lon.values
node_lat = grid.node_lat.values

# Call the function to test
face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
grid.face_node_connectivity.values[0],
grid.face_edge_connectivity.values[0],
grid.edge_node_connectivity.values, grid.node_lon.values,
grid.node_lat.values)
face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
)

# Convert all the values into cartesian coordinates
# Convert the first face's edges to Cartesian coordinates
face_edges_connectivity_lonlat = face_edges_connectivity_lonlat[0]
face_edges_connectivity_cartesian = []
for i in range(len(face_edges_connectivity_lonlat)):
edge = face_edges_connectivity_lonlat[i]
for edge in face_edges_connectivity_lonlat:
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
result = ux.grid.geometry._pole_point_inside_polygon(
'North', np.array(face_edges_connectivity_cartesian))
'North', np.array(face_edges_connectivity_cartesian)
)

# Assert that the result is True
self.assertTrue(result)


def test_get_lonlat_face_edge_nodes_filled_value2(self):
# The face vertices order in counter-clockwise
# face_conn = [[0,1,2],[1,3,4,2]]

#Each vertex is a 2D vector represent the longitude and latitude in degree. Call the node_lonlat_to_xyz to convert it to 3D vector
v0_deg = [10,10]
v1_deg = [15,15]
v2_deg = [5,15]
v3_deg = [15,45]
v4_deg = [5,45]

# First convert them into radians
v0_rad = np.deg2rad(v0_deg)
v1_rad = np.deg2rad(v1_deg)
v2_rad = np.deg2rad(v2_deg)
v3_rad = np.deg2rad(v3_deg)
v4_rad = np.deg2rad(v4_deg)

face_node_conn = np.array([[0, 1, 2, INT_FILL_VALUE],[1, 3, 4, 2]])
n_face = 2
n_max_face_edges = 4
n_nodes_per_face = np.array([len(face) for face in face_node_conn])
node_lon = np.array([v0_rad[0],v1_rad[0],v2_rad[0],v3_rad[0],v4_rad[0]])
node_lat = np.array([v0_rad[1],v1_rad[1],v2_rad[1],v3_rad[1],v4_rad[1]])

# call the function to test
face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
)

# Define correct result
correct_result = np.array([
[
[
[v0_rad[0], v0_rad[1]],
[v1_rad[0], v1_rad[1]]
],
[
[v1_rad[0], v1_rad[1]],
[v2_rad[0], v2_rad[1]]
],
[
[v2_rad[0], v2_rad[1]],
[v0_rad[0], v0_rad[1]]
],
[
[INT_FILL_VALUE, INT_FILL_VALUE],
[INT_FILL_VALUE, INT_FILL_VALUE]
]
],
[
[
[v1_rad[0], v1_rad[1]],
[v3_rad[0], v3_rad[1]]
],
[
[v3_rad[0], v3_rad[1]],
[v4_rad[0], v4_rad[1]]
],
[
[v4_rad[0], v4_rad[1]],
[v2_rad[0], v2_rad[1]]
],
[
[v2_rad[0], v2_rad[1]],
[v1_rad[0], v1_rad[1]]
]
]
])

# Assert that the result is correct
self.assertEqual(face_edges_connectivity_lonlat.shape, correct_result.shape)
Loading

0 comments on commit 3b27d3f

Please sign in to comment.