Skip to content

Commit

Permalink
fix from_cartesian polar cell unfold (#1301)
Browse files Browse the repository at this point in the history
* fix from_cartesian polar cell unfold

* add changelog new fragment
  • Loading branch information
bjlittle authored Jan 29, 2025
1 parent 58bb7ab commit ec65fb5
Show file tree
Hide file tree
Showing 3 changed files with 53 additions and 37 deletions.
2 changes: 2 additions & 0 deletions changelog/1301.bugfix.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
Fixed :func:`geovista.common.from_cartesian` to handle polar points with
no associated cell. (:user:`bjlittle`)
79 changes: 42 additions & 37 deletions src/geovista/common.py
Original file line number Diff line number Diff line change
Expand Up @@ -441,44 +441,49 @@ def from_cartesian(
lons[0] = lons[-1] = lons[1]
else:
pole_submesh = mesh.extract_points(pole_pids)
pole_pids = set(pole_pids)
# get the cids (cell-indices) of mesh cells with polar vertices
pole_cids = np.unique(pole_submesh[VTK_CELL_IDS])
for cid in pole_cids:
# get the pids (point-indices) of the polar cell points
# NOTE: pyvista 0.38.0: cell_point_ids(cid) -> get_cell(cid).point_ids
cell_pids = np.array(mesh.get_cell(cid).point_ids)
# unfold polar quad-cells
if len(cell_pids) == 4:
# identify the pids of the cell on the pole
cell_pole_pids = pole_pids.intersection(cell_pids)
# criterion of exactly two points from the quad-cell
# at the pole to unfold the polar points longitudes
if len(cell_pole_pids) == 2:
# compute the relative offset of the polar points
# within the polar cell connectivity
offset = sorted(
[np.where(cell_pids == pid)[0][0] for pid in cell_pole_pids]
)
if offset == [0, 1]:
lhs = cell_pids[offset]
rhs = cell_pids[[3, 2]]
elif offset == [1, 2]:
lhs = cell_pids[offset]
rhs = cell_pids[[0, 3]]
elif offset == [2, 3]:
lhs = cell_pids[offset]
rhs = cell_pids[[1, 0]]
elif offset == [0, 3]:
lhs = cell_pids[offset]
rhs = cell_pids[[1, 2]]
else:
emsg = (
"Failed to unfold a mesh polar quad-cell. Invalid "
"polar points connectivity detected."
if pole_submesh.n_cells:
pole_pids = set(pole_pids)
# get the cids (cell-indices) of mesh cells with polar vertices
pole_cids = np.unique(pole_submesh[VTK_CELL_IDS])
for cid in pole_cids:
# get the pids (point-indices) of the polar cell points
# NOTE:
# pyvista 0.38.0: cell_point_ids(cid) -> get_cell(cid).point_ids
cell_pids = np.array(mesh.get_cell(cid).point_ids)
# unfold polar quad-cells
if len(cell_pids) == 4:
# identify the pids of the cell on the pole
cell_pole_pids = pole_pids.intersection(cell_pids)
# criterion of exactly two points from the quad-cell
# at the pole to unfold the polar points longitudes
if len(cell_pole_pids) == 2:
# compute the relative offset of the polar points
# within the polar cell connectivity
offset = sorted(
[
np.where(cell_pids == pid)[0][0]
for pid in cell_pole_pids
]
)
raise ValueError(emsg)
lons[lhs] = lons[rhs]
if offset == [0, 1]:
lhs = cell_pids[offset]
rhs = cell_pids[[3, 2]]
elif offset == [1, 2]:
lhs = cell_pids[offset]
rhs = cell_pids[[0, 3]]
elif offset == [2, 3]:
lhs = cell_pids[offset]
rhs = cell_pids[[1, 0]]
elif offset == [0, 3]:
lhs = cell_pids[offset]
rhs = cell_pids[[1, 2]]
else:
emsg = (
"Failed to unfold a mesh polar quad-cell. Invalid "
"polar points connectivity detected."
)
raise ValueError(emsg)
lons[lhs] = lons[rhs]

if closed_interval:
if GV_REMESH_POINT_IDS in mesh.point_data:
Expand Down
9 changes: 9 additions & 0 deletions tests/common/test_from_cartesian.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
to_cartesian,
wrap,
)
from geovista.pantry.meshes import lfric_orog


def test_defaults(lam_uk_sample, lam_uk):
Expand Down Expand Up @@ -125,3 +126,11 @@ def test_lines_closed_interval(closed_interval, lonlat, pids):
mesh.lines = lines
result = from_cartesian(mesh, closed_interval=closed_interval)
assert np.all(np.isclose(result[:, :-1], lonlat)) == closed_interval


def test_convert_edges_pole_cell_unfold():
"""Test edges with no polar cell associated with polar point."""
mesh = lfric_orog()
_, edges = mesh.contour_banded(3)
result = from_cartesian(edges)
assert result.shape == (edges.n_points, 3)

0 comments on commit ec65fb5

Please sign in to comment.