Merge pull request #322 from trhille/generalize_dist_to_edge_and_GL

Create function to calculate distance to ice edge and grounding line

Add new function compass.landice.mesh.get_dist_to_edge_and_GL() that calculates the distance from each grid point in gridded dataset to the ice margin and grounding line, which is then passed to compass.landice.mesh.set_cell_width() as input to the cell spacing functions. Add use case for Humboldt.

compass/landice/mesh.py

@@ -1,6 +1,6 @@
 import numpy as np
 import jigsawpy
-
+import time
 
 def gridded_flood_fill(field):
     """
@@ -18,18 +18,18 @@ def gridded_flood_fill(field):
 
     Returns
     -------
-    floodMask : numpy.ndarray
-        Mask calculated by the flood fill routine,
+    flood_mask : numpy.ndarray
+        _mask calculated by the flood fill routine,
         where cells connected to the ice sheet (or main feature)
         are 1 and everything else is 0.
     """
 
     sz = field.shape
-    searchedMask = np.zeros(sz)
-    floodMask = np.zeros(sz)
+    searched_mask = np.zeros(sz)
+    flood_mask = np.zeros(sz)
     iStart = sz[0] // 2
     jStart = sz[1] // 2
-    floodMask[iStart, jStart] = 1
+    flood_mask[iStart, jStart] = 1
 
     neighbors = np.array([[1, 0], [-1, 0], [0, 1], [0, -1]])
 
@@ -48,19 +48,19 @@ def gridded_flood_fill(field):
                 ii = i + n[0]
                 jj = j + n[1]  # subscripts to neighbor
                 # only consider unsearched neighbors
-                if searchedMask[ii, jj] == 0:
-                    searchedMask[ii, jj] = 1  # mark as searched
+                if searched_mask[ii, jj] == 0:
+                    searched_mask[ii, jj] = 1  # mark as searched
 
                     if field[ii, jj] > 0.0:
-                        floodMask[ii, jj] = 1  # mark as ice
+                        flood_mask[ii, jj] = 1  # mark as ice
                         # add to list of newly found  cells
                         newSearchList = np.append(newSearchList,
                                                   np.ravel_multi_index(
                                                       [[ii], [jj]], sz,
                                                       order='F')[0])
         lastSearchList = newSearchList
 
-    return floodMask
+    return flood_mask
 
 
 def set_rectangular_geom_points_and_edges(xmin, xmax, ymin, ymax):
@@ -204,3 +204,109 @@ def set_cell_width(self, section, thk, vx=None, vy=None,
                    dist_to_edge > (10. * cull_distance))] = max_spac
 
     return cell_width
+
+
+def get_dist_to_edge_and_GL(self, thk, topg, x, y, window_size=1.e5):
+    """
+    Calculate distance from each point to ice edge and grounding line,
+    to be used in mesh density functions in
+    :py:func:`compass.landice.mesh.set_cell_width()`. In future development,
+    this should be updated to use a faster package such as `scikit-fmm`.
+
+    Parameters
+    ----------
+    thk : numpy.ndarray
+        Ice thickness field from gridded dataset,
+        usually after trimming to flood fill mask
+    topg : numpy.ndarray
+        Bed topography field from gridded dataset
+    x : numpy.ndarray
+        x coordinates from gridded dataset
+    y : numpy.ndarray
+        y coordinates from gridded dataset
+    window_size : int or float
+        Size (in meters) of a search 'box' (one-directional) to use
+        to calculate the distance from each cell to the ice margin.
+        Bigger number makes search slower, but if too small, the transition
+        zone could get truncated.
+
+    Returns
+    -------
+    dist_to_edge : numpy.ndarray
+        Distance from each cell to the ice edge
+    dist_to_grounding_line : numpy.ndarray
+        Distance from each cell to the grounding line
+    """
+    logger = self.logger
+    tic = time.time()
+
+    dx = x[1] - x[0]  # assumed constant and equal in x and y
+    nx = len(x)
+    ny = len(y)
+    sz = thk.shape
+
+    # Create masks to define ice edge and grounding line
+    neighbors = np.array([[1, 0], [-1, 0], [0, 1], [0, -1],
+                          [1, 1], [-1, 1], [1, -1], [-1, -1]])
+
+    ice_mask = thk > 0.0
+    grounded_mask = thk > (-1028.0 / 910.0 * topg)
+    floating_mask = np.logical_and(thk < (-1028.0 /
+                                          910.0 * topg), thk > 0.0)
+    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]))
+        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()
+
+    # Calculate dist to margin and grounding line
+    [XPOS, YPOS] = np.meshgrid(x, y)
+    dist_to_edge = np.zeros(sz)
+    dist_to_grounding_line = np.zeros(sz)
+
+    d = int(np.ceil(window_size / dx))
+    rng = np.arange(-1*d, d, dtype='i')
+    max_dist = float(d) * dx
+
+    # just look over areas with ice
+    # ind = np.where(np.ravel(thk, order='F') > 0)[0]
+    ind = np.where(np.ravel(thk, order='F') >= 0)[0]  # do it everywhere
+    for iii in range(len(ind)):
+        [i, j] = np.unravel_index(ind[iii], sz, order='F')
+
+        irng = i + rng
+        jrng = j + rng
+
+        # only keep indices in the grid
+        irng = irng[np.nonzero(np.logical_and(irng >= 0, irng < ny))]
+        jrng = jrng[np.nonzero(np.logical_and(jrng >= 0, jrng < nx))]
+
+        dist_to_here = ((XPOS[np.ix_(irng, jrng)] - x[j])**2 +
+                        (YPOS[np.ix_(irng, jrng)] - y[i])**2)**0.5
+
+        dist_to_here_edge = dist_to_here.copy()
+        dist_to_here_grounding_line = dist_to_here.copy()
+
+        dist_to_here_edge[margin_mask[np.ix_(irng, jrng)] == 0] = max_dist
+        dist_to_here_grounding_line[grounding_line_mask
+                                    [np.ix_(irng, jrng)] == 0] = max_dist
+
+        dist_to_edge[i, j] = dist_to_here_edge.min()
+        dist_to_grounding_line[i, j] = dist_to_here_grounding_line.min()
+
+    toc = time.time()
+    logger.info('compass.landice.mesh.get_dist_to_edge_and_GL() took {:0.2f} '
+                'seconds'.format(toc - tic))
+
+    return dist_to_edge, dist_to_grounding_line

