Coastal Tools

The mpas_tools.ocean.coastal_tools module contains several functions that aid in the construction of meshes refined around coastlines.

Refining a Mesh

The main driver of coastal tools is the function mpas_tools.ocean.coastal_tools.coastal_refined_mesh(). This function is called repeatedly with different dictionaries of params to add refinement in different locations, starting from a background mesh (see Creating a Background Mesh).

coastal_tools provides the following dictionary of default parameters as a starting point for params:

default_params = {

    # Path to bathymetry data and name of file
    "nc_file": "./earth_relief_15s.nc",
    "nc_vars": ["lon","lat","z"],

    # Bounding box of coastal refinement region
    "region_box": Continental_US,
    "origin": np.array([-100, 40]),
    "restrict_box": Empty,

    # Coastline extraction parameters
    "z_contour": 0.0,
    "n_longest": 10,
    "smooth_coastline": 0,
    "point_list": None,

    # Global mesh parameters
    "grd_box": Entire_Globe,
    "ddeg": .1,
    # 'EC' (defaults to 60to30), 'QU' (uses dx_max_global),
    # 'RRS' (uses dx_max_global and dx_min_global)
    "mesh_type": 'EC',
    "dx_max_global": 30.0 * km,
    "dx_min_global": 10.0 * km,

    # Coastal mesh parameters
    "dx_min_coastal": 10.0 * km,
    "trans_width": 600.0 * km,
    "trans_start": 400.0 * km,

    # Bounding box of plotting region
    "plot_box": North_America,

    # Options
    "nn_search": "kdtree",
    "plot_option": True

}

A coastal mesh is defined by:

Creating a Background Mesh
Extracting Coastlines
Computing Distance to Coast
Blending Cell Widths

Steps 2-4 can be repeated with several different regions to add more refinement.

The first time coastal_refined_mesh() is called, no cell width, latitude or longitude coordinates are supplied. In this case, a background mesh resolution is defined through a call to mpas_tools.ocean.coastal_tools.create_background_mesh(). The following entries in params are used as arguments:

'grd_box' - the bounds of the grid defining the mesh resolution
'ddeg' - the resolution in longitude and latitude in degrees
'mesh_type' - one of {'QU', 'EC' or 'RRS'} indicating the type of mesh
'dx_min_global' - the resolution for a QU mesh, the minimum resolution for an RRS mesh, ignored for EC meshes.
'dx_max_global' - the maximum resolution for an RRS mesh, ignored for QU and EC meshes.
'plot_option' - whether to plot the results and save it to a PNG file.
'plot_box' - If plot_option is True, the bounds of the plot in longitude and latitude.

After the background field of cell widths has been created or using the cell widths passed in as function arguments, coastal_refined_mesh() then adds a region of refinement.

First, a set of coastline contours is extracted by calling mpas_tools.ocean.coastal_tools.extract_coastlines() with the following values from params:

'nc_file' - A bathymetry dataset on a lon/lat grid in NetCDF format
'nc_vars' - The names of the lon, lat and bathymetry variables in a list
'region_box' - see Regions
'z_contour' - A contour level from the bathymetry dataset to extract as a region boundary
'n_longest' - The maximum number of contours to retain (sorted from longest to shortest).
'point_list' - An optional list of points to add to the coastline
'plot_option' - Whether to plot the extracted coastline.
'plot_box' - If plot_option is True, the bounds of the plot in longitude and latitude.

Next, the distance to the coastal contours is computed using mpas_tools.ocean.coastal_tools.distance_to_coast() with the following values from params:

'nn_search' - currently, only the 'kdtree' algorithm is supported
'smooth_coastline' - The number of neighboring cells along the coastline over which to average locations to smooth the coastline
'plot_option' - Whether to plot the distance function.
'plot_box' - If plot_option is True, the bounds of the plot in longitude and latitude.

Finally, the distance function is used to blend the background and refined regions using mpas_tools.ocean.coastal_tools.compute_cell_width() with the following values from params:

