From 2adc99be678761463c1b369821f5a166641c26fe Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 7 Apr 2023 12:45:31 -0700 Subject: [PATCH 01/19] Move build_cell_width() to compass.landice.mesh Move build_cell_width() to compass.landice.mesh to greatly reduce code redundancy between all mesh-generating test cases. --- compass/landice/mesh.py | 86 ++++++++++++++++--- .../landice/tests/antarctica/antarctica.cfg | 6 ++ compass/landice/tests/antarctica/mesh.py | 62 +------------ .../greenland/high_res_mesh/high_res_mesh.cfg | 8 ++ compass/landice/tests/greenland/mesh.py | 75 ++-------------- compass/landice/tests/humboldt/humboldt.cfg | 6 ++ compass/landice/tests/humboldt/mesh.py | 68 ++------------- compass/landice/tests/kangerlussuaq/mesh.py | 68 ++------------- .../tests/kangerlussuaq/mesh_gen/mesh_gen.cfg | 6 ++ compass/landice/tests/koge_bugt_s/mesh.py | 78 +++-------------- .../tests/koge_bugt_s/mesh_gen/mesh_gen.cfg | 17 ++++ .../thwaites/high_res_mesh/high_res_mesh.cfg | 19 ++++ compass/landice/tests/thwaites/mesh.py | 67 ++------------- 13 files changed, 180 insertions(+), 386 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index a839525059..c254900d4b 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -2,6 +2,7 @@ import jigsawpy import numpy as np +from netCDF4 import Dataset def gridded_flood_fill(field, iStart=None, jStart=None): @@ -167,16 +168,16 @@ def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, section = self.config[section] # Get config inputs for cell spacing functions - min_spac = float(section.get('min_spac')) - max_spac = float(section.get('max_spac')) - high_log_speed = float(section.get('high_log_speed')) - low_log_speed = float(section.get('low_log_speed')) - high_dist = float(section.get('high_dist')) - low_dist = float(section.get('low_dist')) - high_dist_bed = float(section.get('high_dist_bed')) - low_dist_bed = float(section.get('low_dist_bed')) - low_bed = float(section.get('low_bed')) - high_bed = float(section.get('high_bed')) + min_spac = section.getfloat('min_spac') + max_spac = section.getfloat('max_spac') + high_log_speed = section.getfloat('high_log_speed') + low_log_speed = section.getfloat('low_log_speed') + high_dist = section.getfloat('high_dist') + low_dist = section.getfloat('low_dist') + high_dist_bed = section.getfloat('high_dist_bed') + low_dist_bed = section.getfloat('low_dist_bed') + low_bed = section.getfloat('low_bed') + high_bed = section.getfloat('high_bed') # convert km to m cull_distance = float(section.get('cull_distance')) * 1.e3 @@ -411,3 +412,68 @@ def get_dist_to_edge_and_GL(self, thk, topg, x, y, section, window_size=None): 'seconds'.format(toc - tic)) return dist_to_edge, dist_to_grounding_line + + +def build_cell_width(self, section_name, gridded_dataset, + flood_fill_start=[None, None]): + """ + Determine MPAS mesh cell size based on user-defined density function. + + """ + + section = self.config[section_name] + # get needed fields from gridded dataset + f = Dataset(gridded_dataset, 'r') + f.set_auto_mask(False) # disable masked arrays + + x1 = f.variables['x1'][:] + y1 = f.variables['y1'][:] + thk = f.variables['thk'][0, :, :] + topg = f.variables['topg'][0, :, :] + vx = f.variables['vx'][0, :, :] + vy = f.variables['vy'][0, :, :] + + f.close() + + # Define extent of region to mesh. + xx0 = section.get('x_min') + xx1 = section.get('x_max') + yy0 = section.get('y_min') + yy1 = section.get('y_max') + + if 'None' in [xx0, xx1, yy0, yy1]: + xx0 = np.min(x1) + xx1 = np.max(x1) + yy0 = np.min(y1) + yy1 = np.max(y1) + else: + xx0 = float(xx0) + xx1 = float(xx1) + yy0 = float(yy0) + yy1 = float(yy1) + + geom_points, geom_edges = set_rectangular_geom_points_and_edges( + xx0, xx1, yy0, yy1) + + # Remove ice not connected to the ice sheet. + floodMask = gridded_flood_fill(thk) + thk[floodMask == 0] = 0.0 + vx[floodMask == 0] = 0.0 + vy[floodMask == 0] = 0.0 + + # Calculate distance from each grid point to ice edge + # and grounding line, for use in cell spacing functions. + distToEdge, distToGL = get_dist_to_edge_and_GL( + self, thk, topg, x1, + y1, section=section_name) + + # Set cell widths based on mesh parameters set in config file + cell_width = set_cell_width(self, section=section_name, + thk=thk, bed=topg, vx=vx, vy=vy, + dist_to_edge=distToEdge, + dist_to_grounding_line=distToGL, + flood_fill_iStart=flood_fill_start[0], + flood_fill_jStart=flood_fill_start[1]) + + return (cell_width.astype('float64'), x1.astype('float64'), + y1.astype('float64'), geom_points, geom_edges, floodMask) diff --git a/compass/landice/tests/antarctica/antarctica.cfg b/compass/landice/tests/antarctica/antarctica.cfg index 6bb0910815..f3d7cb7d36 100644 --- a/compass/landice/tests/antarctica/antarctica.cfg +++ b/compass/landice/tests/antarctica/antarctica.cfg @@ -4,6 +4,12 @@ # number of levels in the mesh levels = 5 +# Bounds of Antarctic mesh +x_min = -3333500. +x_max = 3330500. +y_min = -3333500. +y_max = 3330500. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index 5639d3c559..d1d636ef45 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/compass/landice/tests/antarctica/mesh.py @@ -10,12 +10,7 @@ from mpas_tools.mesh.conversion import convert, cull from mpas_tools.mesh.creation import build_planar_mesh -from compass.landice.mesh import ( - get_dist_to_edge_and_GL, - gridded_flood_fill, - set_cell_width, - set_rectangular_geom_points_and_edges, -) +from compass.landice.mesh import build_cell_width from compass.model import make_graph_file from compass.step import Step @@ -67,7 +62,9 @@ def run(self): logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodFillMask = \ - self.build_cell_width() + build_cell_width( + self, section_name='antarctica', + gridded_dataset='antarctica_8km_2020_10_20.nc') logger.info('calling build_planar_mesh') build_planar_mesh(cell_width, x1, y1, geom_points, geom_edges, logger=logger) @@ -172,57 +169,6 @@ def run(self): self.cpus_per_task, tags=['ISMIP6_Basin']) - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. - """ - # get needed fields from AIS dataset - f = netCDF4.Dataset('antarctica_8km_2020_10_20.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Antarctica mesh. - xx0 = -3333500 - xx1 = 3330500 - yy0 = -3333500 - yy1 = 3330500 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL(self, thk, topg, x1, y1, - section='antarctica', - window_size=1.e5) - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, section='antarctica', thk=thk, - vx=vx, vy=vy, dist_to_edge=distToEdge, - dist_to_grounding_line=distToGL) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges, floodMask) - def _make_region_masks(self, mesh_filename, mask_filename, cores, tags): logger = self.logger gf = GeometricFeatures() diff --git a/compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg b/compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg index 8ad755ba30..e77869bd89 100644 --- a/compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg +++ b/compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg @@ -4,6 +4,14 @@ # number of levels in the mesh levels = 10 +# Bounds of GIS mesh. If any bound is set +# to None, the mesh will use the full extent +# of the gridded dataset. +x_min = None +x_max = None +y_min = None +y_max = None + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. diff --git a/compass/landice/tests/greenland/mesh.py b/compass/landice/tests/greenland/mesh.py index d3851b57c2..2f8559f734 100644 --- a/compass/landice/tests/greenland/mesh.py +++ b/compass/landice/tests/greenland/mesh.py @@ -1,17 +1,10 @@ -import netCDF4 -import numpy as np import xarray from mpas_tools.io import write_netcdf 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 compass.landice.mesh import ( - get_dist_to_edge_and_GL, - gridded_flood_fill, - set_cell_width, - set_rectangular_geom_points_and_edges, -) +from compass.landice.mesh import build_cell_width from compass.model import make_graph_file from compass.step import Step @@ -59,8 +52,12 @@ def run(self): config = self.config section = config['high_res_GIS_mesh'] - logger.info('calling build_cell_wdith') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() + logger.info('calling build_cell_width') + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name='high_res_GIS_mesh', + gridded_dataset='greenland_2km_2020_04_20.epsg3413.nc', + flood_fill_start=[100, 700]) logger.info('calling build_planar_mesh') build_planar_mesh(cell_width, x1, y1, geom_points, geom_edges, logger=logger) @@ -136,61 +133,3 @@ def run(self): logger.info('creating graph.info') make_graph_file(mesh_filename='GIS.nc', graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes the definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. - """ - # get needed fields from GIS dataset - f = netCDF4.Dataset('greenland_2km_2020_04_20.epsg3413.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the GIS mesh. - xx0 = np.min(x1) - xx1 = np.max(x1) - yy0 = np.min(y1) - yy1 = np.max(y1) - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL( - self, thk, topg, x1, y1, section='high_res_GIS_mesh') - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file. - # Start flood-fill for bed topography in the ocean at a point - # with bed < low_bed. This is often necessary to remove pockets - # of high resolution that are no longer connected to the rest of - # the high resolution part of the mesh. These indices are specific - # to the mesh and the gridded data product used to set cell spacing. - cell_width = set_cell_width(self, section='high_res_GIS_mesh', thk=thk, - bed=topg, vx=vx, vy=vy, - dist_to_edge=distToEdge, - dist_to_grounding_line=distToGL, - flood_fill_iStart=100, - flood_fill_jStart=700) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges) diff --git a/compass/landice/tests/humboldt/humboldt.cfg b/compass/landice/tests/humboldt/humboldt.cfg index 9bac801e01..59db5bc1b6 100644 --- a/compass/landice/tests/humboldt/humboldt.cfg +++ b/compass/landice/tests/humboldt/humboldt.cfg @@ -4,6 +4,12 @@ # number of levels in the mesh levels = 10 +# Bounds of Humboldt domain +x_min = -630000. +x_max = 84000. +y_min = -1560000. +y_max = -860000. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. diff --git a/compass/landice/tests/humboldt/mesh.py b/compass/landice/tests/humboldt/mesh.py index 5d725aaa54..5ee18a8c46 100644 --- a/compass/landice/tests/humboldt/mesh.py +++ b/compass/landice/tests/humboldt/mesh.py @@ -1,16 +1,10 @@ -import netCDF4 import xarray from mpas_tools.io import write_netcdf 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 compass.landice.mesh import ( - get_dist_to_edge_and_GL, - gridded_flood_fill, - set_cell_width, - set_rectangular_geom_points_and_edges, -) +from compass.landice.mesh import build_cell_width from compass.model import make_graph_file from compass.step import Step @@ -62,8 +56,11 @@ def run(self): config = self.config section = config['humboldt_mesh'] - logger.info('calling build_cell_wdith') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() + logger.info('calling build_cell_width') + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name='humboldt_mesh', + gridded_dataset='greenland_2km_2020_04_20.epsg3413.nc') logger.info('calling build_planar_mesh') build_planar_mesh(cell_width, x1, y1, geom_points, geom_edges, logger=logger) @@ -166,56 +163,3 @@ def run(self): logger.info('creating graph.info') make_graph_file(mesh_filename='Humboldt_1to10km.nc', graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. In the future, this function - and its components will likely be separated into separate generalized - functions to be reusable by multiple test groups. - """ - # get needed fields from GIS dataset - f = netCDF4.Dataset('greenland_2km_2020_04_20.epsg3413.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Humboldt Glacier mesh. - xx0 = -630000 - xx1 = 84000 - yy0 = -1560000 - yy1 = -860000 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL( - self, thk, topg, x1, - y1, section='humboldt_mesh') - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, section='humboldt_mesh', - thk=thk, bed=topg, vx=vx, vy=vy, - dist_to_edge=distToEdge, - dist_to_grounding_line=distToGL) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges) diff --git a/compass/landice/tests/kangerlussuaq/mesh.py b/compass/landice/tests/kangerlussuaq/mesh.py index 8cd4b3c8a0..b126fcc63f 100644 --- a/compass/landice/tests/kangerlussuaq/mesh.py +++ b/compass/landice/tests/kangerlussuaq/mesh.py @@ -1,16 +1,10 @@ -import netCDF4 import xarray from mpas_tools.io import write_netcdf 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 compass.landice.mesh import ( - get_dist_to_edge_and_GL, - gridded_flood_fill, - set_cell_width, - set_rectangular_geom_points_and_edges, -) +from compass.landice.mesh import build_cell_width from compass.model import make_graph_file from compass.step import Step @@ -63,8 +57,11 @@ def run(self): config = self.config section = config['high_res_Kangerlussuaq_mesh'] - logger.info('calling build_cell_wdith') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() + logger.info('calling build_cell_width') + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name='high_res_Kangerlussuaq_mesh', + gridded_dataset='greenland_8km_2020_04_20.epsg3413.nc') logger.info('calling build_planar_mesh') build_planar_mesh(cell_width, x1, y1, geom_points, geom_edges, logger=logger) @@ -162,56 +159,3 @@ def run(self): logger.info('creating graph.info') make_graph_file(mesh_filename='Kangerlussuaq.nc', graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. - """ - # get needed fields from GIS dataset - f = netCDF4.Dataset('greenland_8km_2020_04_20.epsg3413.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Kangerlussuaq Glacier mesh. - xx0 = 215000 - xx1 = 544000 - yy0 = -2370000 - yy1 = -2070000 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL( - self, thk, topg, x1, - y1, section='high_res_Kangerlussuaq_mesh') - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, - section='high_res_Kangerlussuaq_mesh', - thk=thk, bed=topg, vx=vx, vy=vy, - dist_to_edge=distToEdge, - dist_to_grounding_line=distToGL) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges) diff --git a/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg b/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg index 7db6fd63a6..355c4d01d5 100644 --- a/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg @@ -4,6 +4,12 @@ # number of levels in the mesh levels = 10 +# Bounds of Kangerlussuaq domain +x_min = 215000. +x_max = 544000. +y_min = -2370000. +y_max = -2070000. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. diff --git a/compass/landice/tests/koge_bugt_s/mesh.py b/compass/landice/tests/koge_bugt_s/mesh.py index 9a6fa5442b..775762508c 100644 --- a/compass/landice/tests/koge_bugt_s/mesh.py +++ b/compass/landice/tests/koge_bugt_s/mesh.py @@ -1,19 +1,12 @@ -import numpy as np -import netCDF4 import xarray -from matplotlib import pyplot as plt - -from mpas_tools.mesh.creation import build_planar_mesh -from mpas_tools.mesh.conversion import convert, cull -from mpas_tools.planar_hex import make_planar_hex_mesh from mpas_tools.io import write_netcdf 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 compass.step import Step +from compass.landice.mesh import build_cell_width from compass.model import make_graph_file -from compass.landice.mesh import gridded_flood_fill, \ - set_rectangular_geom_points_and_edges, \ - set_cell_width, get_dist_to_edge_and_GL +from compass.step import Step class Mesh(Step): @@ -43,9 +36,9 @@ def __init__(self, test_case): self.add_output_file(filename='graph.info') self.add_output_file(filename='Koge_Bugt_S.nc') self.add_input_file( - filename='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - target='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - database='') + filename='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', + target='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', + database='') self.add_input_file(filename='Koge_Bugt_S.geojson', package='compass.landice.tests.koge_bugt_s', target='Koge_Bugt_S.geojson', @@ -65,7 +58,10 @@ def run(self): section = config['high_res_KogeBugtS_mesh'] logger.info('calling build_cell_width') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name='high_res_KogeBugtS_mesh', + gridded_dataset='greenland_8km_2020_04_20.epsg3413.nc') logger.info('calling build_planar_mesh') build_planar_mesh(cell_width, x1, y1, geom_points, geom_edges, logger=logger) @@ -163,55 +159,3 @@ def run(self): logger.info('creating graph.info') make_graph_file(mesh_filename='Koge_Bugt_S.nc', graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. - """ - # get needed fields from GIS dataset - f = netCDF4.Dataset('greenland_8km_2020_04_20.epsg3413.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Koge_Bugt_S Glacier mesh. - xx0 = 0 - xx1 = 225000 - yy0 = -2830000 - yy1 = -2700000 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL(self, thk, topg, x1, - y1, window_size=1.e5) - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, - section='high_res_KogeBugtS_mesh', - thk=thk, vx=vx, vy=vy, - dist_to_edge=distToEdge, - dist_to_grounding_line=None) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges) diff --git a/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg b/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg index 4631a34063..1b38c73804 100644 --- a/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg @@ -4,6 +4,12 @@ # number of levels in the mesh levels = 10 +# Bounds of Koge Bugt S mesh +x_min = 0. +x_max = 225000. +y_min = -2830000. +y_max = -2700000. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. @@ -23,7 +29,18 @@ high_dist = 1.e5 # distance within which cell spacing = min_spac (meters) low_dist = 1.e4 +# mesh density parameters used if use_bed = True +# distance at which bed topography has no effect +high_dist_bed = 1.e5 +# distance within which bed topography has maximum effect +low_dist_bed = 5.e4 +# Bed elev beneath which cell spacing is minimized +low_bed = 50.0 +# Bed elev above which cell spacing is maximized +high_bed = 100.0 + # mesh density functions use_speed = True use_dist_to_grounding_line = False use_dist_to_edge = True +use_bed = False diff --git a/compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg b/compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg index dcb2ac81a0..b5336c9753 100644 --- a/compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg +++ b/compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg @@ -4,6 +4,12 @@ # number of levels in the mesh levels = 10 +# Bounds of Thwaites regional mesh +x_min = -1864434. +x_max = -975432. +y_min = -901349. +y_max = 0. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. @@ -23,7 +29,20 @@ high_dist = 1.e5 # distance within which cell spacing = min_spac low_dist = 5.e4 +# mesh density parameters used if use_bed = True +# These settings are taken from the Humboldt mesh +# and have not yet been evaluated for Thwaites. +# distance at which bed topography has no effect +high_dist_bed = 1.e5 +# distance within which bed topography has maximum effect +low_dist_bed = 5.e4 +# Bed elev beneath which cell spacing is minimized +low_bed = 50.0 +# Bed elev above which cell spacing is maximized +high_bed = 100.0 + # mesh density functions use_speed = True use_dist_to_grounding_line = True use_dist_to_edge = False +use_bed = False diff --git a/compass/landice/tests/thwaites/mesh.py b/compass/landice/tests/thwaites/mesh.py index e0cb2da985..5a35383151 100644 --- a/compass/landice/tests/thwaites/mesh.py +++ b/compass/landice/tests/thwaites/mesh.py @@ -1,16 +1,12 @@ -import netCDF4 import xarray - -from mpas_tools.mesh.creation import build_planar_mesh -from mpas_tools.mesh.conversion import convert, cull from mpas_tools.io import write_netcdf 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 compass.step import Step +from compass.landice.mesh import build_cell_width from compass.model import make_graph_file -from compass.landice.mesh import gridded_flood_fill, \ - set_rectangular_geom_points_and_edges, \ - set_cell_width, get_dist_to_edge_and_GL +from compass.step import Step class Mesh(Step): @@ -51,7 +47,10 @@ def run(self): section = config['high_res_mesh'] logger.info('calling build_cell_width') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name='high_res_mesh', + gridded_dataset='antarctica_8km_2020_10_20.nc') logger.info('calling build_planar_mesh') build_planar_mesh(cell_width, x1, y1, geom_points, geom_edges, logger=logger) @@ -155,53 +154,3 @@ def run(self): logger.info('creating graph.info') make_graph_file(mesh_filename='Thwaites.nc', graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. In the future, this function - and its components will likely be separated into separate generalized - functions to be reusable by multiple test groups. - """ - # get needed fields from Antarctica dataset - f = netCDF4.Dataset('antarctica_8km_2020_10_20.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - config = self.config - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Thwaites Glacier mesh. - xx0 = -1864434 - xx1 = -975432 - yy0 = -901349 - yy1 = 0 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - flood_mask = gridded_flood_fill(thk) - thk[flood_mask == 0] = 0.0 - vx[flood_mask == 0] = 0.0 - vy[flood_mask == 0] = 0.0 - - # Calculate distances to ice edge and grounding line - dist_to_edge, dist_to_GL = get_dist_to_edge_and_GL(self, thk, - topg, x1, y1, - window_size=1.e5) - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, section='high_res_mesh', thk=thk, - vx=vx, vy=vy, dist_to_edge=dist_to_edge, - dist_to_grounding_line=dist_to_GL) - - return cell_width.astype('float64'), x1.astype('float64'), \ - y1.astype('float64'), geom_points, geom_edges From 83f27c9df4735968f4ab5900a24424c600478d2f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 7 Apr 2023 13:08:18 -0700 Subject: [PATCH 02/19] Add parameters and returns to doc-string for build_cell_width() --- compass/landice/mesh.py | 38 +++++++++++++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 5 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index c254900d4b..eda514ffaf 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -419,6 +419,34 @@ def build_cell_width(self, section_name, gridded_dataset, """ Determine MPAS mesh cell size based on user-defined density function. + Parameters + ---------- + section : str + section of config file used to define mesh parameters + gridded_dataset : str + name of .nc file used to define cell spacing + flood_fill_start : list of ints + i and j indices used to define starting location for flood fill. + Most cases will use [None, None], which will just start the flood + fill in the center of the gridded dataset. + + Returns + ------- + cell_width : numpy.ndarray + Desired width of MPAS cells based on mesh desnity functions to pass to + :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + x1 : float + x coordinates from gridded dataset + y1 : float + y coordinates from gridded dataset + geom_points : jigsawpy.jigsaw_msh_t.VERT2_t + xy node coordinates to pass to build_planar_mesh() + geom_edges : jigsawpy.jigsaw_msh_t.EDGE2_t + xy edge coordinates between nodes to pass to build_planar_mesh() + flood_mask : numpy.ndarray + mask calculated by the flood fill routine, + where cells connected to the ice sheet (or main feature) + are 1 and everything else is 0. """ section = self.config[section_name] @@ -456,10 +484,10 @@ def build_cell_width(self, section_name, gridded_dataset, xx0, xx1, yy0, yy1) # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 + flood_mask = gridded_flood_fill(thk) + thk[flood_mask == 0] = 0.0 + vx[flood_mask == 0] = 0.0 + vy[flood_mask == 0] = 0.0 # Calculate distance from each grid point to ice edge # and grounding line, for use in cell spacing functions. @@ -476,4 +504,4 @@ def build_cell_width(self, section_name, gridded_dataset, flood_fill_jStart=flood_fill_start[1]) return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges, floodMask) + y1.astype('float64'), geom_points, geom_edges, flood_mask) From 73e37e7cbb95999987f6e87d3107ed249cc015f3 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 7 Apr 2023 16:16:53 -0700 Subject: [PATCH 03/19] Move most jigsaw and mpas-tools operations to compass.landice.mesh Move most jigsaw and mpas-tools operations to compass.landice.mesh.build_MALI_mesh() to greatly reduce redundant code between mesh generation test cases. --- compass/landice/mesh.py | 109 ++++++++++++++++++ compass/landice/tests/antarctica/mesh.py | 99 ++-------------- compass/landice/tests/greenland/mesh.py | 91 ++------------- compass/landice/tests/humboldt/mesh.py | 118 ++----------------- compass/landice/tests/kangerlussuaq/mesh.py | 113 ++----------------- compass/landice/tests/koge_bugt_s/mesh.py | 113 ++----------------- compass/landice/tests/thwaites/mesh.py | 119 ++------------------ 7 files changed, 171 insertions(+), 591 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index eda514ffaf..43588c5a5e 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -2,6 +2,11 @@ import jigsawpy import numpy as np +import xarray +from mpas_tools.io import write_netcdf +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 netCDF4 import Dataset @@ -505,3 +510,107 @@ def build_cell_width(self, section_name, gridded_dataset, return (cell_width.astype('float64'), x1.astype('float64'), y1.astype('float64'), geom_points, geom_edges, flood_mask) + + +def build_MALI_mesh(self, cell_width, x1, y1, geom_points, + geom_edges, mesh_name, section_name, + gridded_dataset, projection, geojson_file=None): + + logger = self.logger + section = self.config[section_name] + + logger.info('calling build_planar_mesh') + build_planar_mesh(cell_width, x1, y1, geom_points, + geom_edges, logger=logger) + dsMesh = xarray.open_dataset('base_mesh.nc') + logger.info('culling mesh') + dsMesh = cull(dsMesh, logger=logger) + logger.info('converting to MPAS mesh') + dsMesh = convert(dsMesh, logger=logger) + logger.info('writing grid_converted.nc') + write_netcdf(dsMesh, 'grid_converted.nc') + levels = section.get('levels') + logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') + args = ['create_landice_grid_from_generic_MPAS_grid.py', + '-i', 'grid_converted.nc', + '-o', 'grid_preCull.nc', + '-l', levels, '-v', 'glimmer'] + check_call(args, logger=logger) + + logger.info('calling interpolate_to_mpasli_grid.py') + args = ['interpolate_to_mpasli_grid.py', '-s', + gridded_dataset, '-d', + 'grid_preCull.nc', '-m', 'b', '-t'] + + check_call(args, logger=logger) + + cullDistance = section.get('cull_distance') + if float(cullDistance) > 0.: + logger.info('calling define_cullMask.py') + args = ['define_cullMask.py', '-f', + 'grid_preCull.nc', '-m', + 'distance', '-d', cullDistance] + + check_call(args, logger=logger) + else: + logger.info('cullDistance <= 0 in config file. ' + 'Will not cull by distance to margin. \n') + + if geojson_file is not None: + # This step is only necessary because the GeoJSON region + # is defined by lat-lon. + logger.info('calling set_lat_lon_fields_in_planar_grid.py') + args = ['set_lat_lon_fields_in_planar_grid.py', '-f', + 'grid_preCull.nc', '-p', projection] + + check_call(args, logger=logger) + + logger.info('calling MpasMaskCreator.x') + args = ['MpasMaskCreator.x', 'grid_preCull.nc', + 'mask.nc', '-f', geojson_file] + + check_call(args, logger=logger) + + logger.info('culling to geojson file') + + dsMesh = xarray.open_dataset('grid_preCull.nc') + if geojson_file is not None: + mask = xarray.open_dataset('mask.nc') + else: + mask = None + + dsMesh = cull(dsMesh, dsInverse=mask, logger=logger) + write_netcdf(dsMesh, 'culled.nc') + + logger.info('Marking horns for culling') + args = ['mark_horns_for_culling.py', '-f', 'culled.nc'] + check_call(args, logger=logger) + + logger.info('culling and converting') + dsMesh = xarray.open_dataset('culled.nc') + dsMesh = cull(dsMesh, logger=logger) + dsMesh = convert(dsMesh, logger=logger) + write_netcdf(dsMesh, 'dehorned.nc') + + logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') + args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', + 'dehorned.nc', '-o', + mesh_name, '-l', levels, '-v', 'glimmer', + '--beta', '--thermal', '--obs', '--diri'] + + check_call(args, logger=logger) + + logger.info('calling interpolate_to_mpasli_grid.py') + args = ['interpolate_to_mpasli_grid.py', '-s', + gridded_dataset, '-d', mesh_name, '-m', 'b'] + check_call(args, logger=logger) + + logger.info('Marking domain boundaries dirichlet') + args = ['mark_domain_boundaries_dirichlet.py', + '-f', mesh_name] + 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', + mesh_name, '-p', projection] + check_call(args, logger=logger) diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index d1d636ef45..7e8834523d 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/compass/landice/tests/antarctica/mesh.py @@ -3,14 +3,10 @@ import matplotlib.pyplot as plt import mpas_tools import netCDF4 -import xarray from geometric_features import FeatureCollection, GeometricFeatures -from mpas_tools.io import write_netcdf 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 compass.landice.mesh import build_cell_width +from compass.landice.mesh import build_cell_width, build_MALI_mesh from compass.model import make_graph_file from compass.step import Step @@ -57,31 +53,13 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['antarctica'] + section_name = 'antarctica' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodFillMask = \ build_cell_width( self, section_name='antarctica', gridded_dataset='antarctica_8km_2020_10_20.nc') - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'ais_8km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) # Apply floodFillMask to thickness field to help with culling copyfile('antarctica_8km_2020_10_20.nc', @@ -93,79 +71,26 @@ def run(self): gg.variables['vy'][0, :, :] *= floodFillMask gg.close() - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'antarctica_8km_2020_10_20_floodFillMask.nc', '-d', - 'ais_8km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # Cull a certain distance from the ice margin - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'ais_8km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - dsMesh = xarray.open_dataset('ais_8km_preCull.nc') - dsMesh = cull(dsMesh, logger=logger) - write_netcdf(dsMesh, 'antarctica_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'antarctica_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('antarctica_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'antarctica_dehorned.nc') - - mesh_filename = self.mesh_filename - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'antarctica_dehorned.nc', '-o', - mesh_filename, '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'antarctica_8km_2020_10_20.nc', - '-d', mesh_filename, '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', mesh_filename] - 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', - mesh_filename, '-p', 'ais-bedmap2'] - check_call(args, logger=logger) + build_MALI_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=self.mesh_filename, section_name=section_name, + gridded_dataset='antarctica_8km_2020_10_20.nc', + projection='ais-bedmap2', geojson_file=None) logger.info('creating graph.info') - make_graph_file(mesh_filename=mesh_filename, + make_graph_file(mesh_filename=self.mesh_filename, graph_filename='graph.info') # create a region mask - mask_filename = f'{mesh_filename[:-3]}_imbie_regionMasks.nc' - self._make_region_masks(mesh_filename, mask_filename, + mask_filename = f'{self.mesh_filename[:-3]}_imbie_regionMasks.nc' + self._make_region_masks(self.mesh_filename, mask_filename, self.cpus_per_task, tags=['EastAntarcticaIMBIE', 'WestAntarcticaIMBIE', 'AntarcticPeninsulaIMBIE']) - mask_filename = f'{mesh_filename[:-3]}_ismip6_regionMasks.nc' - self._make_region_masks(mesh_filename, mask_filename, + mask_filename = f'{self.mesh_filename[:-3]}_ismip6_regionMasks.nc' + self._make_region_masks(self.mesh_filename, mask_filename, self.cpus_per_task, tags=['ISMIP6_Basin']) diff --git a/compass/landice/tests/greenland/mesh.py b/compass/landice/tests/greenland/mesh.py index 2f8559f734..c7228bbc4f 100644 --- a/compass/landice/tests/greenland/mesh.py +++ b/compass/landice/tests/greenland/mesh.py @@ -1,10 +1,4 @@ -import xarray -from mpas_tools.io import write_netcdf -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 compass.landice.mesh import build_cell_width +from compass.landice.mesh import build_cell_width, build_MALI_mesh from compass.model import make_graph_file from compass.step import Step @@ -49,87 +43,22 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['high_res_GIS_mesh'] + mesh_name = 'GIS.nc' + section_name = 'high_res_GIS_mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ build_cell_width( - self, section_name='high_res_GIS_mesh', + self, section_name=section_name, gridded_dataset='greenland_2km_2020_04_20.epsg3413.nc', flood_fill_start=[100, 700]) - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'gis_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', '-d', - 'gis_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - cullDistance = section.get('cull_distance') - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'gis_1km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - - dsMesh = xarray.open_dataset('gis_1km_preCull.nc') - dsMesh = cull(dsMesh, logger=logger) - write_netcdf(dsMesh, 'greenland_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'greenland_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('greenland_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'greenland_dehorned.nc') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'greenland_dehorned.nc', '-o', - 'GIS.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - '-d', 'GIS.nc', '-m', 'b'] - 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) - 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) + build_MALI_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa + projection='gis-gimp', geojson_file=None) logger.info('creating graph.info') - make_graph_file(mesh_filename='GIS.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') diff --git a/compass/landice/tests/humboldt/mesh.py b/compass/landice/tests/humboldt/mesh.py index 5ee18a8c46..e14e850818 100644 --- a/compass/landice/tests/humboldt/mesh.py +++ b/compass/landice/tests/humboldt/mesh.py @@ -1,10 +1,4 @@ -import xarray -from mpas_tools.io import write_netcdf -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 compass.landice.mesh import build_cell_width +from compass.landice.mesh import build_cell_width, build_MALI_mesh from compass.model import make_graph_file from compass.step import Step @@ -53,113 +47,21 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['humboldt_mesh'] + section_name = 'humboldt_mesh' + mesh_name = 'Humboldt.nc' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ build_cell_width( - self, section_name='humboldt_mesh', + self, section_name=section_name, gridded_dataset='greenland_2km_2020_04_20.epsg3413.nc') - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - # If no number of levels specified in config file, use 10 - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'gis_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - # This step uses a subset of the whole Greenland dataset trimmed to - # the region around Humboldt Glacier, to speed up interpolation. - # This could also be replaced with the full Greenland Ice Sheet - # dataset. - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', '-d', - 'gis_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # This step is only necessary if you wish to cull a certain - # distance from the ice margin, within the bounds defined by - # the GeoJSON file. - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'gis_1km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - # This step is only necessary because the GeoJSON region - # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'gis_1km_preCull.nc', '-p', 'gis-gimp'] - - check_call(args, logger=logger) - - logger.info('calling MpasMaskCreator.x') - args = ['MpasMaskCreator.x', 'gis_1km_preCull.nc', - 'humboldt_mask.nc', '-f', 'Humboldt.geojson'] - - check_call(args, logger=logger) - - logger.info('culling to geojson file') - dsMesh = xarray.open_dataset('gis_1km_preCull.nc') - humboldtMask = xarray.open_dataset('humboldt_mask.nc') - dsMesh = cull(dsMesh, dsInverse=humboldtMask, logger=logger) - write_netcdf(dsMesh, 'humboldt_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'humboldt_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('humboldt_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'humboldt_dehorned.nc') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'humboldt_dehorned.nc', '-o', - 'Humboldt_1to10km.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', - '-d', 'Humboldt_1to10km.nc', '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'Humboldt_1to10km.nc'] - 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', - 'Humboldt_1to10km.nc', '-p', 'gis-gimp'] - check_call(args, logger=logger) + build_MALI_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', + projection='gis-gimp', geojson_file='Humboldt.geojson') logger.info('creating graph.info') - make_graph_file(mesh_filename='Humboldt_1to10km.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') diff --git a/compass/landice/tests/kangerlussuaq/mesh.py b/compass/landice/tests/kangerlussuaq/mesh.py index b126fcc63f..c319caa6ce 100644 --- a/compass/landice/tests/kangerlussuaq/mesh.py +++ b/compass/landice/tests/kangerlussuaq/mesh.py @@ -1,10 +1,4 @@ -import xarray -from mpas_tools.io import write_netcdf -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 compass.landice.mesh import build_cell_width +from compass.landice.mesh import build_cell_width, build_MALI_mesh from compass.model import make_graph_file from compass.step import Step @@ -54,108 +48,21 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['high_res_Kangerlussuaq_mesh'] + mesh_name = 'Kangerlussuaq.nc' + section_name = 'high_res_Kangerlussuaq_mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ build_cell_width( - self, section_name='high_res_Kangerlussuaq_mesh', + self, section_name=section_name, gridded_dataset='greenland_8km_2020_04_20.epsg3413.nc') - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'gis_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', '-d', - 'gis_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # This step is only necessary if you wish to cull a certain - # distance from the ice margin, within the bounds defined by - # the GeoJSON file. - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'gis_1km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - # This step is only necessary because the GeoJSON region - # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'gis_1km_preCull.nc', '-p', 'gis-gimp'] - - check_call(args, logger=logger) - - logger.info('calling MpasMaskCreator.x') - args = ['MpasMaskCreator.x', 'gis_1km_preCull.nc', - 'kangerlussuaq_mask.nc', '-f', 'Kangerlussuaq.geojson'] - - check_call(args, logger=logger) - - logger.info('culling to geojson file') - dsMesh = xarray.open_dataset('gis_1km_preCull.nc') - kangerMask = xarray.open_dataset('kangerlussuaq_mask.nc') - dsMesh = cull(dsMesh, dsInverse=kangerMask, logger=logger) - write_netcdf(dsMesh, 'kangerlussuaq_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'kangerlussuaq_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('kangerlussuaq_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'kangerlussuaq_dehorned.nc') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'kangerlussuaq_dehorned.nc', '-o', - 'Kangerlussuaq.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - '-d', 'Kangerlussuaq.nc', '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'Kangerlussuaq.nc'] - 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', - 'Kangerlussuaq.nc', '-p', 'gis-gimp'] - check_call(args, logger=logger) + build_MALI_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa + projection='gis-gimp', geojson_file='Kangerlussuaq.geojson') logger.info('creating graph.info') - make_graph_file(mesh_filename='Kangerlussuaq.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') diff --git a/compass/landice/tests/koge_bugt_s/mesh.py b/compass/landice/tests/koge_bugt_s/mesh.py index 775762508c..dc3e1072db 100644 --- a/compass/landice/tests/koge_bugt_s/mesh.py +++ b/compass/landice/tests/koge_bugt_s/mesh.py @@ -1,10 +1,4 @@ -import xarray -from mpas_tools.io import write_netcdf -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 compass.landice.mesh import build_cell_width +from compass.landice.mesh import build_cell_width, build_MALI_mesh from compass.model import make_graph_file from compass.step import Step @@ -54,108 +48,21 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['high_res_KogeBugtS_mesh'] + mesh_name = 'Koge_Bugt_S.nc' + section_name = 'high_res_KogeBugtS_mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ build_cell_width( - self, section_name='high_res_KogeBugtS_mesh', + self, section_name=section_name, gridded_dataset='greenland_8km_2020_04_20.epsg3413.nc') - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'gis_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', '-d', - 'gis_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # This step is only necessary if you wish to cull a certain - # distance from the ice margin, within the bounds defined by - # the GeoJSON file. - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'gis_1km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - # This step is only necessary because the GeoJSON region - # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'gis_1km_preCull.nc', '-p', 'gis-gimp'] - - check_call(args, logger=logger) - - logger.info('calling MpasMaskCreator.x') - args = ['MpasMaskCreator.x', 'gis_1km_preCull.nc', - 'koge_bugt_s_mask.nc', '-f', 'Koge_Bugt_S.geojson'] - - check_call(args, logger=logger) - - logger.info('culling to geojson file') - dsMesh = xarray.open_dataset('gis_1km_preCull.nc') - kangerMask = xarray.open_dataset('koge_bugt_s_mask.nc') - dsMesh = cull(dsMesh, dsInverse=kangerMask, logger=logger) - write_netcdf(dsMesh, 'koge_bugt_s_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'koge_bugt_s_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('koge_bugt_s_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'koge_bugt_s_dehorned.nc') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'koge_bugt_s_dehorned.nc', '-o', - 'Koge_Bugt_S.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - '-d', 'Koge_Bugt_S.nc', '-m', 'b', '-t'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'Koge_Bugt_S.nc'] - 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', - 'Koge_Bugt_S.nc', '-p', 'gis-gimp'] - check_call(args, logger=logger) + build_MALI_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa + projection='gis-gimp', geojson_file='Koge_Bugt_S.geojson') logger.info('creating graph.info') - make_graph_file(mesh_filename='Koge_Bugt_S.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') diff --git a/compass/landice/tests/thwaites/mesh.py b/compass/landice/tests/thwaites/mesh.py index 5a35383151..695762c18b 100644 --- a/compass/landice/tests/thwaites/mesh.py +++ b/compass/landice/tests/thwaites/mesh.py @@ -1,10 +1,4 @@ -import xarray -from mpas_tools.io import write_netcdf -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 compass.landice.mesh import build_cell_width +from compass.landice.mesh import build_cell_width, build_MALI_mesh from compass.model import make_graph_file from compass.step import Step @@ -43,114 +37,21 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['high_res_mesh'] + mesh_name = 'Thwaites.nc' + section_name = 'high_res_mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ build_cell_width( - self, section_name='high_res_mesh', + self, section_name=section_name, gridded_dataset='antarctica_8km_2020_10_20.nc') - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - ds_mesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - ds_mesh = cull(ds_mesh, logger=logger) - logger.info('converting to MPAS mesh') - ds_mesh = convert(ds_mesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(ds_mesh, 'grid_converted.nc') - - levels = section.get('levels') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'ase_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - # This step uses a subset of the whole Antarctica dataset trimmed to - # the Amundsen Sean Embayment, to speed up interpolation. - # This could also be replaced with the full Antarctic Ice Sheet - # dataset. - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'antarctica_1km_2020_10_20_ASE.nc', '-d', - 'ase_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # This step is only necessary if you wish to cull a certain - # distance from the ice margin, within the bounds defined by - # the GeoJSON file. - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'ase_1km_preCull.nc', '-m' - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - # This step is only necessary because the GeoJSON region - # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'ase_1km_preCull.nc', '-p', 'ais-bedmap2'] - - check_call(args, logger=logger) - - logger.info('calling MpasMaskCreator.x') - args = ['MpasMaskCreator.x', 'ase_1km_preCull.nc', - 'thwaites_mask.nc', '-f', 'thwaites_minimal.geojson'] - - check_call(args, logger=logger) - - logger.info('culling to geojson file') - ds_mesh = xarray.open_dataset('ase_1km_preCull.nc') - thwaitesMask = xarray.open_dataset('thwaites_mask.nc') - ds_mesh = cull(ds_mesh, dsInverse=thwaitesMask, logger=logger) - write_netcdf(ds_mesh, 'thwaites_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'thwaites_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - ds_mesh = xarray.open_dataset('thwaites_culled.nc') - ds_mesh = cull(ds_mesh, logger=logger) - ds_mesh = convert(ds_mesh, logger=logger) - write_netcdf(ds_mesh, 'thwaites_dehorned.nc') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'thwaites_dehorned.nc', '-o', - 'Thwaites.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'antarctica_1km_2020_10_20_ASE.nc', - '-d', 'Thwaites.nc', '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'Thwaites.nc'] - 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', - 'Thwaites.nc', '-p', 'ais-bedmap2'] - check_call(args, logger=logger) + build_MALI_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='antarctica_1km_2020_10_20_ASE.nc', + projection='ais-bedmap2', geojson_file='thwaites_minimal.geojson') logger.info('creating graph.info') - make_graph_file(mesh_filename='Thwaites.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') From cd612c4a0078a282f2fa5705bd0f6e2d22ba22e5 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 10 Apr 2023 09:57:14 -0600 Subject: [PATCH 04/19] Standardize mesh generation test case and config section names Standardize mesh generation test case and config section names, using mesh_gen for the test case and 'mesh' for the config section. All .cfg files pertaining only to mesh generation have been moved into the mesh_gen modules. --- compass/landice/tests/antarctica/mesh.py | 4 ++-- .../{antarctica.cfg => mesh_gen/mesh_gen.cfg} | 2 +- compass/landice/tests/greenland/__init__.py | 12 +++++++----- compass/landice/tests/greenland/mesh.py | 2 +- .../{high_res_mesh => mesh_gen}/__init__.py | 6 +++--- .../high_res_mesh.cfg => mesh_gen/mesh_gen.cfg} | 2 +- compass/landice/tests/humboldt/mesh.py | 2 +- .../humboldt/{humboldt.cfg => mesh_gen/mesh_gen.cfg} | 2 +- compass/landice/tests/kangerlussuaq/mesh.py | 2 +- .../tests/kangerlussuaq/mesh_gen/mesh_gen.cfg | 2 +- compass/landice/tests/koge_bugt_s/mesh.py | 2 +- .../landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg | 2 +- compass/landice/tests/thwaites/__init__.py | 6 +++--- compass/landice/tests/thwaites/mesh.py | 2 +- .../thwaites/{high_res_mesh => mesh_gen}/__init__.py | 6 +++--- .../high_res_mesh.cfg => mesh_gen/mesh_gen.cfg} | 2 +- .../landice/test_groups/greenland.rst | 6 +++--- .../landice/test_groups/thwaites.rst | 6 +++--- docs/users_guide/landice/test_groups/greenland.rst | 12 ++++++------ docs/users_guide/landice/test_groups/humboldt.rst | 2 +- .../landice/test_groups/kangerlussuaq.rst | 2 +- docs/users_guide/landice/test_groups/koge_bugt_s.rst | 2 +- docs/users_guide/landice/test_groups/thwaites.rst | 10 +++++----- 23 files changed, 49 insertions(+), 47 deletions(-) rename compass/landice/tests/antarctica/{antarctica.cfg => mesh_gen/mesh_gen.cfg} (98%) rename compass/landice/tests/greenland/{high_res_mesh => mesh_gen}/__init__.py (92%) rename compass/landice/tests/greenland/{high_res_mesh/high_res_mesh.cfg => mesh_gen/mesh_gen.cfg} (98%) rename compass/landice/tests/humboldt/{humboldt.cfg => mesh_gen/mesh_gen.cfg} (98%) rename compass/landice/tests/thwaites/{high_res_mesh => mesh_gen}/__init__.py (92%) rename compass/landice/tests/thwaites/{high_res_mesh/high_res_mesh.cfg => mesh_gen/mesh_gen.cfg} (98%) diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index 7e8834523d..8f465db0f6 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/compass/landice/tests/antarctica/mesh.py @@ -53,12 +53,12 @@ def run(self): Run this step of the test case """ logger = self.logger - section_name = 'antarctica' + section_name = 'mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodFillMask = \ build_cell_width( - self, section_name='antarctica', + self, section_name=section_name, gridded_dataset='antarctica_8km_2020_10_20.nc') # Apply floodFillMask to thickness field to help with culling diff --git a/compass/landice/tests/antarctica/antarctica.cfg b/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg similarity index 98% rename from compass/landice/tests/antarctica/antarctica.cfg rename to compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg index f3d7cb7d36..917f706f52 100644 --- a/compass/landice/tests/antarctica/antarctica.cfg +++ b/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg @@ -1,5 +1,5 @@ # config options for antarctica test cases -[antarctica] +[mesh] # number of levels in the mesh levels = 5 diff --git a/compass/landice/tests/greenland/__init__.py b/compass/landice/tests/greenland/__init__.py index 881438431d..500338d5ff 100644 --- a/compass/landice/tests/greenland/__init__.py +++ b/compass/landice/tests/greenland/__init__.py @@ -1,8 +1,10 @@ -from compass.testgroup import TestGroup -from compass.landice.tests.greenland.smoke_test import SmokeTest -from compass.landice.tests.greenland.decomposition_test import DecompositionTest +from compass.landice.tests.greenland.decomposition_test import ( + DecompositionTest, +) +from compass.landice.tests.greenland.mesh_gen import MeshGen from compass.landice.tests.greenland.restart_test import RestartTest -from compass.landice.tests.greenland.high_res_mesh import HighResMesh +from compass.landice.tests.greenland.smoke_test import SmokeTest +from compass.testgroup import TestGroup class Greenland(TestGroup): @@ -27,4 +29,4 @@ def __init__(self, mpas_core): RestartTest(test_group=self, velo_solver=velo_solver)) self.add_test_case( - HighResMesh(test_group=self)) + MeshGen(test_group=self)) diff --git a/compass/landice/tests/greenland/mesh.py b/compass/landice/tests/greenland/mesh.py index c7228bbc4f..746d70d654 100644 --- a/compass/landice/tests/greenland/mesh.py +++ b/compass/landice/tests/greenland/mesh.py @@ -44,7 +44,7 @@ def run(self): """ logger = self.logger mesh_name = 'GIS.nc' - section_name = 'high_res_GIS_mesh' + section_name = 'mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ diff --git a/compass/landice/tests/greenland/high_res_mesh/__init__.py b/compass/landice/tests/greenland/mesh_gen/__init__.py similarity index 92% rename from compass/landice/tests/greenland/high_res_mesh/__init__.py rename to compass/landice/tests/greenland/mesh_gen/__init__.py index 8b09f9a13d..b8608994e6 100644 --- a/compass/landice/tests/greenland/high_res_mesh/__init__.py +++ b/compass/landice/tests/greenland/mesh_gen/__init__.py @@ -1,8 +1,8 @@ -from compass.testcase import TestCase from compass.landice.tests.greenland.mesh import Mesh +from compass.testcase import TestCase -class HighResMesh(TestCase): +class MeshGen(TestCase): """ The high resolution test case for the greenland test group simply creates the mesh and initial condition. @@ -18,7 +18,7 @@ def __init__(self, test_group): test_group : compass.landice.tests.greenland.Greenland The test group that this test case belongs to """ - name = 'high_res_mesh' + name = 'mesh_gen' subdir = name super().__init__(test_group=test_group, name=name, subdir=subdir) diff --git a/compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg b/compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg similarity index 98% rename from compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg rename to compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg index e77869bd89..cea3358b46 100644 --- a/compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg +++ b/compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg @@ -1,5 +1,5 @@ # config options for high_res_mesh test case -[high_res_GIS_mesh] +[mesh] # number of levels in the mesh levels = 10 diff --git a/compass/landice/tests/humboldt/mesh.py b/compass/landice/tests/humboldt/mesh.py index e14e850818..2f7296c432 100644 --- a/compass/landice/tests/humboldt/mesh.py +++ b/compass/landice/tests/humboldt/mesh.py @@ -47,7 +47,7 @@ def run(self): Run this step of the test case """ logger = self.logger - section_name = 'humboldt_mesh' + section_name = 'mesh' mesh_name = 'Humboldt.nc' logger.info('calling build_cell_width') diff --git a/compass/landice/tests/humboldt/humboldt.cfg b/compass/landice/tests/humboldt/mesh_gen/mesh_gen.cfg similarity index 98% rename from compass/landice/tests/humboldt/humboldt.cfg rename to compass/landice/tests/humboldt/mesh_gen/mesh_gen.cfg index 59db5bc1b6..6927514134 100644 --- a/compass/landice/tests/humboldt/humboldt.cfg +++ b/compass/landice/tests/humboldt/mesh_gen/mesh_gen.cfg @@ -1,5 +1,5 @@ # config options for humboldt test cases -[humboldt_mesh] +[mesh] # number of levels in the mesh levels = 10 diff --git a/compass/landice/tests/kangerlussuaq/mesh.py b/compass/landice/tests/kangerlussuaq/mesh.py index c319caa6ce..cab345cfc2 100644 --- a/compass/landice/tests/kangerlussuaq/mesh.py +++ b/compass/landice/tests/kangerlussuaq/mesh.py @@ -49,7 +49,7 @@ def run(self): """ logger = self.logger mesh_name = 'Kangerlussuaq.nc' - section_name = 'high_res_Kangerlussuaq_mesh' + section_name = 'mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ diff --git a/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg b/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg index 355c4d01d5..72c5a9da2f 100644 --- a/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg @@ -1,5 +1,5 @@ # config options for kangerlussuaq test cases -[high_res_Kangerlussuaq_mesh] +[mesh] # number of levels in the mesh levels = 10 diff --git a/compass/landice/tests/koge_bugt_s/mesh.py b/compass/landice/tests/koge_bugt_s/mesh.py index dc3e1072db..241c829996 100644 --- a/compass/landice/tests/koge_bugt_s/mesh.py +++ b/compass/landice/tests/koge_bugt_s/mesh.py @@ -49,7 +49,7 @@ def run(self): """ logger = self.logger mesh_name = 'Koge_Bugt_S.nc' - section_name = 'high_res_KogeBugtS_mesh' + section_name = 'mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ diff --git a/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg b/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg index 1b38c73804..770432a6d8 100644 --- a/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg @@ -1,5 +1,5 @@ # config options for koge_bugt_s test cases -[high_res_KogeBugtS_mesh] +[mesh] # number of levels in the mesh levels = 10 diff --git a/compass/landice/tests/thwaites/__init__.py b/compass/landice/tests/thwaites/__init__.py index 74343beda3..27a764eaa9 100644 --- a/compass/landice/tests/thwaites/__init__.py +++ b/compass/landice/tests/thwaites/__init__.py @@ -1,7 +1,7 @@ -from compass.testgroup import TestGroup from compass.landice.tests.thwaites.decomposition_test import DecompositionTest +from compass.landice.tests.thwaites.mesh_gen import MeshGen from compass.landice.tests.thwaites.restart_test import RestartTest -from compass.landice.tests.thwaites.high_res_mesh import HighResMesh +from compass.testgroup import TestGroup class Thwaites(TestGroup): @@ -20,4 +20,4 @@ def __init__(self, mpas_core): self.add_test_case(RestartTest(test_group=self)) - self.add_test_case(HighResMesh(test_group=self)) + self.add_test_case(MeshGen(test_group=self)) diff --git a/compass/landice/tests/thwaites/mesh.py b/compass/landice/tests/thwaites/mesh.py index 695762c18b..ce9447376b 100644 --- a/compass/landice/tests/thwaites/mesh.py +++ b/compass/landice/tests/thwaites/mesh.py @@ -38,7 +38,7 @@ def run(self): """ logger = self.logger mesh_name = 'Thwaites.nc' - section_name = 'high_res_mesh' + section_name = 'mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodMask = \ diff --git a/compass/landice/tests/thwaites/high_res_mesh/__init__.py b/compass/landice/tests/thwaites/mesh_gen/__init__.py similarity index 92% rename from compass/landice/tests/thwaites/high_res_mesh/__init__.py rename to compass/landice/tests/thwaites/mesh_gen/__init__.py index 8a44ee7f62..184a64642b 100644 --- a/compass/landice/tests/thwaites/high_res_mesh/__init__.py +++ b/compass/landice/tests/thwaites/mesh_gen/__init__.py @@ -1,8 +1,8 @@ -from compass.testcase import TestCase from compass.landice.tests.thwaites.mesh import Mesh +from compass.testcase import TestCase -class HighResMesh(TestCase): +class MeshGen(TestCase): """ The high resolution test case for the thwaites test group simply creates the mesh and initial condition. @@ -18,7 +18,7 @@ def __init__(self, test_group): test_group : compass.landice.tests.thwaites.Thwaites The test group that this test case belongs to """ - name = 'high_res_mesh' + name = 'mesh_gen' subdir = name super().__init__(test_group=test_group, name=name, subdir=subdir) diff --git a/compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg b/compass/landice/tests/thwaites/mesh_gen/mesh_gen.cfg similarity index 98% rename from compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg rename to compass/landice/tests/thwaites/mesh_gen/mesh_gen.cfg index b5336c9753..98ae96049e 100644 --- a/compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg +++ b/compass/landice/tests/thwaites/mesh_gen/mesh_gen.cfg @@ -1,5 +1,5 @@ # config options for high_res_mesh test case -[high_res_mesh] +[mesh] # number of levels in the mesh levels = 10 diff --git a/docs/developers_guide/landice/test_groups/greenland.rst b/docs/developers_guide/landice/test_groups/greenland.rst index ca464c4d86..5b27d62983 100644 --- a/docs/developers_guide/landice/test_groups/greenland.rst +++ b/docs/developers_guide/landice/test_groups/greenland.rst @@ -44,7 +44,7 @@ mesh The class :py:class:`compass.landice.tests.greenland.mesh.Mesh` defines a step for creating a variable resolution Greenland Ice Sheet mesh. -This is used by the ``high_res_mesh`` test case. +This is used by the ``mesh_gen`` test case. .. _dev_landice_greenland_smoke_test: @@ -83,9 +83,9 @@ the latter perform a 2-day restart run beginning with the end of the first. .. _dev_landice_greenland_high_res_mesh: -high_res_mesh +mesh_gen ------------- -The :py:class:`compass.landice.tests.greenland.high_res_mesh.HighResMesh` +The :py:class:`compass.landice.tests.greenland.mesh_gen.MeshGen` calls the :py:class:`compass.landice.tests.greenland.mesh.Mesh` to create the variable resolution Greenland Ice Sheet mesh. diff --git a/docs/developers_guide/landice/test_groups/thwaites.rst b/docs/developers_guide/landice/test_groups/thwaites.rst index 763e0eabc1..7217d948d4 100644 --- a/docs/developers_guide/landice/test_groups/thwaites.rst +++ b/docs/developers_guide/landice/test_groups/thwaites.rst @@ -68,11 +68,11 @@ mesh The class :py:class:`compass.landice.tests.thwaites.mesh.Mesh` defines a step for creating a variable resolution Thwaites Glacier mesh. -This is used by the ``high_res_mesh`` test case. +This is used by the ``mesh_gen`` test case. -high_res_mesh +mesh_gen ------------- -The :py:class:`compass.landice.tests.thwaites.high_res_mesh.HighResMesh` +The :py:class:`compass.landice.tests.thwaites.mesh_gen.MeshGen` calls the :py:class:`compass.landice.tests.thwaites.mesh.Mesh` to create the variable resolution Thwaites Glacier mesh. diff --git a/docs/users_guide/landice/test_groups/greenland.rst b/docs/users_guide/landice/test_groups/greenland.rst index 7db3806fad..67188d158c 100644 --- a/docs/users_guide/landice/test_groups/greenland.rst +++ b/docs/users_guide/landice/test_groups/greenland.rst @@ -15,24 +15,24 @@ The ``landice/greenland`` test group runs tests with a coarse (20-km) The test group includes 3 test cases, each of which has one or more steps that are variants on ``run_model`` (given other names in the decomposition and restart test cases to distinguish multiple model runs), which performs time -integration of the model. There is a fourth test case, ``high_res_mesh``, that +integration of the model. There is a fourth test case, ``mesh_gen``, that creates a variable resolution Greenland Ice Sheet mesh, with the step ``mesh``. The test cases in this test group can run with either the SIA or FO velocity solvers. Running with the FO solver requires a build of MALI that includes Albany, but the SIA variant of the test can be run without Albany. The FO version uses no-slip basal boundary condition. There is no integration step -for the test case ``high_res_mesh``. +for the test case ``mesh_gen``. config options -------------- -The ``high_res_mesh`` test case uses the default config options below. +The ``mesh_gen`` test case uses the default config options below. The other test cases do not use config options. .. code-block:: cfg - [high_res_mesh] + [mesh] # number of levels in the mesh levels = 10 @@ -85,10 +85,10 @@ from a restart file saved by the first. Prognostic variables are compared between the "full" and "restart" runs at year 2 to make sure they are bit-for-bit identical. -high_res_mesh +mesh_gen ------------- -``landice/greenland/high_res_mesh`` creates a variable resolution mesh based +``landice/greenland/mesh_gen`` creates a variable resolution mesh based on the the config options listed above. This will not be the same as the pre-generated 20 km mesh used in the other three test cases because it uses a newer version of Jigsaw. Note that the basal friction optimization is diff --git a/docs/users_guide/landice/test_groups/humboldt.rst b/docs/users_guide/landice/test_groups/humboldt.rst index ad11ba8d44..40aec04a7d 100644 --- a/docs/users_guide/landice/test_groups/humboldt.rst +++ b/docs/users_guide/landice/test_groups/humboldt.rst @@ -34,7 +34,7 @@ the mesh generation options are adjusted through the config file. .. code-block:: cfg - [humboldt_mesh] + [mesh] # number of levels in the mesh levels = 10 diff --git a/docs/users_guide/landice/test_groups/kangerlussuaq.rst b/docs/users_guide/landice/test_groups/kangerlussuaq.rst index f3b31d5973..5c91bd16a6 100644 --- a/docs/users_guide/landice/test_groups/kangerlussuaq.rst +++ b/docs/users_guide/landice/test_groups/kangerlussuaq.rst @@ -24,7 +24,7 @@ the mesh generation options are adjusted through the config file. .. code-block:: cfg - [high_res_Kangerlussuaq_mesh] + [mesh] # number of levels in the mesh levels = 10 diff --git a/docs/users_guide/landice/test_groups/koge_bugt_s.rst b/docs/users_guide/landice/test_groups/koge_bugt_s.rst index cb6a6a2db6..4f65ad5009 100644 --- a/docs/users_guide/landice/test_groups/koge_bugt_s.rst +++ b/docs/users_guide/landice/test_groups/koge_bugt_s.rst @@ -24,7 +24,7 @@ the mesh generation options are adjusted through the config file. .. code-block:: cfg - [high_res_KogeBugtS_mesh] + [mesh] # number of levels in the mesh levels = 10 diff --git a/docs/users_guide/landice/test_groups/thwaites.rst b/docs/users_guide/landice/test_groups/thwaites.rst index 79b44e87e6..9a7c4bf0b1 100644 --- a/docs/users_guide/landice/test_groups/thwaites.rst +++ b/docs/users_guide/landice/test_groups/thwaites.rst @@ -27,17 +27,17 @@ integration of the model. There is a not an explicit smoke test, but the The ``decomposition_test`` and ``restart_test`` test cases in this test group can only be run with the FO velocity solvers. Running with the FO solver requires a build of MALI that includes Albany. There is no integration step for the test -case ``high_res_mesh``. +case ``mesh_gen``. config options -------------- -The ``high_res_mesh`` test case uses the default config options below. +The ``mesh_gen`` test case uses the default config options below. The other test cases do not use config options. .. code-block:: cfg - [high_res_mesh] + [mesh] # number of levels in the mesh levels = 10 @@ -84,10 +84,10 @@ second begins from a restart file saved by the first. Prognostic variables are compared between the "full" and "restart" runs to make sure they are bit-for-bit identical. -high_res_mesh +mesh_gen ------------- -``landice/thwaites/high_res_mesh`` creates a variable resolution mesh based +``landice/thwaites/mesh_gen`` creates a variable resolution mesh based on the the config options listed above. This will not be the same as the pre-generated 4-14km mesh used in ``decomposition_test`` and ``restart_test`` because it uses a newer version of Jigsaw. Note that the basal friction From b69554d3b51c4889747704c8a8dd18bc39eb4587 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 10 Apr 2023 10:20:20 -0600 Subject: [PATCH 05/19] Update landice api in developers guide Update landice api in developers guide to include renamed methods and classes, and new functions in landice.mesh. --- docs/developers_guide/landice/api.rst | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index bddb48c64a..36e05c843d 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -205,8 +205,8 @@ greenland mesh.Mesh.run mesh.Mesh.build_cell_width - high_res_mesh.HighResMesh - high_res_mesh.HighResMesh.run + mesh_gen.MeshGen + mesh_gen.MeshGen.run humboldt ~~~~~~~~ @@ -400,6 +400,12 @@ thwaites run_model.RunModel.setup run_model.RunModel.run + mesh_gen.MeshGen + mesh_gen.MeshGen.run + + mesh.Mesh + mesh.Mesh.run + Landice Framework ^^^^^^^^^^^^^^^^^ @@ -414,3 +420,5 @@ Landice Framework mesh.set_rectangular_geom_points_and_edges mesh.set_cell_width mesh.get_dist_to_edge_and_GL + mesh.build_cell_width + mesh.build_MALI_mesh From 4532644a8569cf1177cc30fb741a869b3165dcb2 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 10 Apr 2023 10:23:01 -0600 Subject: [PATCH 06/19] Fix output file name in humboldt mesh_gen test case Fix inconsistent usage of Humboldt.nc and Humboldt_1to10km.nc that caused an error at the end of the test case. --- compass/landice/tests/humboldt/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/humboldt/mesh.py b/compass/landice/tests/humboldt/mesh.py index 2f7296c432..ca9f7b58ee 100644 --- a/compass/landice/tests/humboldt/mesh.py +++ b/compass/landice/tests/humboldt/mesh.py @@ -27,7 +27,7 @@ def __init__(self, test_case): super().__init__(test_case=test_case, name='mesh') self.add_output_file(filename='graph.info') - self.add_output_file(filename='Humboldt_1to10km.nc') + self.add_output_file(filename='Humboldt.nc') self.add_input_file( filename='humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', target='humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', From dbf1b87357732b65d6e6f4ce9029f390763bbdc7 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 10 Apr 2023 14:31:50 -0600 Subject: [PATCH 07/19] Add make_region_masks() function to landice.mesh. Add make_region_masks() function to landice.mesh, and add this capability for Greenland as well as Antarctica. --- compass/landice/mesh.py | 43 +++++++++++++++++- compass/landice/tests/antarctica/mesh.py | 56 ++++++------------------ compass/landice/tests/greenland/mesh.py | 19 +++++++- docs/developers_guide/landice/api.rst | 1 + 4 files changed, 74 insertions(+), 45 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 43588c5a5e..4947158cc1 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -3,7 +3,8 @@ import jigsawpy import numpy as np import xarray -from mpas_tools.io import write_netcdf +from geometric_features import FeatureCollection, GeometricFeatures +from mpas_tools.io import default_format, write_netcdf from mpas_tools.logging import check_call from mpas_tools.mesh.conversion import convert, cull from mpas_tools.mesh.creation import build_planar_mesh @@ -614,3 +615,43 @@ def build_MALI_mesh(self, cell_width, x1, y1, geom_points, args = ['set_lat_lon_fields_in_planar_grid.py', '-f', mesh_name, '-p', projection] check_call(args, logger=logger) + + +def make_region_masks(self, mesh_filename, mask_filename, cores, tags): + """ + Create masks for ice-sheet subregions based on data + in MPAS-Dev/geometric_fatures. + + Parameters + ---------- + mesh_filename : str + name of .nc mesh file for which to create region masks + mask_filename : str + name of .nc file to contain region masks + cores : int + number of processors used to create region masks + tags : list of str + Groups of regions for which masks are to be defined + """ + + logger = self.logger + logger.info('creating region masks') + gf = GeometricFeatures() + fcMask = FeatureCollection() + + for tag in tags: + fc = gf.read(componentName='landice', objectType='region', + tags=[tag]) + fcMask.merge(fc) + + geojson_filename = 'regionMask.geojson' + fcMask.to_geojson(geojson_filename) + + args = ['compute_mpas_region_masks', + '-m', mesh_filename, + '-g', geojson_filename, + '-o', mask_filename, + '-t', 'cell', + '--process_count', f'{cores}', + '--format', default_format] + check_call(args, logger=logger) diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index 8f465db0f6..dafe705e35 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/compass/landice/tests/antarctica/mesh.py @@ -1,12 +1,12 @@ from shutil import copyfile -import matplotlib.pyplot as plt -import mpas_tools import netCDF4 -from geometric_features import FeatureCollection, GeometricFeatures -from mpas_tools.logging import check_call -from compass.landice.mesh import build_cell_width, build_MALI_mesh +from compass.landice.mesh import ( + build_cell_width, + build_MALI_mesh, + make_region_masks, +) from compass.model import make_graph_file from compass.step import Step @@ -83,43 +83,13 @@ def run(self): # create a region mask mask_filename = f'{self.mesh_filename[:-3]}_imbie_regionMasks.nc' - self._make_region_masks(self.mesh_filename, mask_filename, - self.cpus_per_task, - tags=['EastAntarcticaIMBIE', - 'WestAntarcticaIMBIE', - 'AntarcticPeninsulaIMBIE']) + make_region_masks(self, self.mesh_filename, mask_filename, + self.cpus_per_task, + tags=['EastAntarcticaIMBIE', + 'WestAntarcticaIMBIE', + 'AntarcticPeninsulaIMBIE']) mask_filename = f'{self.mesh_filename[:-3]}_ismip6_regionMasks.nc' - self._make_region_masks(self.mesh_filename, mask_filename, - self.cpus_per_task, - tags=['ISMIP6_Basin']) - - def _make_region_masks(self, mesh_filename, mask_filename, cores, tags): - logger = self.logger - gf = GeometricFeatures() - fcMask = FeatureCollection() - - for tag in tags: - fc = gf.read(componentName='landice', objectType='region', - tags=[tag]) - fc.plot('southpole') - plt.savefig(f'plot_basins_{tag}.png') - fcMask.merge(fc) - - geojson_filename = 'regionMask.geojson' - fcMask.to_geojson(geojson_filename) - - # these defaults may have been updated from config options -- pass them - # along to the subprocess - netcdf_format = mpas_tools.io.default_format - netcdf_engine = mpas_tools.io.default_engine - - args = ['compute_mpas_region_masks', - '-m', mesh_filename, - '-g', geojson_filename, - '-o', mask_filename, - '-t', 'cell', - '--process_count', f'{cores}', - '--format', netcdf_format, - '--engine', netcdf_engine] - check_call(args, logger=logger) + make_region_masks(self, self.mesh_filename, mask_filename, + self.cpus_per_task, + tags=['ISMIP6_Basin']) diff --git a/compass/landice/tests/greenland/mesh.py b/compass/landice/tests/greenland/mesh.py index 746d70d654..40c3ddb828 100644 --- a/compass/landice/tests/greenland/mesh.py +++ b/compass/landice/tests/greenland/mesh.py @@ -1,4 +1,8 @@ -from compass.landice.mesh import build_cell_width, build_MALI_mesh +from compass.landice.mesh import ( + build_cell_width, + build_MALI_mesh, + make_region_masks, +) from compass.model import make_graph_file from compass.step import Step @@ -62,3 +66,16 @@ def run(self): logger.info('creating graph.info') make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') + + # create region masks + mask_filename = f'{mesh_name[:-3]}_regionMasks.nc' + make_region_masks(self, mesh_name, mask_filename, + self.cpus_per_task, + tags=['eastCentralGreenland', + 'northEastGreenland', + 'northGreenland', + 'northWestGreenland', + 'southEastGreenland', + 'southGreenland', + 'southWestGreenland', + 'westCentralGreenland']) diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 36e05c843d..d832d8a16e 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -422,3 +422,4 @@ Landice Framework mesh.get_dist_to_edge_and_GL mesh.build_cell_width mesh.build_MALI_mesh + mesh.make_region_masks From 14afb416118f67ca323ba470781a1b58d8e91c98 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 10 Apr 2023 17:40:40 -0600 Subject: [PATCH 08/19] Fix netcdf format and engine definitions Fix netcdf format and engine definitions, which for some reason need to be called as mpas_tools.io.default_engine rather than using: from mpas_tools.io import default_engine. The latter approach returned a None type variable, which threw an error when passed to compute_mpas_region_masks. --- compass/landice/mesh.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 4947158cc1..daf469b1e3 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -1,10 +1,11 @@ import time import jigsawpy +import mpas_tools.io import numpy as np import xarray from geometric_features import FeatureCollection, GeometricFeatures -from mpas_tools.io import default_format, write_netcdf +from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call from mpas_tools.mesh.conversion import convert, cull from mpas_tools.mesh.creation import build_planar_mesh @@ -653,5 +654,6 @@ def make_region_masks(self, mesh_filename, mask_filename, cores, tags): '-o', mask_filename, '-t', 'cell', '--process_count', f'{cores}', - '--format', default_format] + '--format', mpas_tools.io.default_format, + '--engine', mpas_tools.io.default_engine] check_call(args, logger=logger) From 2141f32d34ab9aa0f8f4beda8a92876d9ac48d02 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 11:00:10 -0600 Subject: [PATCH 09/19] Fix culling of Antarctic islands Fix path to gridded dataset to use for building MALI mesh, in order to properly remove disconnected islands from the mesh. --- compass/landice/tests/antarctica/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index dafe705e35..fba7e45b89 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/compass/landice/tests/antarctica/mesh.py @@ -74,7 +74,7 @@ def run(self): build_MALI_mesh( self, cell_width, x1, y1, geom_points, geom_edges, mesh_name=self.mesh_filename, section_name=section_name, - gridded_dataset='antarctica_8km_2020_10_20.nc', + gridded_dataset='antarctica_8km_2020_10_20_floodFillMask.nc', projection='ais-bedmap2', geojson_file=None) logger.info('creating graph.info') From 1817dbb1ce584fc30f7fa5b051c998be603185d7 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 12:44:09 -0600 Subject: [PATCH 10/19] Update documentation for compass.landice.mesh Add compass.landice.mesh to framework documentation, and update some doc-strings. --- compass/landice/mesh.py | 39 ++++++- .../tests/antarctica/mesh_gen/mesh_gen.cfg | 8 ++ docs/developers_guide/landice/framework.rst | 93 ++++++++++++++++ docs/users_guide/landice/framework.rst | 105 ++++++++++++++++++ docs/users_guide/landice/framework/index.rst | 6 - docs/users_guide/landice/index.rst | 2 +- 6 files changed, 244 insertions(+), 9 deletions(-) create mode 100644 docs/users_guide/landice/framework.rst delete mode 100644 docs/users_guide/landice/framework/index.rst diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index daf469b1e3..12541e8263 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -130,8 +130,9 @@ def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. Requires the following options to be set in the given config section: ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, - ``high_dist``, ``low_dist``,``cull_distance``, ``use_speed``, - ``use_dist_to_edge``, and ``use_dist_to_grounding_line``. + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. Parameters ---------- @@ -517,6 +518,40 @@ def build_cell_width(self, section_name, gridded_dataset, def build_MALI_mesh(self, cell_width, x1, y1, geom_points, geom_edges, mesh_name, section_name, gridded_dataset, projection, geojson_file=None): + """ + Create the MALI mesh based on final cell widths determined by + py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and + MPAS-Tools functions. Culls the mesh based on config options, interpolates + all available fields from the gridded dataset to the MALI mesh using the + bilinear method, and marks domain boundaries as Dirichlet cells. + + Parameters + ---------- + cell_width : numpy.ndarray + Desired width of MPAS cells calculated by :py:func:`build_cell_width()` + based on mesh density functions define in :py_func:`set_cell_width()` + to pass to + :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + x1 : float + x coordinates from gridded dataset + y1 : float + y coordinates from gridded dataset + geom_points : jigsawpy.jigsaw_msh_t.VERT2_t + xy node coordinates to pass to build_planar_mesh() + geom_edges : jigsawpy.jigsaw_msh_t.EDGE2_t + xy edge coordinates between nodes to pass to build_planar_mesh() + mesh_name : str + Filename to be used for final MALI .nc mesh file. + section_name : str + Name of config section containing mesh creation options. + gridded_dataset : str + Name of gridded dataset file to be used for interpolation to MALI mesh + projection : str + Projection to be used for setting lat-long fields. + Likely 'gis-gimp' or 'ais-bedmap2' + geojson_file : str + Name of geojson file that defines regional domain extent. + """ logger = self.logger section = self.config[section_name] diff --git a/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg b/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg index 917f706f52..bc9d62530d 100644 --- a/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg @@ -28,10 +28,18 @@ low_log_speed = 0.75 high_dist = 1.e5 # distance to ice edge or grounding line within which cell spacing = min_spac (meters) low_dist = 1.e4 + +# These will not be applied unless use_bed = True. +# They are taken from the humboldt.mesh_gen test case +# and have not been evaluated for Antarctica. # distance at which bed topography has no effect high_dist_bed = 1.e5 # distance within which bed topography has maximum effect low_dist_bed = 7.5e4 +# Bed elev beneath which cell spacing is minimized +low_bed = 50.0 +# Bed elev above which cell spacing is maximized +high_bed = 100.0 # mesh density functions use_speed = True diff --git a/docs/developers_guide/landice/framework.rst b/docs/developers_guide/landice/framework.rst index 2bb6433ca7..d3ac1c7dc5 100644 --- a/docs/developers_guide/landice/framework.rst +++ b/docs/developers_guide/landice/framework.rst @@ -10,3 +10,96 @@ extrapolate The landice framework module ``compass/landice/extrapolate.py`` provides a function for extrapolating variables into undefined regions. It is copied from a similar script in MPAS-Tools. + +mesh +~~~ + +The landice framework module :py:mod:`compass.landice.mesh` provides functions +used by all ``mesh_gen`` tests cases, which currently exist within the +``antarctica``, ``greenland``, ``humboldt``, ``kangerlussuaq``, ``koge_bugt_s``, +and ``thwaites`` test groups. These functions include: + +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.set_rectangular_geom_points_and_edges()` sets node +and edge coordinates to pass to py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + +py:func:`compass.landice.mesh.set_cell_width()` sets cell widths based on settings +in config file to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. +Requires the following options to be set in the given config section: ``min_spac``, +``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, +``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, +``use_speed``, ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + +py:func:`compass.landice.mesh.get_dist_to_edge_and_GL()` calculates distance from +each point to ice edge and grounding line, to be used in mesh density functions in +:py:func:`compass.landice.mesh.set_cell_width()`. In future development, +this should be updated to use a faster package such as `scikit-fmm`. + +py:func:`compass.landice.mesh.build_cell_width()` determine final MPAS mesh cell sizes +using desired cell widths calculated by py:func:`compass.landice.mesh.set_cell_width()`, +based on user-defined density functions and config options. + +py:func:`compass.landice.mesh.build_MALI_mesh()` creates the MALI mesh based on final +cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw +and MPAS-Tools functions. Culls the mesh based on config options, interpolates +all available fields from the gridded dataset to the MALI mesh using the bilinear +method, and marks domain boundaries as Dirichlet cells. + +py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions +defined in Geometric Features repository. It is only used by the ``antarctica`` +and ``greenland`` test cases. + +The following config options should be define for all ``mesh_gen`` test cases (although +not necessarily with the same values shown here, which are the defaults for the 1–10km +Humboldt mesh): + +.. code-block:: cfg + + # config options for humboldt test cases + [mesh] + + # number of levels in the mesh + levels = 10 + + # Bounds of Humboldt domain + x_min = -630000. + x_max = 84000. + y_min = -1560000. + y_max = -860000. + + # distance from ice margin to cull (km). + # Set to a value <= 0 if you do not want + # to cull based on distance from margin. + cull_distance = 5.0 + + # mesh density parameters + # minimum cell spacing (meters) + min_spac = 1.e3 + # maximum cell spacing (meters) + max_spac = 1.e4 + # log10 of max speed (m/yr) for cell spacing + high_log_speed = 2.5 + # log10 of min speed (m/yr) for cell spacing + low_log_speed = 0.75 + # distance at which cell spacing = max_spac (meters) + high_dist = 1.e5 + # distance within which cell spacing = min_spac (meters) + low_dist = 1.e4 + + # These *_bed settings are only applied when use_bed = True. + # distance at which bed topography has no effect + high_dist_bed = 1.e5 + # distance within which bed topography has maximum effect + low_dist_bed = 5.e4 + # Bed elev beneath which cell spacing is minimized + low_bed = 50.0 + # Bed elev above which cell spacing is maximized + high_bed = 100.0 + + # mesh density functions + use_speed = True + use_dist_to_grounding_line = False + use_dist_to_edge = True + use_bed = True diff --git a/docs/users_guide/landice/framework.rst b/docs/users_guide/landice/framework.rst new file mode 100644 index 0000000000..d3ac1c7dc5 --- /dev/null +++ b/docs/users_guide/landice/framework.rst @@ -0,0 +1,105 @@ +.. _dev_landice_framework: + +Land-ice Framework +================== + + +extrapolate +~~~~~~~~~~~ + +The landice framework module ``compass/landice/extrapolate.py`` provides a +function for extrapolating variables into undefined regions. It is copied +from a similar script in MPAS-Tools. + +mesh +~~~ + +The landice framework module :py:mod:`compass.landice.mesh` provides functions +used by all ``mesh_gen`` tests cases, which currently exist within the +``antarctica``, ``greenland``, ``humboldt``, ``kangerlussuaq``, ``koge_bugt_s``, +and ``thwaites`` test groups. These functions include: + +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.set_rectangular_geom_points_and_edges()` sets node +and edge coordinates to pass to py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + +py:func:`compass.landice.mesh.set_cell_width()` sets cell widths based on settings +in config file to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. +Requires the following options to be set in the given config section: ``min_spac``, +``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, +``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, +``use_speed``, ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + +py:func:`compass.landice.mesh.get_dist_to_edge_and_GL()` calculates distance from +each point to ice edge and grounding line, to be used in mesh density functions in +:py:func:`compass.landice.mesh.set_cell_width()`. In future development, +this should be updated to use a faster package such as `scikit-fmm`. + +py:func:`compass.landice.mesh.build_cell_width()` determine final MPAS mesh cell sizes +using desired cell widths calculated by py:func:`compass.landice.mesh.set_cell_width()`, +based on user-defined density functions and config options. + +py:func:`compass.landice.mesh.build_MALI_mesh()` creates the MALI mesh based on final +cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw +and MPAS-Tools functions. Culls the mesh based on config options, interpolates +all available fields from the gridded dataset to the MALI mesh using the bilinear +method, and marks domain boundaries as Dirichlet cells. + +py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions +defined in Geometric Features repository. It is only used by the ``antarctica`` +and ``greenland`` test cases. + +The following config options should be define for all ``mesh_gen`` test cases (although +not necessarily with the same values shown here, which are the defaults for the 1–10km +Humboldt mesh): + +.. code-block:: cfg + + # config options for humboldt test cases + [mesh] + + # number of levels in the mesh + levels = 10 + + # Bounds of Humboldt domain + x_min = -630000. + x_max = 84000. + y_min = -1560000. + y_max = -860000. + + # distance from ice margin to cull (km). + # Set to a value <= 0 if you do not want + # to cull based on distance from margin. + cull_distance = 5.0 + + # mesh density parameters + # minimum cell spacing (meters) + min_spac = 1.e3 + # maximum cell spacing (meters) + max_spac = 1.e4 + # log10 of max speed (m/yr) for cell spacing + high_log_speed = 2.5 + # log10 of min speed (m/yr) for cell spacing + low_log_speed = 0.75 + # distance at which cell spacing = max_spac (meters) + high_dist = 1.e5 + # distance within which cell spacing = min_spac (meters) + low_dist = 1.e4 + + # These *_bed settings are only applied when use_bed = True. + # distance at which bed topography has no effect + high_dist_bed = 1.e5 + # distance within which bed topography has maximum effect + low_dist_bed = 5.e4 + # Bed elev beneath which cell spacing is minimized + low_bed = 50.0 + # Bed elev above which cell spacing is maximized + high_bed = 100.0 + + # mesh density functions + use_speed = True + use_dist_to_grounding_line = False + use_dist_to_edge = True + use_bed = True diff --git a/docs/users_guide/landice/framework/index.rst b/docs/users_guide/landice/framework/index.rst deleted file mode 100644 index 04e2038fae..0000000000 --- a/docs/users_guide/landice/framework/index.rst +++ /dev/null @@ -1,6 +0,0 @@ -.. _landice_framework: - -Framework -========= - -So far, the ``landice`` core has no shared framework. diff --git a/docs/users_guide/landice/index.rst b/docs/users_guide/landice/index.rst index 33efd66aec..e76ab0e4aa 100644 --- a/docs/users_guide/landice/index.rst +++ b/docs/users_guide/landice/index.rst @@ -30,4 +30,4 @@ Some helpful external links: test_groups/index suites - framework/index + framework From 4898a44e97110142b328ca440fb64b5a360e2530 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 13:51:24 -0600 Subject: [PATCH 11/19] Fix a few typos after building docs Fix a few typos that only became apparent after building docs. --- compass/landice/mesh.py | 4 ++-- docs/developers_guide/landice/api.rst | 1 - docs/developers_guide/landice/framework.rst | 18 +++++++++--------- docs/users_guide/landice/framework.rst | 20 ++++++++++---------- 4 files changed, 21 insertions(+), 22 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 12541e8263..4e31569f18 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -520,7 +520,7 @@ def build_MALI_mesh(self, cell_width, x1, y1, geom_points, gridded_dataset, projection, geojson_file=None): """ Create the MALI mesh based on final cell widths determined by - py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and + :py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and MPAS-Tools functions. Culls the mesh based on config options, interpolates all available fields from the gridded dataset to the MALI mesh using the bilinear method, and marks domain boundaries as Dirichlet cells. @@ -529,7 +529,7 @@ def build_MALI_mesh(self, cell_width, x1, y1, geom_points, ---------- cell_width : numpy.ndarray Desired width of MPAS cells calculated by :py:func:`build_cell_width()` - based on mesh density functions define in :py_func:`set_cell_width()` + based on mesh density functions define in :py:func:`set_cell_width()` to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. x1 : float diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index d832d8a16e..b51ba690dd 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -203,7 +203,6 @@ greenland mesh.Mesh mesh.Mesh.run - mesh.Mesh.build_cell_width mesh_gen.MeshGen mesh_gen.MeshGen.run diff --git a/docs/developers_guide/landice/framework.rst b/docs/developers_guide/landice/framework.rst index d3ac1c7dc5..d03aaf7641 100644 --- a/docs/developers_guide/landice/framework.rst +++ b/docs/developers_guide/landice/framework.rst @@ -12,46 +12,46 @@ function for extrapolating variables into undefined regions. It is copied from a similar script in MPAS-Tools. mesh -~~~ +~~~~ The landice framework module :py:mod:`compass.landice.mesh` provides functions used by all ``mesh_gen`` tests cases, which currently exist within the ``antarctica``, ``greenland``, ``humboldt``, ``kangerlussuaq``, ``koge_bugt_s``, and ``thwaites`` test groups. These functions include: -py:func:`compass.landice.mesh.gridded_flood_fill()` applies a flood-fill algorithm +: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.set_rectangular_geom_points_and_edges()` sets node +:py:func:`compass.landice.mesh.set_rectangular_geom_points_and_edges()` sets node and edge coordinates to pass to py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. -py:func:`compass.landice.mesh.set_cell_width()` sets cell widths based on settings +:py:func:`compass.landice.mesh.set_cell_width()` sets cell widths based on settings in config file to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. Requires the following options to be set in the given config section: ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. -py:func:`compass.landice.mesh.get_dist_to_edge_and_GL()` calculates distance from +:py:func:`compass.landice.mesh.get_dist_to_edge_and_GL()` calculates distance from each point to ice edge and grounding line, to be used in mesh density functions in :py:func:`compass.landice.mesh.set_cell_width()`. In future development, this should be updated to use a faster package such as `scikit-fmm`. -py:func:`compass.landice.mesh.build_cell_width()` determine final MPAS mesh cell sizes +:py:func:`compass.landice.mesh.build_cell_width()` determine final MPAS mesh cell sizes using desired cell widths calculated by py:func:`compass.landice.mesh.set_cell_width()`, based on user-defined density functions and config options. -py:func:`compass.landice.mesh.build_MALI_mesh()` creates the MALI mesh based on final +:py:func:`compass.landice.mesh.build_MALI_mesh()` creates the MALI mesh based on final cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and MPAS-Tools functions. Culls the mesh based on config options, interpolates all available fields from the gridded dataset to the MALI mesh using the bilinear method, and marks domain boundaries as Dirichlet cells. -py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions +:py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions defined in Geometric Features repository. It is only used by the ``antarctica`` and ``greenland`` test cases. -The following config options should be define for all ``mesh_gen`` test cases (although +The following config options should be defined for all ``mesh_gen`` test cases (although not necessarily with the same values shown here, which are the defaults for the 1–10km Humboldt mesh): diff --git a/docs/users_guide/landice/framework.rst b/docs/users_guide/landice/framework.rst index d3ac1c7dc5..995ded63d5 100644 --- a/docs/users_guide/landice/framework.rst +++ b/docs/users_guide/landice/framework.rst @@ -1,4 +1,4 @@ -.. _dev_landice_framework: +.. _landice_framework: Land-ice Framework ================== @@ -12,46 +12,46 @@ function for extrapolating variables into undefined regions. It is copied from a similar script in MPAS-Tools. mesh -~~~ +~~~~ The landice framework module :py:mod:`compass.landice.mesh` provides functions used by all ``mesh_gen`` tests cases, which currently exist within the ``antarctica``, ``greenland``, ``humboldt``, ``kangerlussuaq``, ``koge_bugt_s``, and ``thwaites`` test groups. These functions include: -py:func:`compass.landice.mesh.gridded_flood_fill()` applies a flood-fill algorithm +: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.set_rectangular_geom_points_and_edges()` sets node +:py:func:`compass.landice.mesh.set_rectangular_geom_points_and_edges()` sets node and edge coordinates to pass to py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. -py:func:`compass.landice.mesh.set_cell_width()` sets cell widths based on settings +:py:func:`compass.landice.mesh.set_cell_width()` sets cell widths based on settings in config file to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. Requires the following options to be set in the given config section: ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. -py:func:`compass.landice.mesh.get_dist_to_edge_and_GL()` calculates distance from +:py:func:`compass.landice.mesh.get_dist_to_edge_and_GL()` calculates distance from each point to ice edge and grounding line, to be used in mesh density functions in :py:func:`compass.landice.mesh.set_cell_width()`. In future development, this should be updated to use a faster package such as `scikit-fmm`. -py:func:`compass.landice.mesh.build_cell_width()` determine final MPAS mesh cell sizes +:py:func:`compass.landice.mesh.build_cell_width()` determine final MPAS mesh cell sizes using desired cell widths calculated by py:func:`compass.landice.mesh.set_cell_width()`, based on user-defined density functions and config options. -py:func:`compass.landice.mesh.build_MALI_mesh()` creates the MALI mesh based on final +:py:func:`compass.landice.mesh.build_MALI_mesh()` creates the MALI mesh based on final cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and MPAS-Tools functions. Culls the mesh based on config options, interpolates all available fields from the gridded dataset to the MALI mesh using the bilinear method, and marks domain boundaries as Dirichlet cells. -py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions +:py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions defined in Geometric Features repository. It is only used by the ``antarctica`` and ``greenland`` test cases. -The following config options should be define for all ``mesh_gen`` test cases (although +The following config options should be defined for all ``mesh_gen`` test cases (although not necessarily with the same values shown here, which are the defaults for the 1–10km Humboldt mesh): From c62cb8d731913f2f8a17a5f5428d5af8d2dab2d0 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 11 Apr 2023 14:17:06 -0600 Subject: [PATCH 12/19] Update file-change-action to 1.2.4 This should reportedly fix: ``` Error: There was an issue sorting changed files from Github. Exception: { "error": "500/TypeError", "from": "sortChangedFiles", "message": "There was an issue sorting files changed.", "payload": "{}" } ``` See: https://github.com/trilom/file-changes-action/issues/104 --- .github/workflows/build_workflow.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build_workflow.yml b/.github/workflows/build_workflow.yml index ae88a9f7eb..ecd3cd4be7 100644 --- a/.github/workflows/build_workflow.yml +++ b/.github/workflows/build_workflow.yml @@ -36,7 +36,7 @@ jobs: python-version: "3.10" - if: ${{ steps.skip_check.outputs.should_skip != 'true' }} id: file_changes - uses: trilom/file-changes-action@v1.2.3 + uses: trilom/file-changes-action@v1.2.4 with: output: ' ' - if: ${{ steps.skip_check.outputs.should_skip != 'true' }} From 5543fa9aa0f2a1aef4d68f066d2683a9aba6b6a7 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 15:35:23 -0600 Subject: [PATCH 13/19] Change function names to all lowercase Change function names to all lowercase in compass.landice.mesh. --- compass/landice/mesh.py | 8 ++++---- compass/landice/tests/antarctica/mesh.py | 4 ++-- compass/landice/tests/greenland/mesh.py | 4 ++-- compass/landice/tests/humboldt/mesh.py | 4 ++-- compass/landice/tests/kangerlussuaq/mesh.py | 4 ++-- compass/landice/tests/koge_bugt_s/mesh.py | 4 ++-- compass/landice/tests/thwaites/mesh.py | 4 ++-- docs/developers_guide/landice/api.rst | 4 ++-- docs/developers_guide/landice/framework.rst | 4 ++-- docs/users_guide/landice/framework.rst | 4 ++-- 10 files changed, 22 insertions(+), 22 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 4e31569f18..2b30324141 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -300,7 +300,7 @@ def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, return cell_width -def get_dist_to_edge_and_GL(self, thk, topg, x, y, section, window_size=None): +def get_dist_to_edge_and_gl(self, thk, topg, x, y, section, window_size=None): """ Calculate distance from each point to ice edge and grounding line, to be used in mesh density functions in @@ -416,7 +416,7 @@ def get_dist_to_edge_and_GL(self, thk, topg, x, y, section, window_size=None): dist_to_grounding_line[i, j] = dist_to_here_grounding_line.min() toc = time.time() - logger.info('compass.landice.mesh.get_dist_to_edge_and_GL() took {:0.2f} ' + logger.info('compass.landice.mesh.get_dist_to_edge_and_gl() took {:0.2f} ' 'seconds'.format(toc - tic)) return dist_to_edge, dist_to_grounding_line @@ -499,7 +499,7 @@ def build_cell_width(self, section_name, gridded_dataset, # Calculate distance from each grid point to ice edge # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL( + distToEdge, distToGL = get_dist_to_edge_and_gl( self, thk, topg, x1, y1, section=section_name) @@ -515,7 +515,7 @@ def build_cell_width(self, section_name, gridded_dataset, y1.astype('float64'), geom_points, geom_edges, flood_mask) -def build_MALI_mesh(self, cell_width, x1, y1, geom_points, +def build_mali_mesh(self, cell_width, x1, y1, geom_points, geom_edges, mesh_name, section_name, gridded_dataset, projection, geojson_file=None): """ diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index fba7e45b89..f36afe0979 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/compass/landice/tests/antarctica/mesh.py @@ -4,7 +4,7 @@ from compass.landice.mesh import ( build_cell_width, - build_MALI_mesh, + build_mali_mesh, make_region_masks, ) from compass.model import make_graph_file @@ -71,7 +71,7 @@ def run(self): gg.variables['vy'][0, :, :] *= floodFillMask gg.close() - build_MALI_mesh( + build_mali_mesh( self, cell_width, x1, y1, geom_points, geom_edges, mesh_name=self.mesh_filename, section_name=section_name, gridded_dataset='antarctica_8km_2020_10_20_floodFillMask.nc', diff --git a/compass/landice/tests/greenland/mesh.py b/compass/landice/tests/greenland/mesh.py index 40c3ddb828..bb7dc702e9 100644 --- a/compass/landice/tests/greenland/mesh.py +++ b/compass/landice/tests/greenland/mesh.py @@ -1,6 +1,6 @@ from compass.landice.mesh import ( build_cell_width, - build_MALI_mesh, + build_mali_mesh, make_region_masks, ) from compass.model import make_graph_file @@ -57,7 +57,7 @@ def run(self): gridded_dataset='greenland_2km_2020_04_20.epsg3413.nc', flood_fill_start=[100, 700]) - build_MALI_mesh( + build_mali_mesh( self, cell_width, x1, y1, geom_points, geom_edges, mesh_name=mesh_name, section_name=section_name, gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa diff --git a/compass/landice/tests/humboldt/mesh.py b/compass/landice/tests/humboldt/mesh.py index ca9f7b58ee..8c168f8c80 100644 --- a/compass/landice/tests/humboldt/mesh.py +++ b/compass/landice/tests/humboldt/mesh.py @@ -1,4 +1,4 @@ -from compass.landice.mesh import build_cell_width, build_MALI_mesh +from compass.landice.mesh import build_cell_width, build_mali_mesh from compass.model import make_graph_file from compass.step import Step @@ -56,7 +56,7 @@ def run(self): self, section_name=section_name, gridded_dataset='greenland_2km_2020_04_20.epsg3413.nc') - build_MALI_mesh( + build_mali_mesh( self, cell_width, x1, y1, geom_points, geom_edges, mesh_name=mesh_name, section_name=section_name, gridded_dataset='humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', diff --git a/compass/landice/tests/kangerlussuaq/mesh.py b/compass/landice/tests/kangerlussuaq/mesh.py index cab345cfc2..a142c34f79 100644 --- a/compass/landice/tests/kangerlussuaq/mesh.py +++ b/compass/landice/tests/kangerlussuaq/mesh.py @@ -1,4 +1,4 @@ -from compass.landice.mesh import build_cell_width, build_MALI_mesh +from compass.landice.mesh import build_cell_width, build_mali_mesh from compass.model import make_graph_file from compass.step import Step @@ -57,7 +57,7 @@ def run(self): self, section_name=section_name, gridded_dataset='greenland_8km_2020_04_20.epsg3413.nc') - build_MALI_mesh( + build_mali_mesh( self, cell_width, x1, y1, geom_points, geom_edges, mesh_name=mesh_name, section_name=section_name, gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa diff --git a/compass/landice/tests/koge_bugt_s/mesh.py b/compass/landice/tests/koge_bugt_s/mesh.py index 241c829996..4a41fdb36b 100644 --- a/compass/landice/tests/koge_bugt_s/mesh.py +++ b/compass/landice/tests/koge_bugt_s/mesh.py @@ -1,4 +1,4 @@ -from compass.landice.mesh import build_cell_width, build_MALI_mesh +from compass.landice.mesh import build_cell_width, build_mali_mesh from compass.model import make_graph_file from compass.step import Step @@ -57,7 +57,7 @@ def run(self): self, section_name=section_name, gridded_dataset='greenland_8km_2020_04_20.epsg3413.nc') - build_MALI_mesh( + build_mali_mesh( self, cell_width, x1, y1, geom_points, geom_edges, mesh_name=mesh_name, section_name=section_name, gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa diff --git a/compass/landice/tests/thwaites/mesh.py b/compass/landice/tests/thwaites/mesh.py index ce9447376b..863bbd6796 100644 --- a/compass/landice/tests/thwaites/mesh.py +++ b/compass/landice/tests/thwaites/mesh.py @@ -1,4 +1,4 @@ -from compass.landice.mesh import build_cell_width, build_MALI_mesh +from compass.landice.mesh import build_cell_width, build_mali_mesh from compass.model import make_graph_file from compass.step import Step @@ -46,7 +46,7 @@ def run(self): self, section_name=section_name, gridded_dataset='antarctica_8km_2020_10_20.nc') - build_MALI_mesh( + build_mali_mesh( self, cell_width, x1, y1, geom_points, geom_edges, mesh_name=mesh_name, section_name=section_name, gridded_dataset='antarctica_1km_2020_10_20_ASE.nc', diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index b51ba690dd..9a51322d0a 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -418,7 +418,7 @@ Landice Framework mesh.gridded_flood_fill mesh.set_rectangular_geom_points_and_edges mesh.set_cell_width - mesh.get_dist_to_edge_and_GL + mesh.get_dist_to_edge_and_gl mesh.build_cell_width - mesh.build_MALI_mesh + mesh.build_mali_mesh mesh.make_region_masks diff --git a/docs/developers_guide/landice/framework.rst b/docs/developers_guide/landice/framework.rst index d03aaf7641..2e5c9161e0 100644 --- a/docs/developers_guide/landice/framework.rst +++ b/docs/developers_guide/landice/framework.rst @@ -32,7 +32,7 @@ Requires the following options to be set in the given config section: ``min_spac ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. -:py:func:`compass.landice.mesh.get_dist_to_edge_and_GL()` calculates distance from +:py:func:`compass.landice.mesh.get_dist_to_edge_and_gl()` calculates distance from each point to ice edge and grounding line, to be used in mesh density functions in :py:func:`compass.landice.mesh.set_cell_width()`. In future development, this should be updated to use a faster package such as `scikit-fmm`. @@ -41,7 +41,7 @@ this should be updated to use a faster package such as `scikit-fmm`. using desired cell widths calculated by py:func:`compass.landice.mesh.set_cell_width()`, based on user-defined density functions and config options. -:py:func:`compass.landice.mesh.build_MALI_mesh()` creates the MALI mesh based on final +:py:func:`compass.landice.mesh.build_mali_mesh()` creates the MALI mesh based on final cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and MPAS-Tools functions. Culls the mesh based on config options, interpolates all available fields from the gridded dataset to the MALI mesh using the bilinear diff --git a/docs/users_guide/landice/framework.rst b/docs/users_guide/landice/framework.rst index 995ded63d5..d6e0e52f81 100644 --- a/docs/users_guide/landice/framework.rst +++ b/docs/users_guide/landice/framework.rst @@ -32,7 +32,7 @@ Requires the following options to be set in the given config section: ``min_spac ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. -:py:func:`compass.landice.mesh.get_dist_to_edge_and_GL()` calculates distance from +:py:func:`compass.landice.mesh.get_dist_to_edge_and_gl()` calculates distance from each point to ice edge and grounding line, to be used in mesh density functions in :py:func:`compass.landice.mesh.set_cell_width()`. In future development, this should be updated to use a faster package such as `scikit-fmm`. @@ -41,7 +41,7 @@ this should be updated to use a faster package such as `scikit-fmm`. using desired cell widths calculated by py:func:`compass.landice.mesh.set_cell_width()`, based on user-defined density functions and config options. -:py:func:`compass.landice.mesh.build_MALI_mesh()` creates the MALI mesh based on final +:py:func:`compass.landice.mesh.build_mali_mesh()` creates the MALI mesh based on final cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and MPAS-Tools functions. Culls the mesh based on config options, interpolates all available fields from the gridded dataset to the MALI mesh using the bilinear From d2ba3d738efd596b9f5538a09f1114283d4411f2 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 16:06:58 -0600 Subject: [PATCH 14/19] Update Users and Developers guides after code review Remove extraneous information from Users Guide. Add explanation of defining x/y bounds in example config file. --- docs/developers_guide/landice/framework.rst | 4 +- docs/users_guide/landice/framework.rst | 41 ++------------------- 2 files changed, 6 insertions(+), 39 deletions(-) diff --git a/docs/developers_guide/landice/framework.rst b/docs/developers_guide/landice/framework.rst index 2e5c9161e0..f8c0fc858c 100644 --- a/docs/developers_guide/landice/framework.rst +++ b/docs/developers_guide/landice/framework.rst @@ -63,7 +63,9 @@ Humboldt mesh): # number of levels in the mesh levels = 10 - # Bounds of Humboldt domain + # Bounds of Humboldt domain. If you want the extent + # of the gridded dataset to determine the extent of + # the MALI domain, set these to None. x_min = -630000. x_max = 84000. y_min = -1560000. diff --git a/docs/users_guide/landice/framework.rst b/docs/users_guide/landice/framework.rst index d6e0e52f81..89222b9a51 100644 --- a/docs/users_guide/landice/framework.rst +++ b/docs/users_guide/landice/framework.rst @@ -14,43 +14,6 @@ from a similar script in MPAS-Tools. mesh ~~~~ -The landice framework module :py:mod:`compass.landice.mesh` provides functions -used by all ``mesh_gen`` tests cases, which currently exist within the -``antarctica``, ``greenland``, ``humboldt``, ``kangerlussuaq``, ``koge_bugt_s``, -and ``thwaites`` test groups. These functions include: - -: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.set_rectangular_geom_points_and_edges()` sets node -and edge coordinates to pass to py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. - -:py:func:`compass.landice.mesh.set_cell_width()` sets cell widths based on settings -in config file to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. -Requires the following options to be set in the given config section: ``min_spac``, -``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, -``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, -``use_speed``, ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. - -:py:func:`compass.landice.mesh.get_dist_to_edge_and_gl()` calculates distance from -each point to ice edge and grounding line, to be used in mesh density functions in -:py:func:`compass.landice.mesh.set_cell_width()`. In future development, -this should be updated to use a faster package such as `scikit-fmm`. - -:py:func:`compass.landice.mesh.build_cell_width()` determine final MPAS mesh cell sizes -using desired cell widths calculated by py:func:`compass.landice.mesh.set_cell_width()`, -based on user-defined density functions and config options. - -:py:func:`compass.landice.mesh.build_mali_mesh()` creates the MALI mesh based on final -cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw -and MPAS-Tools functions. Culls the mesh based on config options, interpolates -all available fields from the gridded dataset to the MALI mesh using the bilinear -method, and marks domain boundaries as Dirichlet cells. - -:py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions -defined in Geometric Features repository. It is only used by the ``antarctica`` -and ``greenland`` test cases. - The following config options should be defined for all ``mesh_gen`` test cases (although not necessarily with the same values shown here, which are the defaults for the 1–10km Humboldt mesh): @@ -63,7 +26,9 @@ Humboldt mesh): # number of levels in the mesh levels = 10 - # Bounds of Humboldt domain + # Bounds of Humboldt domain. If you want the extent + # of the gridded dataset to determine the extent of + # the MALI domain, set these to None. x_min = -630000. x_max = 84000. y_min = -1560000. From ae88c40091f02f847496bec58315031bc12a99f5 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 16:27:00 -0600 Subject: [PATCH 15/19] Clean up doc-strings and add required config options Clean up doc-strings and add a list of required config options for each function in compass/landice/mesh.py that requires section_name as a parameter. Refer to the Users and Developers guides for more information on the meaning of the config options. --- compass/landice/mesh.py | 139 +++++++++++++++++++++++++++++----------- 1 file changed, 101 insertions(+), 38 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 2b30324141..f7ed2376a8 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -25,9 +25,11 @@ def gridded_flood_fill(field, iStart=None, jStart=None): field : numpy.ndarray Array from gridded dataset to use for flood-fill. Usually ice thickness. + iStart : int x index from which to start flood fill for field. Defaults to the center x coordinate. + jStart : int y index from which to start flood fill. Defaults to the center y coordinate. @@ -35,7 +37,7 @@ def gridded_flood_fill(field, iStart=None, jStart=None): Returns ------- flood_mask : numpy.ndarray - _mask calculated by the flood fill routine, + mask calculated by the flood fill routine, where cells connected to the ice sheet (or main feature) are 1 and everything else is 0. """ @@ -90,19 +92,23 @@ def set_rectangular_geom_points_and_edges(xmin, xmax, ymin, ymax): ---------- xmin : int or float Left-most x-coordinate in region to mesh + xmax : int or float Right-most x-coordinate in region to mesh + ymin : int or float Bottom-most y-coordinate in region to mesh + ymax : int or float Top-most y-coordinate in region to mesh Returns ------- geom_points : jigsawpy.jigsaw_msh_t.VERT2_t - xy node coordinates to pass to build_planar_mesh() + xy node coordinates to pass to ``build_planar_mesh()`` + geom_edges : jigsawpy.jigsaw_msh_t.EDGE2_t - xy edge coordinates between nodes to pass to build_planar_mesh() + xy edge coordinates between nodes to pass to ``build_planar_mesh()`` """ geom_points = np.array([ # list of xy "node" coordinates @@ -122,46 +128,56 @@ def set_rectangular_geom_points_and_edges(xmin, xmax, ymin, ymax): return geom_points, geom_edges -def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, +def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None, dist_to_edge=None, dist_to_grounding_line=None, flood_fill_iStart=None, flood_fill_jStart=None): """ Set cell widths based on settings in config file to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. - Requires the following options to be set in the given config section: - ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, - ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, - ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, - ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. Parameters ---------- - section : str - Section of the config file from which to read parameters + section_name : str + Section of the config file from which to read parameters. The + following options to be set in the given config section: + ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + See the Land-Ice Framework section of the Users or Developers guide + for more information about these options and their uses. + thk : numpy.ndarray Ice thickness field from gridded dataset, usually after trimming to flood fill mask + bed : numpy.ndarray Bed topography from gridded dataset + vx : numpy.ndarray, optional x-component of ice velocity from gridded dataset, usually after trimming to flood fill mask. Can be set to ``None`` if ``use_speed == False`` in config file. + vy : numpy.ndarray, optional y-component of ice velocity from gridded dataset, usually after trimming to flood fill mask. Can be set to ``None`` if ``use_speed == False`` in config file. + dist_to_edge : numpy.ndarray, optional Distance from each cell to ice edge, calculated in separate function. Can be set to ``None`` if ``use_dist_to_edge == False`` in config file and you do not want to set large ``cell_width`` where cells will be culled anyway, but this is not recommended. + dist_to_grounding_line : numpy.ndarray, optional Distance from each cell to grounding line, calculated in separate function. Can be set to ``None`` if ``use_dist_to_grounding_line == False`` in config file. + flood_fill_iStart : int, optional x-index location to start flood-fill when using bed topography + flood_fill_jStart : int, optional y-index location to start flood-fill when using bed topography @@ -173,7 +189,7 @@ def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, """ logger = self.logger - section = self.config[section] + section = self.config[section_name] # Get config inputs for cell spacing functions min_spac = section.getfloat('min_spac') @@ -300,44 +316,58 @@ def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, return cell_width -def get_dist_to_edge_and_gl(self, thk, topg, x, y, section, window_size=None): +def get_dist_to_edge_and_gl(self, thk, topg, x, y, + section_name, window_size=None): """ Calculate distance from each point to ice edge and grounding line, to be used in mesh density functions in :py:func:`compass.landice.mesh.set_cell_width()`. In future development, - this should be updated to use a faster package such as `scikit-fmm`. + this should be updated to use a faster package such as ``scikit-fmm``. Parameters ---------- thk : numpy.ndarray Ice thickness field from gridded dataset, usually after trimming to flood fill mask + topg : numpy.ndarray Bed topography field from gridded dataset + x : numpy.ndarray x coordinates from gridded dataset + y : numpy.ndarray y coordinates from gridded dataset - section : str - section of config file used to define mesh parameters + + section_name : str + Section of the config file from which to read parameters. The + following options to be set in the given config section: + ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + See the Land-Ice Framework section of the Users or Developers guide + for more information about these options and their uses. + window_size : int or float Size (in meters) of a search 'box' (one-directional) to use to calculate the distance from each cell to the ice margin. Bigger number makes search slower, but if too small, the transition zone could get truncated. We usually want this calculated as the - maximum of high_dist and high_dist_bed, but there may be cases in - which it is useful to set it manually. However, it should never be - smaller than either high_dist or high_dist_bed. + maximum of ``high_dist`` and ``high_dist_bed``, but there may be cases + in which it is useful to set it manually. However, it should never be + smaller than either ``high_dist`` or ``high_dist_bed``. Returns ------- dist_to_edge : numpy.ndarray Distance from each cell to the ice edge + dist_to_grounding_line : numpy.ndarray Distance from each cell to the grounding line """ logger = self.logger - section = self.config[section] + section = self.config[section_name] tic = time.time() high_dist = float(section.get('high_dist')) @@ -429,28 +459,42 @@ def build_cell_width(self, section_name, gridded_dataset, Parameters ---------- - section : str - section of config file used to define mesh parameters + section_name : str + Section of the config file from which to read parameters. The + following options to be set in the given config section: + ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + See the Land-Ice Framework section of the Users or Developers guide + for more information about these options and their uses. + gridded_dataset : str - name of .nc file used to define cell spacing + name of NetCDF file used to define cell spacing + flood_fill_start : list of ints - i and j indices used to define starting location for flood fill. - Most cases will use [None, None], which will just start the flood - fill in the center of the gridded dataset. + ``i`` and ``j`` indices used to define starting location for flood + fill. Most cases will use ``[None, None]``, which will just start the + flood fill in the center of the gridded dataset. Returns ------- cell_width : numpy.ndarray Desired width of MPAS cells based on mesh desnity functions to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + x1 : float x coordinates from gridded dataset + y1 : float y coordinates from gridded dataset + geom_points : jigsawpy.jigsaw_msh_t.VERT2_t - xy node coordinates to pass to build_planar_mesh() + xy node coordinates to pass to ``build_planar_mesh()`` + geom_edges : jigsawpy.jigsaw_msh_t.EDGE2_t - xy edge coordinates between nodes to pass to build_planar_mesh() + xy edge coordinates between nodes to pass to ``build_planar_mesh()`` + flood_mask : numpy.ndarray mask calculated by the flood fill routine, where cells connected to the ice sheet (or main feature) @@ -501,10 +545,10 @@ def build_cell_width(self, section_name, gridded_dataset, # and grounding line, for use in cell spacing functions. distToEdge, distToGL = get_dist_to_edge_and_gl( self, thk, topg, x1, - y1, section=section_name) + y1, section_name=section_name) # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, section=section_name, + cell_width = set_cell_width(self, section_name=section_name, thk=thk, bed=topg, vx=vx, vy=vy, dist_to_edge=distToEdge, dist_to_grounding_line=distToGL, @@ -532,23 +576,39 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, based on mesh density functions define in :py:func:`set_cell_width()` to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + x1 : float x coordinates from gridded dataset + y1 : float y coordinates from gridded dataset + geom_points : jigsawpy.jigsaw_msh_t.VERT2_t - xy node coordinates to pass to build_planar_mesh() + xy node coordinates to pass to ``build_planar_mesh()`` + geom_edges : jigsawpy.jigsaw_msh_t.EDGE2_t - xy edge coordinates between nodes to pass to build_planar_mesh() + xy edge coordinates between nodes to pass to ``build_planar_mesh()`` + mesh_name : str - Filename to be used for final MALI .nc mesh file. + Filename to be used for final MALI NetCDF mesh file. + section_name : str - Name of config section containing mesh creation options. + Section of the config file from which to read parameters. The + following options to be set in the given config section: + ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + See the Land-Ice Framework section of the Users or Developers guide + for more information about these options and their uses. + gridded_dataset : str Name of gridded dataset file to be used for interpolation to MALI mesh + projection : str Projection to be used for setting lat-long fields. - Likely 'gis-gimp' or 'ais-bedmap2' + Likely ``'gis-gimp'`` or ``'ais-bedmap2'`` + geojson_file : str Name of geojson file that defines regional domain extent. """ @@ -656,16 +716,19 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, def make_region_masks(self, mesh_filename, mask_filename, cores, tags): """ Create masks for ice-sheet subregions based on data - in MPAS-Dev/geometric_fatures. + in ``MPAS-Dev/geometric_fatures``. Parameters ---------- mesh_filename : str - name of .nc mesh file for which to create region masks + name of NetCDF mesh file for which to create region masks + mask_filename : str - name of .nc file to contain region masks + name of NetCDF file to contain region masks + cores : int number of processors used to create region masks + tags : list of str Groups of regions for which masks are to be defined """ From 95af2fff603062fd4f32e02758e6c58a8fc06851 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 16:41:23 -0600 Subject: [PATCH 16/19] Add a few more required config parameters to doc-strings. Add a few more required config parameters to doc-strings. --- compass/landice/mesh.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index f7ed2376a8..4c0bfaac7a 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -140,6 +140,7 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None, section_name : str Section of the config file from which to read parameters. The following options to be set in the given config section: + ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``, ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, @@ -342,6 +343,7 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y, section_name : str Section of the config file from which to read parameters. The following options to be set in the given config section: + ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``, ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, @@ -462,6 +464,7 @@ def build_cell_width(self, section_name, gridded_dataset, section_name : str Section of the config file from which to read parameters. The following options to be set in the given config section: + ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``, ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, @@ -595,6 +598,7 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, section_name : str Section of the config file from which to read parameters. The following options to be set in the given config section: + ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``, ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, From 60c8d6a06b07126a80ac77f3dcf4d0f53f13c9b9 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 16:55:51 -0600 Subject: [PATCH 17/19] Allow a mix of user-defined and dataset-defined bounds Allow the user to define the bounds around a regional domain with any combination of floats and 'None'. For example, a config file that contained ``x_min = 'None'; x_max = 1.0e6; y_min = -1.0e6; y_max = 'None'`` would use you left and top boundaries of the gridded dataset, but also the user-defined values for the right and bottom boundaries. --- compass/landice/mesh.py | 28 +++++++++------------------- 1 file changed, 9 insertions(+), 19 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 4c0bfaac7a..5aa7780cf2 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -518,25 +518,15 @@ def build_cell_width(self, section_name, gridded_dataset, f.close() - # Define extent of region to mesh. - xx0 = section.get('x_min') - xx1 = section.get('x_max') - yy0 = section.get('y_min') - yy1 = section.get('y_max') - - if 'None' in [xx0, xx1, yy0, yy1]: - xx0 = np.min(x1) - xx1 = np.max(x1) - yy0 = np.min(y1) - yy1 = np.max(y1) - else: - xx0 = float(xx0) - xx1 = float(xx1) - yy0 = float(yy0) - yy1 = float(yy1) - - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) + # Get bounds defined by user, or use bound of gridded dataset + bnds = [np.min(x1), np.max(x1), np.min(y1), np.max(y1)] + bnds_options = ['x_min', 'x_max', 'y_min', 'y_max'] + for index, option in enumerate(bnds_options): + bnd = section.get(option) + if bnd != 'None': + bnds[index] = float(bnd) + + geom_points, geom_edges = set_rectangular_geom_points_and_edges(*bnds) # Remove ice not connected to the ice sheet. flood_mask = gridded_flood_fill(thk) From 8b327bd224c27d3590e247002588469544701c51 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 11 Apr 2023 17:02:52 -0600 Subject: [PATCH 18/19] Small clean-up from code review Remove redundant use of ``flood_mask == 0`` Co-authored-by: Xylar Asay-Davis --- compass/landice/mesh.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 5aa7780cf2..12f0bcf12b 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -529,10 +529,10 @@ def build_cell_width(self, section_name, gridded_dataset, geom_points, geom_edges = set_rectangular_geom_points_and_edges(*bnds) # Remove ice not connected to the ice sheet. - flood_mask = gridded_flood_fill(thk) - thk[flood_mask == 0] = 0.0 - vx[flood_mask == 0] = 0.0 - vy[flood_mask == 0] = 0.0 + flood_mask = gridded_flood_fill(thk) == 0 + thk[flood_mask] = 0.0 + vx[flood_mask] = 0.0 + vy[flood_mask] = 0.0 # Calculate distance from each grid point to ice edge # and grounding line, for use in cell spacing functions. From 342856f7fb82c8de71007c22b7ca73ea3e0ffa2e Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 11 Apr 2023 17:12:32 -0600 Subject: [PATCH 19/19] Remove unnecessary calls to logger.info() Remove unnecessary calls to logger.info(), as check_call already takes care of writing the commands to the log file. --- compass/landice/mesh.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index 12f0bcf12b..942405e346 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -621,14 +621,13 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, logger.info('writing grid_converted.nc') write_netcdf(dsMesh, 'grid_converted.nc') levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', 'grid_converted.nc', '-o', 'grid_preCull.nc', '-l', levels, '-v', 'glimmer'] + check_call(args, logger=logger) - logger.info('calling interpolate_to_mpasli_grid.py') args = ['interpolate_to_mpasli_grid.py', '-s', gridded_dataset, '-d', 'grid_preCull.nc', '-m', 'b', '-t'] @@ -637,7 +636,6 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, cullDistance = section.get('cull_distance') if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') args = ['define_cullMask.py', '-f', 'grid_preCull.nc', '-m', 'distance', '-d', cullDistance] @@ -650,13 +648,11 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, if geojson_file is not None: # This step is only necessary because the GeoJSON region # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') args = ['set_lat_lon_fields_in_planar_grid.py', '-f', 'grid_preCull.nc', '-p', projection] check_call(args, logger=logger) - logger.info('calling MpasMaskCreator.x') args = ['MpasMaskCreator.x', 'grid_preCull.nc', 'mask.nc', '-f', geojson_file] @@ -675,6 +671,7 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, logger.info('Marking horns for culling') args = ['mark_horns_for_culling.py', '-f', 'culled.nc'] + check_call(args, logger=logger) logger.info('culling and converting') @@ -683,7 +680,6 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, dsMesh = convert(dsMesh, logger=logger) write_netcdf(dsMesh, 'dehorned.nc') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', 'dehorned.nc', '-o', mesh_name, '-l', levels, '-v', 'glimmer', @@ -691,9 +687,9 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, check_call(args, logger=logger) - logger.info('calling interpolate_to_mpasli_grid.py') args = ['interpolate_to_mpasli_grid.py', '-s', gridded_dataset, '-d', mesh_name, '-m', 'b'] + check_call(args, logger=logger) logger.info('Marking domain boundaries dirichlet') @@ -701,7 +697,6 @@ def build_mali_mesh(self, cell_width, x1, y1, geom_points, '-f', mesh_name] 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', mesh_name, '-p', projection] check_call(args, logger=logger)