compass/landice/tests/humboldt/mesh.py

@@ -13,7 +13,7 @@
 from compass.model import make_graph_file
 from compass.landice.mesh import gridded_flood_fill, \
                                  set_rectangular_geom_points_and_edges, \
-                                 set_cell_width
+                                 set_cell_width, get_dist_to_edge_and_GL
 
 
 class Mesh(Step):
@@ -189,12 +189,6 @@ def build_cell_width(self):
         vx = f.variables['vx'][0, :, :]
         vy = f.variables['vy'][0, :, :]
 
-        dx = x1[1] - x1[0]  # assumed constant and equal in x and y
-        nx = len(x1)
-        ny = len(y1)
-
-        sz = thk.shape
-
         # Define extent of region to mesh.
         # These coords are specific to the Humboldt Glacier mesh.
         xx0 = -630000
@@ -210,68 +204,17 @@ def build_cell_width(self):
         vx[floodMask == 0] = 0.0
         vy[floodMask == 0] = 0.0
 
-        # make masks -------------------
-        neighbors = np.array([[1, 0], [-1, 0], [0, 1], [0, -1],
-                              [1, 1], [-1, 1], [1, -1], [-1, -1]])
-
-        iceMask = thk > 0.0
-        # groundedMask = thk > (-1028.0 / 910.0 * topg)
-        # floatingMask = np.logical_and(thk < (-1028.0 /
-        #                               910.0 * topg), thk > 0.0)
-        marginMask = np.zeros(sz, dtype='i')
-        for n in neighbors:
-            marginMask = np.logical_or(marginMask,
-                                       np.logical_not(
-                                           np.roll(iceMask, n, axis=[0, 1])))
-        # where ice exists and neighbors non-ice locations
-        marginMask = np.logical_and(marginMask, iceMask)
-        # optional - plot mask
-        # plt.pcolor(marginMask); plt.show()
-
-        # calc dist to margin -------------------
-        [XPOS, YPOS] = np.meshgrid(x1, y1)
-        distToEdge = np.zeros(sz)
-
-        # -- KEY PARAMETER: how big of a search 'box' (one-directional) to use
-        # to calculate the distance from each cell to the ice margin.
-        # Bigger number makes search slower, but if too small, the transition
-        # zone could get truncated. Could automatically set this from maxDist
-        # variables used in next section. Currently, this is only used to
-        # determine mesh spacing in ice-free areas in order to keep mesh
-        # density low in areas that will be culled.
-        windowSize = 100.0e3
-        # ---
-
-        d = int(np.ceil(windowSize / dx))
-        # logger.info(windowSize, d)
-        rng = np.arange(-1*d, d, dtype='i')
-        maxdist = float(d) * dx
-
-        # just look over areas with ice
-        # ind = np.where(np.ravel(thk, order='F') > 0)[0]
-        ind = np.where(np.ravel(thk, order='F') >= 0)[0]  # do it everywhere
-        for iii in range(len(ind)):
-            [i, j] = np.unravel_index(ind[iii], sz, order='F')
-
-            irng = i + rng
-            jrng = j + rng
-
-            # only keep indices in the grid
-            irng = irng[np.nonzero(np.logical_and(irng >= 0, irng < ny))]
-            jrng = jrng[np.nonzero(np.logical_and(jrng >= 0, jrng < nx))]
-
-            dist2Here = ((XPOS[np.ix_(irng, jrng)] - x1[j])**2 +
-                         (YPOS[np.ix_(irng, jrng)] - y1[i])**2)**0.5
-            dist2Here[marginMask[np.ix_(irng, jrng)] == 0] = maxdist
-            distToEdge[i, j] = dist2Here.min()
+        # Calculate distance from each grid point to ice edge
+        # and grounding line, for use in cell spacing functions.
+        distToEdge, distToGL = get_dist_to_edge_and_GL(self, thk, topg, x1,
+                                                       y1, window_size=1.e5)
         # optional - plot distance calculation
         # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show()
 
         # Set cell widths based on mesh parameters set in config file
         cell_width = set_cell_width(self, section='humboldt', thk=thk,
                                     vx=vx, vy=vy, dist_to_edge=distToEdge,
                                     dist_to_grounding_line=None)
-
         # plt.pcolor(cell_width); plt.colorbar(); plt.show()
 
         return (cell_width.astype('float64'), x1.astype('float64'),