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.

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.

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.

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.

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.

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

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.

`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.

`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`

.

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.)
```

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.

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)
```

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.))
```

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)
```

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.))
```