diff --git a/compass/ocean/mesh/cull.py b/compass/ocean/mesh/cull.py index 81965f2fcc..2ff60f8496 100644 --- a/compass/ocean/mesh/cull.py +++ b/compass/ocean/mesh/cull.py @@ -11,6 +11,7 @@ from mpas_tools.logging import LoggingContext, check_call from mpas_tools.mesh.conversion import cull from mpas_tools.mesh.creation.sort_mesh import sort_mesh +from mpas_tools.mesh.cull import cull_dataset, map_culled_to_base from mpas_tools.mesh.mask import compute_mpas_flood_fill_mask from mpas_tools.ocean import inject_bathymetry from mpas_tools.ocean.coastline_alteration import ( @@ -126,7 +127,8 @@ def setup(self): work_dir_target=target) self.add_output_file('topography_culled.nc') self.add_output_file('map_culled_to_base.nc') - + self.add_input_file(filename='south_pole.geojson', + package='compass.ocean.mesh') config = self.config self.cpus_per_task = config.getint('spherical_mesh', 'cull_mesh_cpus_per_task') @@ -268,39 +270,6 @@ def cull_mesh(with_cavities=False, with_critical_passages=False, convert_to_cdf5, latitude_threshold, sweep_count) -def write_map_culled_to_base(base_mesh_filename, culled_mesh_filename, - out_filename): - ds_base = xr.open_dataset(base_mesh_filename) - ncells_base = ds_base.sizes['nCells'] - lon_base = ds_base.lonCell.values - lat_base = ds_base.latCell.values - - ds_culled = xr.open_dataset(culled_mesh_filename) - ncells_culled = ds_culled.sizes['nCells'] - lon_culled = ds_culled.lonCell.values - lat_culled = ds_culled.latCell.values - - # create a map from lat-lon pairs to base-mesh cell indices - map_base = dict() - for cell_index in range(ncells_base): - lon = lon_base[cell_index] - lat = lat_base[cell_index] - map_base[(lon, lat)] = cell_index - - # create a map from culled- to base-mesh cell indices - map_culled_to_base = list() - for cell_index in range(ncells_culled): - lon = lon_culled[cell_index] - lat = lat_culled[cell_index] - # each (lon, lat) for a culled cell *must* be in the base mesh - map_culled_to_base.append(map_base[(lon, lat)]) - - ds_map_culled_to_base = xr.Dataset() - ds_map_culled_to_base['mapCulledToBase'] = (('nCells',), - map_culled_to_base) - write_netcdf(ds_map_culled_to_base, out_filename) - - def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages, custom_critical_passages, custom_land_blockages, preserve_floodplain, use_progress_bar, @@ -321,6 +290,9 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages, netcdf_engine = mpas_tools.io.default_engine has_remapped_topo = os.path.exists('topography.nc') + if with_cavities and not has_remapped_topo: + raise ValueError('Mesh culling with caviites must be from ' + 'remapped topography.') if has_remapped_topo: _land_mask_from_topo(with_cavities, @@ -463,24 +435,11 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages, check_call(args, logger=logger) if has_remapped_topo: - write_map_culled_to_base(base_mesh_filename='base_mesh.nc', - culled_mesh_filename='culled_mesh.nc', - out_filename='map_culled_to_base.nc') - _cull_topo() + _cull_topo(with_cavities, process_count, logger) if with_cavities: - if has_remapped_topo: - _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, - logger=logger, - mesh_filename='culled_mesh.nc', - geojson_filename='ice_coverage.geojson', - mask_filename='ice_coverage.nc') - - dsMask = xr.open_dataset('ice_coverage.nc') + dsMask = xr.open_dataset('topography_culled.nc') + dsMask = dsMask[['regionCellMasks',]] dsMask = add_land_locked_cells_to_mask( dsMask, dsCulledMesh, latitude_threshold=latitude_threshold, @@ -511,13 +470,26 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages, use_progress_bar=use_progress_bar) -def _cull_topo(): - ds_map = xr.open_dataset('map_culled_to_base.nc') - map_culled_to_base = ds_map.mapCulledToBase.values +def _cull_topo(with_cavities, process_count, logger): ds_topo = xr.open_dataset('topography.nc') - ds_topo = ds_topo.isel(nCells=map_culled_to_base) - write_netcdf(ds_topo, 'topography_culled.nc') + ds_base = xr.open_dataset('base_mesh.nc') + if with_cavities: + _add_land_ice_mask_and_mask_draft(ds_topo, ds_base, logger) + write_netcdf(ds_topo, 'topography_with_land_ice_mask.nc') + + ds_culled = xr.open_dataset('culled_mesh.nc') + ds_map_culled_to_base = map_culled_to_base(ds_base=ds_base, + ds_culled=ds_culled, + workers=process_count) + logger.info('Culling topography') + write_netcdf(ds_map_culled_to_base, 'map_culled_to_base.nc') + ds_culled_topo = cull_dataset(ds=ds_topo, ds_base_mesh=ds_base, + ds_culled_mesh=ds_culled, + ds_map_culled_to_base=ds_map_culled_to_base, + logger=logger) + + write_netcdf(ds_culled_topo, 'topography_culled.nc') def _land_mask_from_topo(with_cavities, topo_filename, mask_filename): @@ -544,18 +516,49 @@ 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) +def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh, logger): 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) + ocean_mask = 1 - land_ice_mask - ds_mask = xr.Dataset() - ds_mask['regionCellMasks'] = land_ice_mask - write_netcdf(ds_mask, mask_filename) + gf = GeometricFeatures() + fc_ocean_seed = gf.read(componentName='ocean', objectType='point', + tags=['seed_point']) + + # flood fill the ocean portion to fill in any holes in the land ice + ds_mask = compute_mpas_flood_fill_mask(dsMesh=ds_base_mesh, + daGrow=ocean_mask, + fcSeed=fc_ocean_seed, + logger=logger) + land_ice_mask = 1 - ds_mask.cellSeedMask + + fc_south_pole_seed = read_feature_collection('south_pole.geojson') + + # now flood fill the ice portion to remove isolated land ice + ds_mask = compute_mpas_flood_fill_mask(dsMesh=ds_base_mesh, + daGrow=land_ice_mask, + fcSeed=fc_south_pole_seed, + logger=logger) + land_ice_mask = ds_mask.cellSeedMask + ds_topo['landIceMask'] = land_ice_mask + 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_mask * ds_topo.landIceDraftObserved) def _land_mask_from_geojson(with_cavities, process_count, logger, diff --git a/compass/ocean/mesh/south_pole.geojson b/compass/ocean/mesh/south_pole.geojson new file mode 100644 index 0000000000..aa70859c34 --- /dev/null +++ b/compass/ocean/mesh/south_pole.geojson @@ -0,0 +1,22 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "name": "South Pole", + "tags": "seed_point", + "object": "point", + "component": "landice", + "author": "Xylar Asay-Davis" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 0.0, + -90.0 + ] + } + } + ] +} diff --git a/compass/ocean/tests/global_ocean/init/remap_ice_shelf_melt.py b/compass/ocean/tests/global_ocean/init/remap_ice_shelf_melt.py index f92f2fd54a..4980b81acf 100644 --- a/compass/ocean/tests/global_ocean/init/remap_ice_shelf_melt.py +++ b/compass/ocean/tests/global_ocean/init/remap_ice_shelf_melt.py @@ -4,10 +4,10 @@ import xarray as xr from mpas_tools.cime.constants import constants from mpas_tools.io import write_netcdf +from mpas_tools.mesh.cull import write_map_culled_to_base from pyremap import MpasCellMeshDescriptor, ProjectionGridDescriptor, Remapper from scipy.spatial import KDTree -from compass.ocean.mesh.cull import write_map_culled_to_base from compass.step import Step @@ -445,7 +445,7 @@ def _land_ice_mask_on_base_mesh(base_mesh_filename, land_ice_mask_filename, """ Map the land-ice mask back to the base mesh """ ds_map = xr.open_dataset(map_culled_to_base_filename) - map_culled_to_base = ds_map.mapCulledToBase.values + map_culled_to_base = ds_map.mapCulledToBaseCell.values ds_base = xr.open_dataset(base_mesh_filename) ncells_base = ds_base.sizes['nCells'] @@ -513,7 +513,7 @@ def _reroute_missing_flux(base_mesh_filename, map_culled_to_base_filename, _, indices = tree.query(reroute_xyz, workers=-1) ds_map = xr.open_dataset(map_culled_to_base_filename) - map_culled_to_base = ds_map.mapCulledToBase.values + map_culled_to_base = ds_map.mapCulledToBaseCell.values fwf_rerouted = land_ice_mask.values * fwf.where(fwf.notnull(), 0.).values hf_rerouted = land_ice_mask.values * hf.where(hf.notnull(), 0.).values diff --git a/compass/version.py b/compass/version.py index 864ac90c1a..30df5b3418 100644 --- a/compass/version.py +++ b/compass/version.py @@ -1 +1 @@ -__version__ = '1.3.0-alpha.1' +__version__ = '1.3.0-alpha.2' diff --git a/conda/compass_env/spec-file.template b/conda/compass_env/spec-file.template index 7435830a3e..90c6d6d074 100644 --- a/conda/compass_env/spec-file.template +++ b/conda/compass_env/spec-file.template @@ -22,7 +22,7 @@ matplotlib-base metis moab >=5.5.1 moab=*={{ mpi_prefix }}_tempest_* -mpas_tools=0.32.0 +mpas_tools=0.33.0 nco netcdf4=*=nompi_* numpy diff --git a/conda/recipe/meta.yaml b/conda/recipe/meta.yaml index 7cef5e05d6..0029f52ec0 100644 --- a/conda/recipe/meta.yaml +++ b/conda/recipe/meta.yaml @@ -1,5 +1,5 @@ {% set name = "compass" %} -{% set version = "1.2.0alpha.8" %} +{% set version = "1.3.0alpha.2" %} {% set build = 0 %} {% if mpi == "nompi" %} @@ -57,7 +57,7 @@ requirements: - mache 1.17.0 - matplotlib-base - metis - - mpas_tools 0.26.0 + - mpas_tools 0.33.0 - nco - netcdf4 * nompi_* - numpy