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

@@ -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()