From 1a9ec0e866afad10aeee52a619c970781ad7049b Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 19 Nov 2023 15:09:15 -0700 Subject: [PATCH] Fix a bug in how `landIceMask` is computed Before this fix, cells that were more than half land according to the remapped topography dataset but had been included in the ocean (through critical passages) were mistakenly being categorized as land ice simply because they were "not ocean". With this fix, the `landIceMask` is computed based on where the `landIceFracObserved` from the remapped topography dataset is greater than 0.5. A new function is needed to compute this mask because this definition is not the same as the mask used to cull the mesh when ice-shelf cavities are excluded from the ocean domain. --- compass/ocean/mesh/cull.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/compass/ocean/mesh/cull.py b/compass/ocean/mesh/cull.py index d9989ac598..b7a015425f 100644 --- a/compass/ocean/mesh/cull.py +++ b/compass/ocean/mesh/cull.py @@ -413,9 +413,8 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages, if with_cavities: if has_remapped_topo: - _land_mask_from_topo(with_cavities=False, - topo_filename='topography_culled.nc', - mask_filename='ice_coverage.nc') + _land_ice_mask_from_topo(topo_filename='topography_culled.nc', + mask_filename='ice_coverage.nc') else: _land_mask_from_geojson(with_cavities=False, process_count=process_count, @@ -506,6 +505,20 @@ def _land_mask_from_topo(with_cavities, topo_filename, mask_filename): write_netcdf(ds_mask, mask_filename) +def _land_ice_mask_from_topo(topo_filename, mask_filename): + ds_topo = xr.open_dataset(topo_filename) + 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_mask = land_ice_mask.expand_dims(dim='nRegions', axis=1) + + ds_mask = xr.Dataset() + ds_mask['regionCellMasks'] = land_ice_mask + write_netcdf(ds_mask, mask_filename) + + def _land_mask_from_geojson(with_cavities, process_count, logger, mesh_filename, geojson_filename, mask_filename): gf = GeometricFeatures()