Merge pull request #876 from andrewdnolan/moderate_refinement

Enable MALI Greenland mesh overlap with ocean mesh, through a buffer defined by a geojson file

Author: trhille

compass/landice/mesh.py

@@ -189,6 +189,61 @@ def set_rectangular_geom_points_and_edges(xmin, xmax, ymin, ymax):
     return geom_points, geom_edges
 
 
+def clip_mesh_to_bounding_box(mask_ds, base_ds, bounding_box):
+    """
+    Set cells to culled if they lay outside the bounding box
+
+    Parameters
+    ----------
+    mask_ds: xr.Dataset
+        mask dataset, generated by ``compute_mpas_region_mask``
+    base_ds: xr.Dataset
+        unculled mesh dataset
+    bounding_box: list of 4 ints
+        Bounding box [x_min, x_max, y_min, y_max] to cull mesh outside of
+
+    Returns
+    -------
+    mask_ds: xarray.Dataset
+        mask dataset with updated masks based on bounding box
+    """
+
+    if len(bounding_box) != 4:
+        msg = f"bounding box must be len 4, instead is len {len(bounding_box)}"
+        raise ValueError(msg)
+
+    x_min, x_max, y_min, y_max = bounding_box
+
+    if (x_max < x_min) or (y_max < y_min):
+        msg = "Bounding box must be ordered: [x_min, x_max, y_min, y_max]"
+        msg += (
+            f"\n x_max < x_min ({x_max:.2f} < {x_min:.2f})"
+            if x_max < x_min else ""
+        )
+        msg += (
+            f"\n y_max < y_min ({y_max:.2f} < {y_min:.2f})"
+            if y_max < y_min else ""
+        )
+        raise ValueError(msg)
+
+    for loc in ["Cell", "Edge", "Vertex"]:
+        mask_var = f"region{loc}Masks"
+
+        if mask_var not in mask_ds:
+            continue
+
+        mask = (
+            (base_ds[f"x{loc}"] < x_min) |
+            (base_ds[f"x{loc}"] > x_max) |
+            (base_ds[f"y{loc}"] < y_min) |
+            (base_ds[f"y{loc}"] > y_max)
+        )
+
+        mask_ds[mask_var] = xarray.where(mask, 0, mask_ds[mask_var])
+
+    return mask_ds
+
+
 def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
                    dist_to_edge=None, dist_to_grounding_line=None,
                    flood_fill_iStart=None, flood_fill_jStart=None):
@@ -618,7 +673,7 @@ def build_cell_width(self, section_name, gridded_dataset,
 def build_mali_mesh(self, cell_width, x1, y1, geom_points,
                     geom_edges, mesh_name, section_name,
                     gridded_dataset, projection, geojson_file=None,
-                    cores=1):
+                    cores=1, bounding_box=None):
     """
     Create the MALI mesh based on final cell widths determined by
     :py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and
@@ -672,8 +727,18 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points,
 
     cores : int, optional
         The number of cores to use for mask creation
+
+    bounding_box : array_like of float, shape (4,), optional
+        Bounding box [x_min, x_max, y_min, y_max] to cull mesh outside of
     """
 
+    if bounding_box is not None and geojson_file is None:
+        msg = (
+            "Bounding box clipping can only be applied to an existing cull"
+            "mask. You must provide a geojson file for this to work."
+        )
+        raise ValueError(msg)
+
     logger = self.logger
     section = self.config[section_name]
 
@@ -735,6 +800,10 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points,
     dsMesh = xarray.open_dataset('grid_preCull.nc')
     if geojson_file is not None:
         mask = xarray.open_dataset('mask.nc')
+
+        if bounding_box is not None:
+            mask = clip_mesh_to_bounding_box(mask, dsMesh, bounding_box)
+
     else:
         mask = None
 

compass/landice/tests/greenland/greenland_only_outline_45km_buffer_latlon_singlepart.geojson

@@ -54,7 +54,17 @@ def __init__(self, test_case):
                             target='greenland_2km_2024_01_29.epsg3413.nc',
                             database='')
 