'dx_min_coastal' - the resolution in meters of the refined region
'trans_start' - the distance in meters from the coast at which the transition in resolution begins—the center of the transition is half a trans_width farther from the coastline
'trans_width' - the distance in meters over which the transition occurs
'restrict_box' - A region of made up of quadrilaterals to include and exclude that defines where resolution may be altered. Outside of the restrict_box, the resolution remains unchanged. See Regions.
'plot_option' - Whether to plot the cell widths and transition function.
'plot_box' - If plot_option is True, the bounds of the plot in longitude and latitude.

Here is an example of multiple calls to coastal_refined_mesh() in action, taken from the hurricane/USDEQU120at30cr10rr2 test case from COMPASS. This workflow refines a background uniform mesh with 120-km resolution with successively higher and higher resolution down to the Delaware Bay at 2-km resolution.

import mpas_tools.ocean.coastal_tools as ct


km = 1000.0

params = ct.default_params

print("****QU 120 background mesh and enhanced Atlantic (30km)****")
params["mesh_type"] = "QU"
params["dx_max_global"] = 120.0 * km
params["region_box"] = ct.Atlantic
params["restrict_box"] = ct.Atlantic_restrict
params["plot_box"] = ct.Western_Atlantic
params["dx_min_coastal"] = 30.0 * km
params["trans_width"] = 5000.0 * km
params["trans_start"] = 500.0 * km

cell_width, lon, lat = ct.coastal_refined_mesh(params)

print("****Northeast refinement (10km)***")
params["region_box"] = ct.Delaware_Bay
params["plot_box"] = ct.Western_Atlantic
params["dx_min_coastal"] = 10.0 * km
params["trans_width"] = 600.0 * km
params["trans_start"] = 400.0 * km

cell_width, lon, lat = ct.coastal_refined_mesh(
    params, cell_width, lon, lat)

print("****Delaware regional refinement (5km)****")
params["region_box"] = ct.Delaware_Region
params["plot_box"] = ct.Delaware
params["dx_min_coastal"] = 5.0 * km
params["trans_width"] = 175.0 * km
params["trans_start"] = 75.0 * km

cell_width, lon, lat = ct.coastal_refined_mesh(
    params, cell_width, lon, lat)

print("****Delaware Bay high-resolution (2km)****")
params["region_box"] = ct.Delaware_Bay
params["plot_box"] = ct.Delaware
params["restrict_box"] = ct.Delaware_restrict
params["dx_min_coastal"] = 2.0 * km
params["trans_width"] = 100.0 * km
params["trans_start"] = 17.0 * km

cell_width, lon, lat = ct.coastal_refined_mesh(
    params, cell_width, lon, lat)

Creating a Background Mesh

A background mesh is typically created by calling mpas_tools.ocean.coastal_tools.coastal_refined_mesh() without providing an input mesh but can also be created by calling mpas_tools.ocean.coastal_tools.create_background_mesh() directly.

The user must define the bounds of the grid in longitude and latitude (typically -180 to 180 and -90 to 90, respectively) and the resolution in degrees. The mesh can be any of three types: {'QU', 'EC' or 'RRS'}. For Quasi-Uniform (QU) meshes, the resulting cell width will be a constant equal to dx_min_global. For Eddy-Closure (EC) meshes, the default parameters are always used (see Defining an Eddy-closure Mesh). For Rossby-Radius Scaled (RRS) meshes, dx_min_global is the resolution at the poles while dx_max_global is the resolution at the equator (see Defining a Rossby-radius Mesh).

Extracting Coastlines

coastal_tools extracts points along a coastline using mpas_tools.ocean.coastal_tools.extract_coastlines(). The default parameters are set up to use the earth_relief_15s.nc dataset, but any bathymetry data set on a lon/lat grid could be used as long as params['nc_file'] is modified to point to the new name of the dataset and params['nc_vars'] is set to the appropriate variable names.

By default, the coastline is extracted using the z = 0.0 contour of the bathyemtry but other values can be selected (e.g. to use distance from the continental shelf break) by defining params['z_contour']. By default, only the 10 longest contours are retained to reduce computational cost but more (or fewer) contours can be retained by setting params['n_longest'] to another number.

