Reduce number of empty cells in extended high-resolution Greenland Ice Sheet mesh

compass/landice/mesh.py

@@ -224,8 +224,9 @@ def clip_mesh_to_bounding_box(mask_ds, base_ds, bounding_box):
     return mask_ds
 
 
-def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
+def set_cell_width(self, section_name, thk, bed, vx=None, vy=None,
                    dist_to_edge=None, dist_to_grounding_line=None,
+                   dist_to_coast=None,
                    flood_fill_iStart=None, flood_fill_jStart=None):
     """
     Set cell widths based on settings in config file to pass to
@@ -238,9 +239,10 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
         following options to be set in the given config section:
         ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``,
         ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``,
-        ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``,
-        ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``,
-        ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``.
+        ``high_dist``, ``low_dist``, ``high_dist_coast``, ``low_dist_coast``,
+        ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``,
+        ``cull_distance``, ``use_speed``, ``use_dist_to_edge``,
+        ``use_dist_to_grounding_line``, ``use_dist_to_coast``, and ``use_bed``.
         See the Land-Ice Framework section of the Users or Developers guide
         for more information about these options and their uses.
 
@@ -272,6 +274,11 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
         function.  Can be set to ``None`` if
         ``use_dist_to_grounding_line == False`` in config file.
 
+    dist_to_coast : numpy.ndarray, optional
+        Distance from each cell to coast, calculated in separate
+        function.  Can be set to ``None`` if
+        ``use_dist_to_coast == False`` in config file.
+
     flood_fill_iStart : int, optional
         x-index location to start flood-fill when using bed topography
 
@@ -291,21 +298,19 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
     # Get config inputs for cell spacing functions
     min_spac = section.getfloat('min_spac')
     max_spac = section.getfloat('max_spac')
-    high_log_speed = section.getfloat('high_log_speed')
-    low_log_speed = section.getfloat('low_log_speed')
-    high_dist = section.getfloat('high_dist')
-    low_dist = section.getfloat('low_dist')
-    high_dist_bed = section.getfloat('high_dist_bed')
-    low_dist_bed = section.getfloat('low_dist_bed')
-    low_bed = section.getfloat('low_bed')
-    high_bed = section.getfloat('high_bed')
-
     # convert km to m
     cull_distance = section.getfloat('cull_distance') * 1.e3
 
+    land_mask = np.logical_and(thk == 0.0, bed >= 0.)
+    ocean_mask = np.logical_and(thk == 0.0, bed < 0.)
+
     # Cell spacing function based on union of masks
     if section.get('use_bed') == 'True':
         logger.info('Using bed elevation for spacing.')
+        high_dist_bed = section.getfloat('high_dist_bed')
+        low_dist_bed = section.getfloat('low_dist_bed')
+        low_bed = section.getfloat('low_bed')
+        high_bed = section.getfloat('high_bed')
         if flood_fill_iStart is not None and flood_fill_jStart is not None:
             logger.info('calling gridded_flood_fill to find ' +
                         'bedTopography <= low_bed connected to the ocean.')
@@ -355,6 +360,8 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
     # Make cell spacing function mapping from log speed to cell spacing
     if section.get('use_speed') == 'True':
         logger.info('Using speed for cell spacing')
+        high_log_speed = section.getfloat('high_log_speed')
+        low_log_speed = section.getfloat('low_log_speed')
         speed = (vx ** 2 + vy ** 2) ** 0.5
         lspd = np.log10(speed)
         spacing_speed = np.interp(lspd, [low_log_speed, high_log_speed],
@@ -372,33 +379,53 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
                     f'dataset with missing velocity values. Setting '
                     f'velocity-based spacing to maximum value.')
 
-        spacing_speed[thk == 0.0] = min_spac
     else:
         spacing_speed = max_spac * np.ones_like(thk)
 
     # Make cell spacing function mapping from distance to ice edge
     if section.get('use_dist_to_edge') == 'True':
         logger.info('Using distance to ice edge for cell spacing')
+        high_dist = section.getfloat('high_dist')
+        low_dist = section.getfloat('low_dist')
         spacing_edge = np.interp(dist_to_edge, [low_dist, high_dist],
                                  [min_spac, max_spac], left=min_spac,
                                  right=max_spac)
-        spacing_edge[thk == 0.0] = min_spac
     else:
         spacing_edge = max_spac * np.ones_like(thk)
 
     # Make cell spacing function mapping from distance to grounding line
     if section.get('use_dist_to_grounding_line') == 'True':
         logger.info('Using distance to grounding line for cell spacing')
+        high_dist = section.getfloat('high_dist')
+        low_dist = section.getfloat('low_dist')
         spacing_gl = np.interp(dist_to_grounding_line, [low_dist, high_dist],
                                [min_spac, max_spac], left=min_spac,
                                right=max_spac)
-        spacing_gl[thk == 0.0] = min_spac
     else:
         spacing_gl = max_spac * np.ones_like(thk)
 
+    # Make cell spacing function mapping from distance to coast
+    if section.get('use_dist_to_coast') == 'True':
+        logger.info('Using distance to coast for cell spacing')
+        high_dist_coast = section.getfloat('high_dist_coast')
+        low_dist_coast = section.getfloat('low_dist_coast')
+        spacing_coast = np.interp(dist_to_coast,
+                                  [low_dist_coast, high_dist_coast],
+                                  [min_spac, max_spac], left=min_spac,
+                                  right=max_spac)
+        # distance from coast is only used to set spacing in ocean
+        spacing_coast[bed > 0.] = max_spac
+    else:
+        spacing_coast = max_spac * np.ones_like(thk)
+        # If not using distance from coast, use finest spacing
+        # in the ocean as well as ice-free land.
+        spacing_coast[land_mask] = min_spac
+        spacing_coast[ocean_mask] = min_spac
+
     # Merge cell spacing methods
     cell_width = max_spac * np.ones_like(thk)
-    for width in [spacing_bed, spacing_speed, spacing_edge, spacing_gl]:
+    for width in [spacing_bed, spacing_speed, spacing_edge,
+                  spacing_gl, spacing_coast]:
         cell_width = np.minimum(cell_width, width)
 
     # Set large cell_width in areas we are going to cull anyway (speeds up
@@ -420,8 +447,7 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y, section_name,
     """
     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``.
