From f50b00d28189daa99797574e171071d0bb32eb37 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 21 May 2024 12:36:50 -0700 Subject: [PATCH] Unify gridded dataset preprocessing to work for boht GIS/AIS. --- compass/landice/mesh.py | 24 ++++++++++--------- compass/landice/tests/antarctica/mesh.py | 6 ++--- compass/landice/tests/greenland/mesh.py | 30 +++++++----------------- 3 files changed, 24 insertions(+), 36 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 56736113bc..769697c29e 100644 --- a/compass/landice/mesh.py +++ b/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 diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index 4b0042b92e..d5c825bea8 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/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 diff --git a/compass/landice/tests/greenland/mesh.py b/compass/landice/tests/greenland/mesh.py index aa89600242..86583c7aed 100644 --- a/compass/landice/tests/greenland/mesh.py +++ b/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()