Unify gridded dataset preprocessing to work for boht GIS/AIS.

compass/landice/mesh.py

@@ -889,16 +889,15 @@ def add_bedmachine_thk_to_ais_gridded_data(self, source_gridded_dataset,
     return gridded_dataset_with_bm_thk
 
 
-def preprocess_ais_data(self, source_gridded_dataset,
-                        floodFillMask):
+def preprocess_gridded_data(self, source_gridded_dataset, floodFillMask):
     """
-    Perform adjustments to gridded AIS datasets needed
+    Perform adjustments to gridded AIS/GIS datasets needed
     for rest of compass workflow to utilize them
 
     Parameters
     ----------
     source_gridded_dataset : str
-        name of NetCDF file containing original AIS gridded datasets
+        name of NetCDF file containing original AIS/GIS gridded datasets
 
     floodFillMask : numpy.ndarray
         0/1 mask of flood filled ice region
@@ -934,10 +933,16 @@ def preprocess_ais_data(self, source_gridded_dataset,
     data.createVariable('iceMask', 'f', ('time', 'y1', 'x1'))
     data.variables['iceMask'][:] = data.variables["thk"][:] > 0.
 
+    # rename variables to common/shared names
+    if ('dHdt' not in data.variables) & ('dhdt' in data.variables):
+        data.renameVariable('dhdt', 'dHdt')
+    if ('topgerr' not in data.variables) & ('thkerr' in data.variables):
+        data.renameVariable('thkerr', 'topgerr')
+
     # Note: dhdt is only reported over grounded ice, so we will have to
     # either update the dataset to include ice shelves or give them values of
     # 0 with reasonably large uncertainties.
-    dHdt = data.variables["dhdt"][:]
+    dHdt = data.variables["dHdt"][:]
     dHdtErr = 0.05 * dHdt  # assign arbitrary uncertainty of 5%
     # Where dHdt data are missing, set large uncertainty
     dHdtErr[dHdt > 1.e30] = 1.
@@ -952,16 +957,16 @@ def preprocess_ais_data(self, source_gridded_dataset,
     yy = yGrid.ravel()
     bigTic = time.perf_counter()
     for field in ['thk', 'bheatflx', 'vx', 'vy',
-                  'ex', 'ey', 'thkerr', 'dhdt']:
+                  'ex', 'ey', 'topgerr', 'dHdt']:
         tic = time.perf_counter()
         logger.info(f"Beginning building interpolator for {field}")
-        if field in ['thk', 'thkerr']:
+        if field in ['thk', 'topgerr']:
             mask = cellsWithIce.ravel()
         elif field == 'bheatflx':
             mask = np.logical_and(
                 data.variables[field][:].ravel() < 1.0e9,
                 data.variables[field][:].ravel() != 0.0)
-        elif field in ['vx', 'vy', 'ex', 'ey', 'dhdt']:
+        elif field in ['vx', 'vy', 'ex', 'ey', 'dHdt']:
             mask = np.logical_and(
                 data.variables[field][:].ravel() < 1.0e9,
                 cellsWithIce.ravel() > 0)
@@ -996,9 +1001,6 @@ def preprocess_ais_data(self, source_gridded_dataset,
     data.variables['subm'][:] *= -1.0  # correct basal melting sign
     data.variables['subm_ss'][:] *= -1.0
 
-    data.renameVariable('dhdt', 'dHdt')
-    data.renameVariable('thkerr', 'topgerr')
-
     data.createVariable('x', 'f', ('x1'))
     data.createVariable('y', 'f', ('y1'))
     data.variables['x'][:] = x1

compass/landice/tests/antarctica/mesh.py

@@ -11,7 +11,7 @@
     clean_up_after_interp,
     interp_gridded2mali,
     make_region_masks,
-    preprocess_ais_data,
+    preprocess_gridded_data,
 )
 from compass.model import make_graph_file
 from compass.step import Step
@@ -88,8 +88,8 @@ def run(self):
 
         # Preprocess the gridded AIS source datasets to work
         # with the rest of the workflow
-        logger.info('calling preprocess_ais_data')
-        preprocessed_gridded_dataset = preprocess_ais_data(
+        logger.info('calling preprocess_gridded_data')
+        preprocessed_gridded_dataset = preprocess_gridded_data(
             self, bm_updated_gridded_dataset, floodFillMask)
 
         # Now build the base mesh and perform the standard interpolation

compass/landice/tests/greenland/mesh.py

@@ -1,7 +1,5 @@
 import os
 
-import netCDF4
-import numpy as np
 from mpas_tools.scrip.from_mpas import scrip_from_mpas
 
 from compass.landice.mesh import (
@@ -10,6 +8,7 @@
     clean_up_after_interp,
     interp_gridded2mali,
     make_region_masks,
+    preprocess_griddedd_data,
 )
 from compass.model import make_graph_file
 from compass.step import Step
@@ -84,11 +83,17 @@ def run(self):
                 gridded_dataset=source_gridded_dataset_2km,
                 flood_fill_start=[100, 700])
 
+        # Preprocess the gridded GrIS source datasets to work
+        # with the rest of the workflow
+        logger.info('calling preprocess_gridded_data')
+        preprocessed_gridded_dataset = preprocess_griddedd_data(
+            self, source_gridded_dataset_1km, floodMask)
+
         # Now build the base mesh and perform the standard interpolation
         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,
+            gridded_dataset=preprocessed_gridded_dataset,
             projection=src_proj, geojson_file=None)
 
         # Create scrip file for the newly generated mesh
@@ -129,22 +134,3 @@ def run(self):
                                 'southGreenland',
                                 'southWestGreenland',
                                 'westCentralGreenland'])
-
-        # is there a way to encompass this in an existing framework function?
-        data = netCDF4.Dataset('GIS.nc', 'r+')
-        data.set_auto_mask(False)
-        # Ensure basalHeatFlux is positive
-        data.variables['basalHeatFlux'][:] = np.abs(
-            data.variables['basalHeatFlux'][:])
-        # Ensure reasonable dHdt values
-        dHdt = data.variables["observedThicknessTendency"][:]
-        dHdtErr = data.variables["observedThicknessTendencyUncertainty"][:]
-        # Arbitrary 5% uncertainty; improve this later
-        dHdtErr = np.abs(dHdt) * 0.05
-        # large uncertainty where data is missing
-        dHdtErr[np.abs(dHdt) > 1.0] = 1.0
-        dHdt[np.abs(dHdt) > 1.0] = 0.0  # Remove ridiculous values
-        data.variables["observedThicknessTendency"][:] = dHdt
-        data.variables["observedThicknessTendencyUncertainty"][:] = dHdtErr
-
-        data.close()