+    :py:func:`compass.landice.mesh.set_cell_width()`.
 
     Parameters
     ----------
@@ -465,6 +491,11 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y, section_name,
 
     dist_to_grounding_line : numpy.ndarray
         Distance from each cell to the grounding line
+
+    dist_to_coast : numpy.ndarray
+        Distance from each cell to the coast, defined as the last cell
+        adjacent to ocean (ice free with bed < 0) whether ice-filled
+        or ice-free.
     """
 
     logger = self.logger
@@ -495,12 +526,14 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y, section_name,
     ])
 
     ice_mask = thk > 0.0
+    ocean_mask = np.logical_and(thk == 0., topg < 0.)
     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(thk.shape, dtype=bool)
     grounding_line_mask = np.zeros(thk.shape, dtype=bool)
+    coast_mask = np.zeros(thk.shape, dtype=bool)
 
     for n in neighbors:
         shifted_ice = np.roll(ice_mask, n, axis=[0, 1])
@@ -512,17 +545,24 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y, section_name,
         grounding_line_mask = np.logical_or(grounding_line_mask,
                                             not_grounded_mask)
 
+        shifted_ocean = np.roll(ocean_mask, n, axis=[0, 1])
+        coast_mask = np.logical_or(coast_mask, shifted_ocean)
+
     # where ice exists and neighbors non-ice locations
     margin_mask = np.logical_and(margin_mask, ice_mask)
     # where grounded ice exists and neighbors floating ice
     grounding_line_mask = np.logical_and(grounding_line_mask, grounded_mask)
+    # where non-ocean exists and neighbors ocean locations
+    coast_mask = np.logical_and(coast_mask, np.logical_not(ocean_mask))
 
-    fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(6, 3))
+    fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(9, 3))
     margin_plot = ax[0].pcolor(margin_mask)
     gl_plot = ax[1].pcolor(grounding_line_mask)  # noqa F841
+    coast_plot = ax[2].pcolor(coast_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[2].set_title("coast mask")
+    plt.colorbar(margin_plot, ax=ax, shrink=0.7)
     [ax.set_aspect('equal') for ax in ax]
     fig.savefig("masks.png", dpi=400)
 
@@ -534,19 +574,23 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y, section_name,
     dist_to_grounding_line = distance_transform_edt(
         ~grounding_line_mask, sampling=(dy, dx)
     )
+    dist_to_coast = distance_transform_edt(
+        ~coast_mask, sampling=(dy, dx)
+    )
 
     # Limit maximum distance
     max_dist = float(np.ceil(window_size / max(dx, dy))) * max(dx, dy)
     dist_to_edge = np.minimum(dist_to_edge, max_dist)
     dist_to_grounding_line = np.minimum(dist_to_grounding_line, max_dist)
+    dist_to_coast = np.minimum(dist_to_coast, max_dist)
 
     toc = time.time()
     logger.info(
         'compass.landice.mesh.get_dist_to_edge_and_gl() took '
         f'{toc - tic:0.2f} seconds'
     )
 
-    return dist_to_edge, dist_to_grounding_line
+    return dist_to_edge, dist_to_grounding_line, dist_to_coast
 
 
 def build_cell_width(self, section_name, gridded_dataset,
@@ -628,7 +672,7 @@ def build_cell_width(self, section_name, gridded_dataset,
 
     # 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(
+    distToEdge, distToGL, distToCoast = get_dist_to_edge_and_gl(
         self, thk, topg, x1,
         y1, section_name=section_name)
 
@@ -637,6 +681,7 @@ def build_cell_width(self, section_name, gridded_dataset,
                                 thk=thk, bed=topg, vx=vx, vy=vy,
                                 dist_to_edge=distToEdge,
                                 dist_to_grounding_line=distToGL,
+                                dist_to_coast=distToCoast,
                                 flood_fill_iStart=flood_fill_start[0],
                                 flood_fill_jStart=flood_fill_start[1])
 
@@ -683,9 +728,10 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points,
         following options to be set in the given config section:
         ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``,
         ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``,
-        ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``,
-        ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``,
-        ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``.
+        ``high_dist``, ``low_dist``, ``high_dist_coast``, ``low_dist_coast``,
+        ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``,
+        ``cull_distance``, ``use_speed``, ``use_dist_to_edge``,
+        ``use_dist_to_grounding_line``, ``use_dist_to_coast``, and ``use_bed``.
         See the Land-Ice Framework section of the Users or Developers guide
         for more information about these options and their uses.
 

compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg

@@ -15,13 +15,13 @@ y_max = None
 # distance from ice margin to cull (km).
 # Set to a value <= 0 if you do not want
 # to cull based on distance from margin.
-cull_distance = 10.0
+cull_distance = -100.0
 
 # mesh density parameters
 # minimum cell spacing (meters)
-min_spac = 3.e3
+min_spac = 4.e3
 # maximum cell spacing (meters)
-max_spac = 30.e3
+max_spac = 40.e3
 # log10 of max speed for cell spacing
 high_log_speed = 2.5
 # log10 of min speed for cell spacing
@@ -30,6 +30,10 @@ low_log_speed = 0.75
 high_dist = 1.e5
 # distance within which cell spacing = min_spac
 low_dist = 1.e4
+# distance at which cell spacing in ocean = max_spac
+high_dist_coast = 16.e3
+# distance within which cell spacing in ocean = min_spac
+low_dist_coast = 8.e3
 # distance at which bed topography has no effect
 high_dist_bed = 1.e5
 # distance within which bed topography has maximum effect
@@ -43,6 +47,7 @@ high_bed = 100.0
 use_speed = True
 use_dist_to_grounding_line = False
 use_dist_to_edge = True
+use_dist_to_coast = True
 use_bed = True
 
 [greenland]

compass/landice/tests/greenland/mesh_gen/mesh_gen_1to10km.cfg

@@ -0,0 +1,73 @@
+# config options for high_res_mesh test case
+[mesh]
+
+# number of levels in the mesh
+levels = 10
+
+# Bounds of GIS mesh. If any bound is set
+# to None, the mesh will use the full extent
+# of the gridded dataset.
+x_min = None
+x_max = None
+y_min = None
+y_max = None
+
+# distance from ice margin to cull (km).
+# Set to a value <= 0 if you do not want
+# to cull based on distance from margin.
+cull_distance = -100.0
+
+# mesh density parameters
+# minimum cell spacing (meters)
+min_spac = 1.e3
+# maximum cell spacing (meters)
+max_spac = 10.e3
+# log10 of max speed for cell spacing
+high_log_speed = 2.5
+# log10 of min speed for cell spacing
+low_log_speed = 0.75
+# distance at which cell spacing = max_spac
+high_dist = 1.e5
+# distance within which cell spacing = min_spac
+low_dist = 1.e4
+# distance at which cell spacing in ocean = max_spac
+high_dist_coast = 10.e3
+# distance within which cell spacing in ocean = min_spac
+low_dist_coast = 2.e3
+# distance at which bed topography has no effect
+high_dist_bed = 1.e5
+# distance within which bed topography has maximum effect
+low_dist_bed = 7.5e4
+# Bed elev beneath which cell spacing is minimized
+low_bed = 50.0
+# Bed elev above which cell spacing is maximized
+high_bed = 100.0
+
+# mesh density functions
+use_speed = True
+use_dist_to_grounding_line = False
+use_dist_to_edge = True
+use_dist_to_coast = True
+use_bed = True
+
+[greenland]
+# Whether to interpolate data (controls run_optional_interpolation)
+interpolate_data = True
+# path to directory containing BedMachine and Measures datasets
+# (default value is for Perlmutter)
+data_path = /global/cfs/cdirs/fanssie/standard_datasets/GIS_datasets/
+
+# filename of the BedMachine thickness and bedTopography dataset
+# (default value is for Perlmutter)
+bedmachine_filename = BedMachineGreenland-v6_edits_floodFill_extrap.nc
+
+# filename of the MEaSUREs ice velocity dataset
+# (default value is for Perlmutter)
+measures_filename = greenland_vel_mosaic500_extrap.nc
+
+# projection of the source datasets, according to the dictionary keys
+# create_scrip_file_from_planar_rectangular_grid from MPAS_Tools
+src_proj = gis-gimp
+
+# number of processors to use for ESMF_RegridWeightGen
+nProcs = 2048