Optionally, the results can be plotted withing the given “plot box” and saved to a file.

Computing Distance to Coast

A key ingredient in defining resolution in coastal meshes is a field containing the distance from each location in the field to the nearest point on the coastline. This distance field D is computed with mpas_tools.ocean.coastal_tools.distance_to_coast() The user could previouly control the search algorithm used via params['nn_search'] but 'kdtree' is now the only option. They can also decide to smooth the coastline as long as there is a single coastline contour—with multiple contours, the current algorithm will average the end of one contour with the start fo the next—by specifying an integer number of neighbors as params['smooth_coastline']. The default is no smoothing (0 neighbors).

Blending Cell Widths

The final step in each iteration of coastal refinement is to blend the new, refined resolution into the previous grid of cell widths with the function mpas_tools.ocean.coastal_tools.compute_cell_width().

The most important parameters to set are params['dx_min_coastal'], the resolution in meters of the refined region; params['trans_start'], the distance from the coast in meters where the transition in resolution should start; and params['trans_width'], the width of the transition itself in meters.

The resolution refinement can be confined to a region using params['restrict_box'] to supply a region of made up of quadrilaterals to include and exclude from the restricted region. include boxes specify regions where the distance-based coastal refinement function should be used to update the desired cell width values. exclude boxes eliminate areas within the include regions so they are not affected by the coastal refinement function. This is useful for preventing resolution from appearing on one side of a narrow piece of land. For example, coastal refinement in the Gulf of Mexico should not appear on the other Pacific Ocean side of Central America. The params['restrict_box'] regions can be used to enforce this.

Here is an example of a restriction box for the region around Delaware, used in the example above:

Delaware_restrict = {"include": [np.array([[-75.853, 39.732],
                                           [-74.939, 36.678],
                                           [-71.519, 40.156],
                                           [-75.153, 40.077]]),
                                 np.array([[-76.024, 37.188],
                                           [-75.214, 36.756],
                                           [-74.512, 37.925],
                                           [-75.274, 38.318]])],
                     "exclude": []}

Regions

Extracting Coastlines requires a set of bounding regions to be defined. These regions are made up of a list of quadrilaterals to include and another list to exclude. The quadrilaterals are either bounding rectangles (min lon, max lon, min lat, max lat) or lists of 4 (lon, lat) points numbered counter clockwise.

include boxes specify areas where coastline contours will be extracted from the underlying bathymetry dataset. exclude boxes can be used to select regions within include regions where coastline extraction should not be done. For example, a large include box covering the U.S. East Coast may contain small islands in the Atlantic, which may need to be ignored when placing coastal refinement. In this case, exclude boxes can be specified to eliminate the small coastlines.

An example of such a region is:

Greenland = {"include":[np.array([-81.5, -12.5, 60, 85])],
             "exclude":[np.array([[-87.6, 58.7],
                                  [-51.9, 56.6],
                                  [-68.9, 75.5],
                                  [-107.0, 73.3]]),
                        np.array([[-119.0, 74.5],
                                  [-92.7, 75.9],
                                  [-52.84, 83.25],
                                  [-100.8, 84.0]]),
                        np.array([[-101.3, 68.5],
                                  [-82.4, 72.3],
                                  [-68.7, 81.24],
                                  [-117.29, 77.75]]),
                        np.array([-25.0, -10.0, 62.5, 67.5])]}

Note

This example includes both bounding rectangles (e.g. np.array([-81.5, -12.5, 60, 85])) and more general quadrilaterals (e.g. np.array([[-101.3, 68.5], [-82.4, 72.3],...)

coastal_tools defines 16 regions via dictionaries of this type. The defined regions are:

Delaware_Bay
Galveston_Bay
Delaware_Region
US_East_Coast
US_Gulf_Coast
Caribbean
US_West_Coast
Hawaii
Alaska
Bering_Sea_E
Bering_Sea_W
Aleutian_Islands_E
Aleutian_Islands_W
Greenland
CONUS
Continental_US