diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 33fcc9931f..0427e47556 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -1,4 +1,5 @@ import os +import re import sys import time from shutil import copyfile @@ -12,6 +13,7 @@ from mpas_tools.logging import check_call from mpas_tools.mesh.conversion import convert, cull from mpas_tools.mesh.creation import build_planar_mesh +from mpas_tools.mesh.creation.sort_mesh import sort_mesh from netCDF4 import Dataset from scipy.interpolate import NearestNDInterpolator, interpn @@ -442,7 +444,7 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y, logger.info('WARNING: window_size was set to a value smaller' ' than high_dist and/or high_dist_bed. Resetting' f' window_size to {max(high_dist, high_dist_bed)},' - ' which is max(high_dist, high_dist_bed)') + ' which is max(high_dist, high_dist_bed)') window_size = max(high_dist, high_dist_bed) dx = x[1] - x[0] # assumed constant and equal in x and y @@ -744,10 +746,11 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, check_call(args, logger=logger) - logger.info('culling and converting') + logger.info('culling, converting, and sorting') dsMesh = xarray.open_dataset('culled.nc') dsMesh = cull(dsMesh, logger=logger) dsMesh = convert(dsMesh, logger=logger) + dsMesh = sort_mesh(dsMesh) write_netcdf(dsMesh, 'dehorned.nc') args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', @@ -968,6 +971,8 @@ def preprocess_ais_data(self, source_gridded_dataset, tic = time.perf_counter() logger.info(f"Beginning interpolation for {field}") + # NOTE: Do not need to evaluate the extrapolator at all grid cells. + # Only needed for ice-free grid cells, since is NN extrapolation data.variables[field][0, :] = interp(xGrid, yGrid) toc = time.perf_counter() logger.info(f"Interpolation completed in {toc - tic} seconds") @@ -1002,15 +1007,15 @@ def preprocess_ais_data(self, source_gridded_dataset, return preprocessed_gridded_dataset -def interp_ais_bedmachine(self, data_path, mali_scrip, nProcs, dest_file): +def interp_gridded2mali(self, source_file, mali_scrip, nProcs, dest_file, proj, + variables="all"): """ - Interpolates BedMachine thickness and bedTopography dataset - to a MALI mesh + Interpolate gridded dataset (e.g. MEASURES, BedMachine) onto a MALI mesh Parameters ---------- - data_path : str - path to AIS datasets, including BedMachine + source_file : str + filepath to the source gridded datatset to be interpolated mali_scrip : str name of scrip file corresponding to destination MALI mesh @@ -1020,104 +1025,69 @@ def interp_ais_bedmachine(self, data_path, mali_scrip, nProcs, dest_file): dest_file: str MALI input file to which data should be remapped - """ - - logger = self.logger - - logger.info('creating scrip file for BedMachine dataset') - # Note: writing scrip file to workdir - args = ['create_SCRIP_file_from_planar_rectangular_grid.py', - '-i', - os.path.join(data_path, - 'BedMachineAntarctica_2020-07-15_v02_edits_floodFill_extrap_fillVostok.nc'), # noqa - '-s', - 'BedMachineAntarctica_2020-07-15_v02.scrip.nc', - '-p', 'ais-bedmap2', - '-r', '2'] - check_call(args, logger=logger) - - # Generate remapping weights - # Testing shows 5 badger/grizzly nodes works well. - # 2 nodes is too few. I have not tested anything in between. - logger.info('generating gridded dataset -> MPAS weights') - args = ['srun', '-n', nProcs, 'ESMF_RegridWeightGen', - '--source', - 'BedMachineAntarctica_2020-07-15_v02.scrip.nc', - '--destination', mali_scrip, - '--weight', 'BedMachine_to_MPAS_weights.nc', - '--method', 'conserve', - "--netcdf4", - "--dst_regional", "--src_regional", '--ignore_unmapped'] - check_call(args, logger=logger) - - # Perform actual interpolation using the weights - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - os.path.join(data_path, - 'BedMachineAntarctica_2020-07-15_v02_edits_floodFill_extrap_fillVostok.nc'), # noqa - '-d', dest_file, - '-m', 'e', - '-w', 'BedMachine_to_MPAS_weights.nc'] - check_call(args, logger=logger) + proj: str + projection of the source dataset -def interp_ais_measures(self, data_path, mali_scrip, nProcs, dest_file): + variables: "all" or list of strings + either the string "all" or a list of strings """ - Interpolates MEASURES ice velocity dataset - to a MALI mesh - Parameters - ---------- - data_path : str - path to AIS datasets, including BedMachine + def __guess_scrip_name(filename): - mali_scrip : str - name of scrip file corresponding to destination MALI mesh + # try searching for string followed by a version number + match = re.search(r'(^.*[_-]v\d*[_-])+', filename) - nProcs : int - number of processors to use for generating remapping weights + if match: + # slice string to end of match minus one to leave of final _ or - + base_fn = filename[:match.end() - 1] + else: + # no matches were found, just use the filename (minus extension) + base_fn = os.path.splitext(filename)[0] - dest_file: str - MALI input file to which data should be remapped - """ + return f"{base_fn}.scrip.nc" logger = self.logger - logger.info('creating scrip file for velocity dataset') + source_scrip = __guess_scrip_name(os.path.basename(source_file)) + weights_filename = "gridded_to_MPAS_weights.nc" + + # make sure variables is a list, encompasses the variables="all" case + if isinstance(variables, str): + variables = [variables] + if not isinstance(variables, list): + raise TypeError("Arugment 'variables' is of incorrect type, must" + " either the string 'all' or a list of strings") + + logger.info('creating scrip file for source dataset') # Note: writing scrip file to workdir args = ['create_SCRIP_file_from_planar_rectangular_grid.py', - '-i', - os.path.join(data_path, - 'antarctica_ice_velocity_450m_v2_edits_extrap.nc'), - '-s', - 'antarctica_ice_velocity_450m_v2.scrip.nc', - '-p', 'ais-bedmap2', + '-i', source_file, + '-s', source_scrip, + '-p', proj, '-r', '2'] check_call(args, logger=logger) # Generate remapping weights logger.info('generating gridded dataset -> MPAS weights') args = ['srun', '-n', nProcs, 'ESMF_RegridWeightGen', - '--source', - 'antarctica_ice_velocity_450m_v2.scrip.nc', + '--source', source_scrip, '--destination', mali_scrip, - '--weight', 'measures_to_MPAS_weights.nc', + '--weight', weights_filename, '--method', 'conserve', "--netcdf4", "--dst_regional", "--src_regional", '--ignore_unmapped'] check_call(args, logger=logger) + # Perform actual interpolation using the weights logger.info('calling interpolate_to_mpasli_grid.py') args = ['interpolate_to_mpasli_grid.py', - '-s', - os.path.join(data_path, - 'antarctica_ice_velocity_450m_v2_edits_extrap.nc'), + '-s', source_file, '-d', dest_file, '-m', 'e', - '-w', 'measures_to_MPAS_weights.nc', - '-v', 'observedSurfaceVelocityX', - 'observedSurfaceVelocityY', - 'observedSurfaceVelocityUncertainty'] + '-w', weights_filename, + '-v'] + variables + check_call(args, logger=logger) diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index 609bb7275c..fcff251304 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/compass/landice/tests/antarctica/mesh.py @@ -9,8 +9,7 @@ build_cell_width, build_mali_mesh, clean_up_after_interp, - interp_ais_bedmachine, - interp_ais_measures, + interp_gridded2mali, make_region_masks, preprocess_ais_data, ) @@ -44,10 +43,10 @@ def __init__(self, test_case): self.mesh_filename = 'Antarctica.nc' self.add_output_file(filename='graph.info') self.add_output_file(filename=self.mesh_filename) - self.add_output_file(filename=f'{self.mesh_filename[:-3]}_' - f'imbie_regionMasks.nc') - self.add_output_file(filename=f'{self.mesh_filename[:-3]}_' - f'ismip6_regionMasks.nc') + self.add_output_file( + filename=f'{self.mesh_filename[:-3]}_imbie_regionMasks.nc') + self.add_output_file( + filename=f'{self.mesh_filename[:-3]}_ismip6_regionMasks.nc') self.add_input_file( filename='antarctica_8km_2024_01_29.nc', target='antarctica_8km_2024_01_29.nc', @@ -61,37 +60,38 @@ def run(self): """ logger = self.logger config = self.config + section_ais = config['antarctica'] - data_path = section_ais.get('data_path') + nProcs = section_ais.get('nProcs') + src_proj = section_ais.get("src_proj") + data_path = section_ais.get('data_path') + measures_filename = section_ais.get("measures_filename") + bedmachine_filename = section_ais.get("bedmachine_filename") + + measures_dataset = os.path.join(data_path, measures_filename) + bedmachine_dataset = os.path.join(data_path, bedmachine_filename) section_name = 'mesh' + # TODO: do we want to add this to the config file? source_gridded_dataset = 'antarctica_8km_2024_01_29.nc' - bedmachine_path = os.path.join( - data_path, - 'BedMachineAntarctica_2020-07-15_v02_edits_floodFill_extrap_fillVostok.nc') # noqa bm_updated_gridded_dataset = add_bedmachine_thk_to_ais_gridded_data( - self, source_gridded_dataset, bedmachine_path) + self, source_gridded_dataset, bedmachine_dataset) + logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodFillMask = \ build_cell_width( self, section_name=section_name, gridded_dataset=bm_updated_gridded_dataset) - # 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( - self, bm_updated_gridded_dataset, floodFillMask) - # 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=bm_updated_gridded_dataset, - projection='ais-bedmap2', geojson_file=None) + projection=src_proj, geojson_file=None) # Now that we have base mesh with standard interpolation # perform advanced interpolation for specific fields @@ -107,12 +107,19 @@ def run(self): data.variables['iceMask'][:] = 0. data.close() - # interpolate fields from composite dataset - # Note: this was already done in build_mali_mesh() using - # bilinear interpolation. Redoing it here again is likely - # not needed. Also, it should be assessed if bilinear or - # barycentric used here is preferred for this application. - # Current thinking is they are both equally appropriate. + # 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( + self, bm_updated_gridded_dataset, floodFillMask) + + # interpolate fields from *preprocessed* composite dataset + # NOTE: while this has already been done in `build_mali_mesh()` + # we are using an updated version of the gridded dataset here, + # which has had unit conversion and extrapolation done. + # Also, it should be assessed if bilinear or + # barycentric used here is preferred for this application. + # Current thinking is they are both equally appropriate. logger.info('calling interpolate_to_mpasli_grid.py') args = ['interpolate_to_mpasli_grid.py', '-s', preprocessed_gridded_dataset, @@ -131,10 +138,16 @@ def run(self): # Now perform bespoke interpolation of geometry and velocity data # from their respective sources - interp_ais_bedmachine(self, data_path, dst_scrip_file, nProcs, - self.mesh_filename) - interp_ais_measures(self, data_path, dst_scrip_file, nProcs, - self.mesh_filename) + 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) diff --git a/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg b/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg index 3ce52b0ee7..74cbae28c0 100644 --- a/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg @@ -52,5 +52,17 @@ use_bed = False # (default value is for Perlmutter) data_path = /global/cfs/cdirs/fanssie/standard_datasets/AIS_datasets +# filename of the BedMachine thickness and bedTopography dataset +# (default value is for Perlmutter) +bedmachine_filename = BedMachineAntarctica_2020-07-15_v02_edits_floodFill_extrap_fillVostok.nc + +# filename of the MEASURES ice velocity dataset +# (default value is for Perlmutter) +measures_filename = antarctica_ice_velocity_450m_v2_edits_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 = ais-bedmap2 + # number of processors to use for ESMF_RegridWeightGen nProcs = 128 diff --git a/compass/landice/tests/greenland/mesh.py b/compass/landice/tests/greenland/mesh.py index 37b66aeb96..9742398a4c 100644 --- a/compass/landice/tests/greenland/mesh.py +++ b/compass/landice/tests/greenland/mesh.py @@ -1,14 +1,14 @@ -from os.path import exists -from shutil import copyfile +import os -import netCDF4 import numpy as np -from mpas_tools.logging import check_call +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, make_region_masks, ) from compass.model import make_graph_file @@ -21,8 +21,8 @@ class Mesh(Step): Attributes ---------- - mesh_type : str - The resolution or mesh type of the test case + mesh_filename : str + File name of the MALI mesh """ def __init__(self, test_case): """ @@ -33,14 +33,19 @@ def __init__(self, test_case): test_case : compass.TestCase The test case this step belongs to - mesh_type : str - The resolution or mesh type of the test case """ super().__init__(test_case=test_case, name='mesh', cpus_per_task=128, min_cpus_per_task=1) + # output files + self.mesh_filename = 'GIS.nc' self.add_output_file(filename='graph.info') - self.add_output_file(filename='GIS.nc') + self.add_output_file(filename=self.mesh_filename) + self.add_output_file( + filename=f'{self.mesh_filename[:-3]}_ismip6_regionMasks.nc') + self.add_output_file( + filename=f'{self.mesh_filename[:-3]}_zwally_regionMasks.nc') + # input files self.add_input_file( filename='greenland_1km_2024_01_29.epsg3413.icesheetonly.nc', target='greenland_1km_2024_01_29.epsg3413.icesheetonly.nc', @@ -56,142 +61,73 @@ def run(self): Run this step of the test case """ logger = self.logger - mesh_name = 'GIS.nc' - section_name = 'mesh' config = self.config - section = config[section_name] - data_path = section.get('data_path') - nProcs = section.get('nProcs') + + section_gis = 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) + + section_name = 'mesh' + + source_gridded_dataset_1km = 'greenland_1km_2024_01_29.epsg3413.icesheetonly.nc' # noqa: E501 + source_gridded_dataset_2km = 'greenland_2km_2024_01_29.epsg3413.nc' 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=source_gridded_dataset_2km, flood_fill_start=[100, 700]) + # Now build the base mesh and perform the standard interpolation build_mali_mesh( self, cell_width, x1, y1, geom_points, geom_edges, - mesh_name=mesh_name, section_name=section_name, - gridded_dataset='greenland_1km_2024_01_29.epsg3413.icesheetonly.nc', # noqa - projection='gis-gimp', geojson_file=None) - - # Create scrip files if they don't already exist - if exists(data_path + '/BedMachineGreenland-v5.scrip.nc'): - logger.info('BedMachine script file exists;' - ' skipping file creation') - else: - logger.info('creating scrip file for BedMachine dataset') - args = ['create_SCRIP_file_from_planar_rectangular_grid.py', - '-i', data_path + '/BedMachineGreenland-v5_edits_floodFill_extrap.nc', # noqa - '-s', data_path + '/BedMachineGreenland-v5.scrip.nc', - '-p', 'gis-gimp', '-r', '2'] - check_call(args, logger=logger) - if exists(data_path + '/greenland_vel_mosaic500.scrip.nc'): - logger.info('Measures script file exists; skipping file creation') - else: - logger.info('creating scrip file for 2006-2010 velocity dataset') - args = ['create_SCRIP_file_from_planar_rectangular_grid.py', - '-i', data_path + '/greenland_vel_mosaic500_extrap.nc', - '-s', data_path + '/greenland_vel_mosaic500.scrip.nc', - '-p', 'gis-gimp', '-r', '2'] - check_call(args, logger=logger) - - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'GIS.nc', '-p', 'gis-gimp'] - check_call(args, logger=logger) + mesh_name=self.mesh_filename, section_name=section_name, + gridded_dataset=source_gridded_dataset_1km, + projection=src_proj, geojson_file=None) + # Create scrip file for the newly generated mesh logger.info('creating scrip file for destination mesh') - scrip_from_mpas('GIS.nc', 'GIS.scrip.nc') - args = ['create_SCRIP_file_from_MPAS_mesh.py', - '-m', 'GIS.nc', - '-s', 'GIS.scrip.nc'] - check_call(args, logger=logger) - - # Create weight files from datasets to mesh - if exists('BedMachine_to_MPAS_weights.nc'): - logger.info('BedMachine_to_MPAS_weights.nc exists; skipping') - else: - logger.info('generating gridded dataset -> MPAS weights') - args = ['srun', '-n', nProcs, 'ESMF_RegridWeightGen', '--source', - data_path + 'BedMachineGreenland-v5.scrip.nc', - '--destination', - 'GIS.scrip.nc', - '--weight', 'BedMachine_to_MPAS_weights.nc', - '--method', 'conserve', - "-i", "-64bit_offset", - "--dst_regional", "--src_regional", '--netcdf4'] - check_call(args, logger=logger) - - if exists('measures_to_MPAS_weights.nc'): - logger.info('measures_to_MPAS_weights.nc exists; skipping') - else: - logger.info('generating gridded dataset -> MPAS weights') - args = ['srun', '-n', nProcs, 'ESMF_RegridWeightGen', '--source', - data_path + 'greenland_vel_mosaic500.scrip.nc', - '--destination', - 'GIS.scrip.nc', - '--weight', 'measures_to_MPAS_weights.nc', - '--method', 'conserve', - "-i", "-64bit_offset", '--netcdf4', - "--dst_regional", "--src_regional", '--ignore_unmapped'] - check_call(args, logger=logger) - - # interpolate fields from BedMachine and Measures - # Using conservative remapping - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - data_path + '/BedMachineGreenland-v5_edits_floodFill_extrap.nc', # noqa - '-d', 'GIS.nc', '-m', 'e', - '-w', 'BedMachine_to_MPAS_weights.nc'] - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - data_path + '/greenland_vel_mosaic500_extrap.nc', - '-d', 'GIS.nc', '-m', 'e', - '-w', 'measures_to_MPAS_weights.nc', - '-v', 'observedSurfaceVelocityX', - 'observedSurfaceVelocityY', - 'observedSurfaceVelocityUncertainty'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'GIS.nc'] - check_call(args, logger=logger) - + 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) + + # create graph file logger.info('creating graph.info') - make_graph_file(mesh_filename=mesh_name, + make_graph_file(mesh_filename=self.mesh_filename, graph_filename='graph.info') - # Create a backup in case clean-up goes awry - copyfile('GIS.nc', 'GIS_backup.nc') - - # Clean up: trim to iceMask and set large velocity - # uncertainties where appropriate. - data = netCDF4.Dataset('GIS.nc', 'r+') - data.set_auto_mask(False) - data.variables['thickness'][:] *= (data.variables['iceMask'][:] > 1.5) - - mask = np.logical_or( - np.isnan( - data.variables['observedSurfaceVelocityUncertainty'][:]), - data.variables['thickness'][:] < 1.0) - mask = np.logical_or( - mask, - data.variables['observedSurfaceVelocityUncertainty'][:] == 0.0) - data.variables['observedSurfaceVelocityUncertainty'][0, mask[0, :]] = 1.0 # noqa # create region masks - mask_filename = f'{mesh_name[:-3]}_ismip6_regionMasks.nc' - make_region_masks(self, mesh_name, mask_filename, + mask_filename = f'{self.mesh_filename[:-3]}_ismip6_regionMasks.nc' + make_region_masks(self, self.mesh_filename, mask_filename, self.cpus_per_task, tags=["Greenland", "ISMIP6", "Shelf"], component='ocean') - mask_filename = f'{mesh_name[:-3]}_zwally_regionMasks.nc' - make_region_masks(self, mesh_name, mask_filename, + mask_filename = f'{self.mesh_filename[:-3]}_zwally_regionMasks.nc' + make_region_masks(self, self.mesh_filename, mask_filename, self.cpus_per_task, tags=['eastCentralGreenland', 'northEastGreenland', @@ -203,18 +139,22 @@ def run(self): 'westCentralGreenland'], all_tags=False) + # Do some final validation of the mesh + ds = xr.open_dataset(self.mesh_filename) # Ensure basalHeatFlux is positive - data.variables['basalHeatFlux'][:] = np.abs( - data.variables['basalHeatFlux'][:]) + ds["basalHeatFlux"] = np.abs(ds.basalHeatFlux) # Ensure reasonable dHdt values - dHdt = data.variables["observedThicknessTendency"][:] - dHdtErr = data.variables["observedThicknessTendencyUncertainty"][:] + dHdt = ds["observedThicknessTendency"] # 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() + # 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') diff --git a/compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg b/compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg index 77ca7e42b5..a285a1cb7c 100644 --- a/compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg @@ -1,9 +1,6 @@ # config options for high_res_mesh test case [mesh] -# path to directory containing BedMachine and Measures datasets -data_path = /global/cfs/cdirs/fanssie/standard_datasets/GIS_datasets/ - # number of levels in the mesh levels = 10 @@ -50,3 +47,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 diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 456dce84ae..ae9736f688 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -492,8 +492,7 @@ Landice Framework mesh.add_bedmachine_thk_to_ais_gridded_data mesh.clean_up_after_interp mesh.gridded_flood_fill - mesh.interp_ais_bedmachine - mesh.interp_ais_measures + mesh.interp_gridded2mali mesh.mpas_flood_fill mesh.preprocess_ais_data mesh.set_rectangular_geom_points_and_edges diff --git a/docs/developers_guide/landice/framework.rst b/docs/developers_guide/landice/framework.rst index 02123c5b22..e3fddd6e17 100644 --- a/docs/developers_guide/landice/framework.rst +++ b/docs/developers_guide/landice/framework.rst @@ -56,13 +56,10 @@ clean up steps after interpolation for the AIS mesh case. :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. -:py:func:`compass.landice.mesh.interp_ais_bedmachine()` interpolates BedMachine -thickness and bedTopography dataset to a MALI mesh, accounting for masking of -the ice extent to avoid interpolation ramps. - -:py:func:`compass.landice.mesh.interp_ais_interp_ais_measures()` interpolates -MEASURES ice velocity dataset to a MALI mesh, accounting for masking at the ice -edge and extrapolation. +:py:func:`compass.landice.mesh.interp_gridded2mali()` interpolates gridded data +(e.g. BedMachine thickness or MEaSUREs ice velocity) to a MALI mesh, accounting +for masking of the ice extent to avoid interpolation ramps. This functions works +for both Antarctica and Greenland. :py:func:`compass.landice.mesh.preprocess_ais_data()` performs adjustments to gridded AIS datasets needed for rest of compass workflow to utilize them. diff --git a/docs/users_guide/landice/test_groups/antarctica.rst b/docs/users_guide/landice/test_groups/antarctica.rst index f0a0629e80..cef797f9ee 100644 --- a/docs/users_guide/landice/test_groups/antarctica.rst +++ b/docs/users_guide/landice/test_groups/antarctica.rst @@ -73,7 +73,19 @@ the mesh generation options are adjusted through the config file. # path to directory containing BedMachine and Measures datasets # (default value is for Perlmutter) data_path = /global/cfs/cdirs/fanssie/standard_datasets/AIS_datasets - + + # filename of the BedMachine thickness and bedTopography dataset + # (default value is for Perlmutter) + bedmachine_filename = BedMachineAntarctica_2020-07-15_v02_edits_floodFill_extrap_fillVostok.nc + + # filename of the MEASURES ice velocity dataset + # (default value is for Perlmutter) + measures_filename = antarctica_ice_velocity_450m_v2_edits_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 = ais-bedmap2 + # number of processors to use for ESMF_RegridWeightGen nProcs = 128 diff --git a/docs/users_guide/landice/test_groups/greenland.rst b/docs/users_guide/landice/test_groups/greenland.rst index 902baa0161..146119b909 100644 --- a/docs/users_guide/landice/test_groups/greenland.rst +++ b/docs/users_guide/landice/test_groups/greenland.rst @@ -86,6 +86,26 @@ The other test cases do not use config options. 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 + smoke_test ----------