From 1d6ef70e9dd7ecf4fb01fd627cfd05cac4f8d466 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 30 Jul 2024 04:01:22 -0500 Subject: [PATCH] Floodfill landIceFraction, not landIceMask During the `cull_mesh` step of creating a global ocean mesh, if ice-shelf cavities are present, we need to flood fill to make sure all land ice is connected (no islands of isolated land ice). This merge changes the flood fill to be on all locations where `landIceFraction > 0`, not all locations where `landIceFraction > 0.5`. It then masks all `landIce*` fields to only be nonzero within the flood-filled region. This prevents islands of nonzero `landIceFraction`, `landIcePressure`, `landIceDraft`, etc. --- compass/ocean/mesh/cull.py | 38 +++++++++++++++++--------------------- 1 file changed, 17 insertions(+), 21 deletions(-) diff --git a/compass/ocean/mesh/cull.py b/compass/ocean/mesh/cull.py index a974499dcd..425cfc7cd1 100644 --- a/compass/ocean/mesh/cull.py +++ b/compass/ocean/mesh/cull.py @@ -488,7 +488,7 @@ def _cull_topo(with_cavities, process_count, logger, latitude_threshold, write_netcdf(ds_map_culled_to_base, 'map_culled_to_base.nc') if with_cavities: - _add_land_ice_mask_and_mask_draft(ds_topo, ds_base, + _flood_fill_and_add_land_ice_mask(ds_topo, ds_base, ds_map_culled_to_base, ds_preserve, logger, latitude_threshold, sweep_count) @@ -525,7 +525,7 @@ def _land_mask_from_topo(with_cavities, topo_filename, mask_filename): write_netcdf(ds_mask, mask_filename) -def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh, +def _flood_fill_and_add_land_ice_mask(ds_topo, ds_base_mesh, ds_map_culled_to_base, ds_perserve, logger, latitude_threshold, sweep_count): @@ -534,7 +534,8 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh, not_preserve = ds.transectCellMasks.sum(dim='nTransects') == 0 for var in ['landIceFracObserved', 'landIcePressureObserved', - 'landIceDraftObserved', 'landIceGroundedFracObserved', + 'landIceDraftObserved', 'ssh', + 'landIceGroundedFracObserved', 'landIceFloatingFracObserved', 'landIceThkObserved']: ds_topo[var] = ds_topo[var].where(not_preserve, 0.0) @@ -543,8 +544,7 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh, land_ice_frac = ds_topo.landIceFracObserved - # we want the mask to be 1 where there's at least half land-ice - land_ice_mask = xr.where(land_ice_frac > 0.5, 1, 0) + land_ice_present = xr.where(land_ice_frac > 0.0, 1, 0) gf = GeometricFeatures() fc_ocean_seed = gf.read(componentName='ocean', objectType='point', @@ -552,13 +552,22 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh, fc_south_pole_seed = read_feature_collection('south_pole.geojson') - # flood fill the ice portion to remove isolated land ice + # flood fill anywhere land ice is present to remove isolated land ice ds_mask = compute_mpas_flood_fill_mask(dsMesh=ds_base_mesh, - daGrow=land_ice_mask, + daGrow=land_ice_present, fcSeed=fc_south_pole_seed, logger=logger) - land_ice_mask = ds_mask.cellSeedMask + land_ice_present = ds_mask.cellSeedMask + + # update land-ice variables and ocean fraction accordingly + for var in ['landIceFracObserved', 'landIcePressureObserved', + 'landIceDraftObserved', 'ssh', 'landIceGroundedFracObserved', + 'landIceFloatingFracObserved', 'landIceThkObserved']: + ds_topo[var] = ds_topo[var].where(land_ice_present, 0.0) + + land_ice_frac = ds_topo.landIceFracObserved + land_ice_mask = xr.where(land_ice_frac > 0.5, 1, 0) # now, remove land-locked or land-ice-locked cells @@ -600,19 +609,6 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh, region_cell_mask = land_ice_mask.expand_dims(dim='nRegions', axis=1) ds_topo['regionCellMasks'] = region_cell_mask - # we also want to mask out the land-ice draft for cells detatched from the - # ice sheet, this time anywhere the land-ice fraction > 0 - land_ice_draft_mask = xr.where(land_ice_frac > 0.0, 1, 0) - ds_mask = compute_mpas_flood_fill_mask(dsMesh=ds_base_mesh, - daGrow=land_ice_draft_mask, - fcSeed=fc_south_pole_seed, - logger=logger) - - land_ice_draft_mask = ds_mask.cellSeedMask - ds_topo['landIceDraftMask'] = land_ice_draft_mask - ds_topo['landIceDraftObserved'] = ( - land_ice_draft_mask * ds_topo.landIceDraftObserved) - def _land_mask_from_geojson(with_cavities, process_count, logger, mesh_filename, geojson_filename, mask_filename):