Interpolate BedMachine v5 and MEaSUREs data onto Humboldt mesh

compass/landice/tests/humboldt/mesh.py

@@ -1,4 +1,15 @@
-from compass.landice.mesh import build_cell_width, build_mali_mesh
+import os
+
+import numpy as np
+import xarray as xr
+from mpas_tools.scrip.from_mpas import scrip_from_mpas
+
+from compass.landice.mesh import (
+    build_cell_width,
+    build_mali_mesh,
+    clean_up_after_interp,
+    interp_gridded2mali,
+)
 from compass.model import make_graph_file
 from compass.step import Step
 
@@ -27,8 +38,9 @@ def __init__(self, test_case):
         super().__init__(test_case=test_case, name='mesh', cpus_per_task=128,
                          min_cpus_per_task=1)
 
+        self.mesh_filename = 'Humboldt.nc'
         self.add_output_file(filename='graph.info')
-        self.add_output_file(filename='Humboldt.nc')
+        self.add_output_file(filename=self.mesh_filename)
         self.add_input_file(
             filename='humboldt_1km_2024_01_29.epsg3413.icesheetonly.nc',
             target='humboldt_1km_2024_01_29.epsg3413.icesheetonly.nc',
@@ -49,21 +61,71 @@ def run(self):
         """
         logger = self.logger
         section_name = 'mesh'
-        mesh_name = 'Humboldt.nc'
+        section_gis = self.config['greenland']
+
+        nProcs = section_gis.get('nProcs')
+        src_proj = section_gis.get("src_proj")
+        data_path = section_gis.get('data_path')
+        measures_filename = section_gis.get("measures_filename")
+        bedmachine_filename = section_gis.get("bedmachine_filename")
+
+        measures_dataset = os.path.join(data_path, measures_filename)
+        bedmachine_dataset = os.path.join(data_path, bedmachine_filename)
 
         logger.info('calling build_cell_width')
         cell_width, x1, y1, geom_points, geom_edges, floodMask = \
             build_cell_width(
                 self, section_name=section_name,
-                gridded_dataset='greenland_2km_2024_01_29.epsg3413.nc')
+                gridded_dataset='greenland_2km_2024_01_29.epsg3413.nc',
+                flood_fill_start=[100, 700])
 
         build_mali_mesh(
             self, cell_width, x1, y1, geom_points, geom_edges,
-            mesh_name=mesh_name, section_name=section_name,
+            mesh_name=self.mesh_filename, section_name=section_name,
             gridded_dataset='humboldt_1km_2024_01_29.epsg3413.icesheetonly.nc',
             projection='gis-gimp', geojson_file='Humboldt.geojson',
             cores=self.cpus_per_task)
 
+        # Create scrip file for the newly generated mesh
+        logger.info('creating scrip file for destination mesh')
+        dst_scrip_file = f"{self.mesh_filename.split('.')[:-1][0]}_scrip.nc"
+        scrip_from_mpas(self.mesh_filename, dst_scrip_file)
+
+        # Now perform bespoke interpolation of geometry and velocity data
+        # from their respective sources
+        interp_gridded2mali(self, bedmachine_dataset, dst_scrip_file, nProcs,
+                            self.mesh_filename, src_proj, variables="all")
+
+        # only interpolate a subset of MEaSUREs variables onto the MALI mesh
+        measures_vars = ['observedSurfaceVelocityX',
+                         'observedSurfaceVelocityY',
+                         'observedSurfaceVelocityUncertainty']
+        interp_gridded2mali(self, measures_dataset, dst_scrip_file, nProcs,
+                            self.mesh_filename, src_proj,
+                            variables=measures_vars)
+
+        # perform some final cleanup details
+        clean_up_after_interp(self.mesh_filename)
         logger.info('creating graph.info')
-        make_graph_file(mesh_filename=mesh_name,
+        make_graph_file(mesh_filename=self.mesh_filename,
                         graph_filename='graph.info')
+
+        # Do some final validation of the mesh
+        ds = xr.open_dataset(self.mesh_filename)
+        # Ensure basalHeatFlux is positive
+        ds["basalHeatFlux"] = np.abs(ds.basalHeatFlux)
+        # Ensure reasonable dHdt values
+        dHdt = ds["observedThicknessTendency"]
+        # Arbitrary 5% uncertainty; improve this later
+        dHdtErr = np.abs(dHdt) * 0.05
+        # Use threshold of |dHdt| > 1.0 to determine invalid data
+        mask = np.abs(dHdt) > 1.0
+        # Assign very large uncertainty where data is missing
+        dHdtErr = dHdtErr.where(~mask, 1.0)
+        # Remove ridiculous values
+        dHdt = dHdt.where(~mask, 0.0)
+        # Put the updated fields back in the dataset
+        ds["observedThicknessTendency"] = dHdt
+        ds["observedThicknessTendencyUncertainty"] = dHdtErr
+        # Write the data to disk
+        ds.to_netcdf(self.mesh_filename, 'a')

compass/landice/tests/humboldt/mesh_gen/mesh_gen.cfg

@@ -42,3 +42,23 @@ use_speed = True
 use_dist_to_grounding_line = False
 use_dist_to_edge = True
 use_bed = True
+
+[greenland]
+# 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-v5_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.py from MPAS_Tools
+src_proj = gis-gimp
+
+# number of processors to use for ESMF_RegridWeightGen
+nProcs = 128