-    # no setup() method is needed
+    def setup(self):
+        """
+        Add the configured geojson cull-mask file as an input.
+        """
+        section = self.config['greenland']
+        geojson_filename = section.get('geojson_filename')
+
+        self.add_input_file(filename=geojson_filename,
+                            package='compass.landice.tests.greenland',
+                            target=geojson_filename,
+                            database=None)
 
     def run(self):
         """
@@ -71,10 +81,13 @@ def run(self):
         data_path = section_gis.get('data_path')
         measures_filename = section_gis.get("measures_filename")
         bedmachine_filename = section_gis.get("bedmachine_filename")
+        geojson_filename = section_gis.get('geojson_filename')
 
         measures_dataset = os.path.join(data_path, measures_filename)
         bedmachine_dataset = os.path.join(data_path, bedmachine_filename)
 
+        bounding_box = self._get_bedmachine_bounding_box(bedmachine_dataset)
+
         section_name = 'mesh'
 
         source_gridded_dataset_1km = 'greenland_1km_2024_01_29.epsg3413.icesheetonly.nc'  # noqa: E501
@@ -91,8 +104,10 @@ def run(self):
         build_mali_mesh(
             self, cell_width, x1, y1, geom_points, geom_edges,
             mesh_name=self.mesh_filename, section_name=section_name,
-            gridded_dataset=source_gridded_dataset_1km,
-            projection=src_proj, geojson_file=None)
+            gridded_dataset=source_gridded_dataset_1km, projection=src_proj,
+            geojson_file=geojson_filename,
+            bounding_box=bounding_box,
+        )
 
         # Create scrip file for the newly generated mesh
         logger.info('creating scrip file for destination mesh')
@@ -161,3 +176,14 @@ def run(self):
         ds["observedThicknessTendencyUncertainty"] = dHdtErr
         # Write the data to disk
         ds.to_netcdf(self.mesh_filename, 'a')
+
+    def _get_bedmachine_bounding_box(self, bedmachine_filepath):
+
+        ds = xr.open_dataset(bedmachine_filepath)
+
+        x_min = ds.x1.min()
+        x_max = ds.x1.max()
+        y_min = ds.y1.min()
+        y_max = ds.y1.max()
+
+        return [x_min, x_max, y_min, y_max]

compass/landice/tests/greenland/mesh.py

@@ -53,6 +53,9 @@ use_bed = True
 # (default value is for Perlmutter)
 data_path = /global/cfs/cdirs/fanssie/standard_datasets/GIS_datasets/
 
+# geojson used to create the cull mask in mesh generation
+geojson_filename = greenland_only_outline_45km_buffer_latlon_singlepart.geojson
+
 # filename of the BedMachine thickness and bedTopography dataset
 # (default value is for Perlmutter)
 bedmachine_filename = BedMachineGreenland-v6_edits_floodFill_extrap.nc

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

@@ -222,6 +222,7 @@ greenland
    run_model.RunModel.run
 
    mesh.Mesh
+   mesh.Mesh.setup
    mesh.Mesh.run
 
    mesh_gen.MeshGen
@@ -491,6 +492,7 @@ Landice Framework
 
    mesh.add_bedmachine_thk_to_ais_gridded_data
    mesh.clean_up_after_interp
+   mesh.clip_mesh_to_bounding_box
    mesh.gridded_flood_fill
    mesh.interp_gridded2mali
    mesh.mpas_flood_fill

docs/developers_guide/landice/api.rst

@@ -53,6 +53,13 @@ thickness to the final MALI mesh.
 :py:func:`compass.landice.mesh.clean_up_after_interp()` performs some final
 clean up steps after interpolation for the AIS mesh case.
 
+:py:func:`compass.landice.mesh.clip_mesh_to_bounding_box()` updates an
+existing MPAS region mask so cells, edges, and vertices outside a bounding
+box ``[x_min, x_max, y_min, y_max]`` are culled.  This is useful when a
+``geojson_file`` supplies the initial cull mask but the final culled mesh
+should also be limited to the valid extent of an input gridded dataset.
+The bounding box must be ordered as ``[x_min, x_max, y_min, y_max]``.
+
 :py:func:`compass.landice.mesh.gridded_flood_fill()` applies a flood-fill algorithm
 to the gridded dataset in order to separate the ice sheet from peripheral ice.
 
@@ -87,7 +94,12 @@ based on user-defined density functions and config options.
 cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw
 and MPAS-Tools functions. Culls the mesh based on config options, interpolates
 all available fields from the gridded dataset to the MALI mesh using the bilinear
-method, and marks domain boundaries as Dirichlet cells.
+method, and marks domain boundaries as Dirichlet cells.  When a ``geojson_file``
+is supplied, :py:func:`compass.landice.mesh.build_mali_mesh()` can also take an
+optional ``bounding_box`` argument to clip the cull mask to
+``[x_min, x_max, y_min, y_max]`` before culling.  This additional clipping is
+applied to the mask, not to the rectangular geometry passed to
+:py:func:`compass.landice.mesh.set_rectangular_geom_points_and_edges()`.
 
 :py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions
 defined in Geometric Features repository. It is only used by the ``antarctica``