Marching Front algorithm

This Jupyter notebook can be found in this collection of Clawpack apps as the file $CLAW/apps/notebooks/geoclaw/MarchingFront.ipynb. To run this notebook, install Clawpack, and clone the apps repository. A static view of this and other notebooks can be found in the Clawpack Gallery of Jupyter notebooks.

The module clawpack.geoclaw.marching_front.py defines a function select_by_flooding that takes as input an array Ztopo containing topography elevations on a rectangular grid and returns an array pt_chosen of the same shape with values 1 (if chosen) or 0 (if not chosen). Other inputs specify the criteria used to choose points, as described below.

The basic idea is that chosen points satisfy certain elevation requirements along with connectivity to the coast. This was originally developed to identify points in a topography DEM where the topography value $Z$ is below MHW, but that should be initialized as dry land because they are in regions protected by dikes or levies. In this situation, the marching algorithm is used by initializing points well offshore (e.g. with $Z < -5$ meters) as wet and other points to unset. Then the front between wet/unset points is advanced by marking neighboring points as dry if $Z\geq0$ or as wet if $Z<0$. This is repeated iteratively for each new front until there are no more wet points with unset neighbors. At this point any points still unset are entirely buffered by dry points and must lie behind dikes, so these are also set to dry. See the notebook ForceDry.ipynb for more about this application.

Other applications are also described below along with some examples.

Function arguments

The main function in the marching_front.py module is select_by_flooding.

The function is defined by:

def select_by_flooding(Ztopo, mask=None, prev_pts_chosen=None,
                       Z1=-5., Z2=0., max_iters=None, verbose=False):

If Z1 <= Z2 then points where Ztopo <= Z1 are first selected and then a marching algorithm is used to select neighboring points that satisfy Ztopo <= Z2.

Think of chosen points as "wet" and unchosen points as "dry" and new points are flooded if at least one neighbor is "wet" and the topography is low enough. Starting from deep water (e.g. Z1 = -5) this allows flooding up to MHW (Z2 = 0) without going past dikes that protect dry land with lower elevation.

However, this can also be called with Z1 > Z2, in which case points where Ztopo >= Z1 are first selected and then the marching algorithm is used to select neighboring points that satisfy Ztopo >= Z2. This is useful to select offshore points where the water is shallow, (and that are connected to the shore by a path of shallow points) e.g. to identify points on the continental shelf to set up a flag region for refinement.

If max_iters=None then iterate to convergence. Otherwise, at most max_iters iterations are taken, so setting this to a small value only selects points within this many grid points of where Ztopo <= Z1, useful for buffering or selecting a few onshore points along the coast regardless of elevation.

If prev_pts_chosen is None we are starting from scratch, otherwise we possibly add additional chosen points to an existing array. Points where prev_pts_chosen[i,j]==1 won't change but those ==0 may be changed to 1 based on the new criteria.

If mask==None or mask==False then the entire array is subject to selection. If mask is an array with mask[i,j]==True at some points, then either:

  • these masked points are marked as not selected (pts_chosen[i,j]=0) if prev_pts_chosen==None, or
  • these masked points are not touched if prev_pts_chosen is an array (so previous 0/1 values are preserved).

These arguments are described in more detail below with examples of how they might be used.

output array

The function returns an array pt_chosen with the same shape as Ztopo and has the value 1 at points chosen and 0 at points not chosen.

creating a masked array

This array can be used to define a masked array from Ztopo that masks out the points not chosen via:

from numpy import ma # masked array module
Zmasked = ma.masked_array(Ztopo, mask=logical_not(pt_chosen))

This could be used to plot only the points chosen using the matplotlib function pcolormesh, for example, or the function we have defined in plottools.pcolorcells that better plots finite volume grid cell data with proper alignment and boundary cells.

topofile mask for initializing dry points

The array pt_chosen can be used to create a file that is read in by GeoClaw to identify points where Ztopo is below MHW but where there is dry land because of protection by dikes or levies. This is done by defining a geoclaw.topotools.Topography object and setting its Z attribute based on pt_chosen, and then writing this as a topofile with topo_type==3. Then in setrun.py this file can be specified as a mask that is read in and used when initializing grid cells. There are some subtleties in how this is done, described in more detail in The documentation page Force Cells to be Dry Initially.

fgmax points

The chosen points might be used as fgmax points, as described below.

flag regions

The output array could also be used to define an AMR flag region as a ruled rectangle, using the function region_tools.ruledrectangle_covering_selected_points described in documentation pages Ruled Rectangles and Specifying flagregions for adaptive refinement. This would give a minimal ruled rectangle covering all the chosen points. An example is given below.

The mask argument

mask can be False or None, or else must be an array of the same shape as Ztopo.

