Skip to content

Fix grounding line mask definition when creating landice meshes #935

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
trhille wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fix grounding line mask definition when creating landice meshes #935

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 20 additions & 5 deletions compass/landice/mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
from shutil import copyfile

import jigsawpy
import matplotlib.pyplot as plt
import mpas_tools.io
import numpy as np
import xarray
Expand Down Expand Up @@ -457,23 +458,37 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y,
[1, 1], [-1, 1], [1, -1], [-1, -1]])

ice_mask = thk > 0.0
grounded_mask = thk > (-1028.0 / 910.0 * topg)
grounded_mask = np.logical_and(thk > (-1028.0 / 910.0 * topg),
ice_mask)
float_mask = np.logical_and(thk < (-1028.0 / 910.0 * topg),
ice_mask)
margin_mask = np.zeros(sz, dtype='i')
grounding_line_mask = np.zeros(sz, dtype='i')

for n in neighbors:
not_ice_mask = np.logical_not(np.roll(ice_mask, n, axis=[0, 1]))
margin_mask = np.logical_or(margin_mask, not_ice_mask)

not_grounded_mask = np.logical_not(np.roll(grounded_mask,
n, axis=[0, 1]))
not_grounded_mask = np.logical_and(
Copy link
Collaborator Author

@trhille trhille Feb 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure this is the most straightforward way to define this. I have trouble thinking through the logic with these rolled arrays. If @matthewhoffman or @xylar has ideas about making this simpler, I'd be happy to hear them.

Copy link
Member

@matthewhoffman matthewhoffman Feb 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also find these things confusing. ChatGPT suggested this as an alternative:

shifted_grounded = np.roll(grounded_mask, n, axis=[0, 1])
shifted_float = np.roll(float_mask, n, axis=[0, 1])

not_grounded_mask = (~shifted_grounded) & shifted_float

Side question: would it be more accurate to rename not_grounded_mask to neighbor_not_grounded_mask?

np.logical_not(np.roll(grounded_mask,
n, axis=[0, 1])),
np.roll(float_mask, n, axis=[0, 1]))
grounding_line_mask = np.logical_or(grounding_line_mask,
not_grounded_mask)

# where ice exists and neighbors non-ice locations
margin_mask = np.logical_and(margin_mask, ice_mask)
# optional - plot mask
# plt.pcolor(margin_mask); plt.show()
# where grounded ice exists and neighbors floating ice
grounding_line_mask = np.logical_and(grounding_line_mask, grounded_mask)

fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(6, 3))
margin_plot = ax[0].pcolor(margin_mask)
gl_plot = ax[1].pcolor(grounding_line_mask) # noqa F841
ax[0].set_title("margin mask")
ax[1].set_title("grounding line mask")
plt.colorbar(margin_plot, ax=[ax[0], ax[1]], shrink=0.7)
[ax.set_aspect('equal') for ax in ax]
fig.savefig("masks.png", dpi=400)

# Calculate dist to margin and grounding line
[XPOS, YPOS] = np.meshgrid(x, y)
Expand Down
Loading