The Ztopo array input must be a rectangular array, but sometimes we want to select points covering only a subset of these points, e.g. when defining fgmax points along some stretch of coastline. In this case mask can be used to mask out values we do not want to select. Set mask[i,j] = True at points that should not be chosen.

To mask out points that lie outside some ruled rectangle that has been defined as rr, you can use theRuled RectanglesRuledRectangles.ipynb.

The previous_pts_chosen argument

This argument is useful if you want to apply a sequence of different criteria to choose points. For example, suppose you first want to choose all grid points within 5 points of the coast (as can be done using max_iters) and then supplement this will all grid points below a specified elevation that are farther inland from the coast.

Examples are given below, also of how the mask array works in conjunction with previous_pts_chosen.

Examples

First import some needed modules and set up color maps.

In [1]:
%matplotlib inline
In [2]:
from pylab import *
import os,sys
from numpy import ma # masked arrays
from clawpack.visclaw import colormaps, plottools
from clawpack.geoclaw import topotools, marching_front
from clawpack.amrclaw import region_tools
In [3]:
zmin = -60.
zmax = 40.

land_cmap = colormaps.make_colormap({ 0.0:[0.1,0.4,0.0],
                                     0.25:[0.0,1.0,0.0],
                                      0.5:[0.8,1.0,0.5],
                                      1.0:[0.8,0.5,0.2]})

sea_cmap = colormaps.make_colormap({ 0.0:[0,0,1], 1.:[.8,.8,1]})

cmap, norm = colormaps.add_colormaps((land_cmap, sea_cmap),
                                     data_limits=(zmin,zmax),
                                     data_break=0.)
                                     
sea_cmap_dry = colormaps.make_colormap({ 0.0:[1.0,0.7,0.7], 1.:[1.0,0.7,0.7]})
cmap_dry, norm_dry = colormaps.add_colormaps((land_cmap, sea_cmap_dry),
                                     data_limits=(zmin,zmax),
                                     data_break=0.)

Sample topography from a 1/3 arcsecond DEM

We consider a small region on the SW coast of Whidbey Island north of Maxwelton Beach as an example:

In [4]:
Whidbey1_png = imread('http://www.clawpack.org/gallery/_static/figures/Whidbey1.png')

extent = [-122.46, -122.38, 47.93, 47.96]

figure(figsize=(12,6))
imshow(Whidbey1_png, extent=extent)
gca().set_aspect(1./cos(48*pi/180.))

We read this small portion of the 1/3 arcsecond Puget Sound DEM, available from the NCEI thredds server:

In [5]:
path = 'https://www.ngdc.noaa.gov/thredds/dodsC/regional/puget_sound_13_mhw_2014.nc'
topo = topotools.read_netcdf(path, extent=extent)
In [6]:
figure(figsize=(12,6))
plottools.pcolorcells(topo.X, topo.Y, topo.Z, cmap=cmap, norm=norm)
colorbar(extend='both')
gca().set_aspect(1./cos(48*pi/180.))

This plot shows that there is a region with elevation below MHW (0 in the DEM) where the Google Earth image shows wetland that should not be initialized as a lake. We will tackle this problem below, but first we show how to choose only the DEM points that are close to the shore and/or below a given elevation.

Finding all points below a given elevation

First we choose all points with elevation below 15 m and that are connected to the coast over this topography extent. Note the latter requirement will eliminate the low lying area at the bottom of the figure above near longitude -122.4 (which is connected to the coast through Cultus Bay, but not by points in this extent).

In [7]:
pts_chosen = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=15., max_iters=None) 
Selecting points with Z1 = 0, Z2 = 15, max_iters=279936
Done after 183 iterations with 89871 points chosen
In [8]:
Zmasked = ma.masked_array(topo.Z, logical_not(pts_chosen))
In [9]:
figure(figsize=(12,6))
plottools.pcolorcells(topo.X, topo.Y, Zmasked, cmap=cmap, norm=norm)
colorbar(extend='both')
gca().set_aspect(1./cos(48*pi/180.))

Create a buffer zone along shore

To select more points along the shore where the topography is steep, we could have first used max_iters.

To illustrate this, we start again and fist use max_iters = 20 so that at least 20 grid points are selected near the coast, also setting Z2 = 1e6 (a huge value) so that the arbitrarily high regions will be included if they are within 20 points of the coast:

In [10]:
pts_chosen = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=1e6, max_iters=20) 
Selecting points with Z1 = 0, Z2 = 1e+06, max_iters=20
Done after 20 iterations with 84800 points chosen

Plot what we have so far:

In [11]:
Zmasked = ma.masked_array(topo.Z, logical_not(pts_chosen))

figure(figsize=(12,6))
plottools.pcolorcells(topo.X, topo.Y, Zmasked, cmap=cmap, norm=norm)
colorbar(extend='both')
gca().set_aspect(1./cos(48*pi/180.))