.. _marching_front: Marching Front algorithm ======================== **New in Version 5.7.0.** Adapted from `this notebook `_. The module `clawpack.geoclaw.marching_front` 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 :math:`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 :math:`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 :math:`Z\geq0` or as *wet* if :math:`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*. The use of such a `force_dry_array` in GeoClaw is explained in :ref:`force_dry`. Other applications are also described below along with some examples. Contents -------- - `Function arguments <#mf-args>`__ - `Output array <#mf-output>`__ - `The mask argument <#mf-mask>`__ - `The previous_pts_chosen argument <#mf-prev>`__ - `Examples <#mf-examples>`__ - `Finding points below a given elevation <#mf-find-Zlow>`__ - `Choose points only near shore <#mf-find-nearshore>`__ - `Write out the masked array indicating fgmax points <#mf-fgmax-file>`__ - `Creating an AMR flag region <#mf-amr-flag>`__ - `Determining dry areas below MHW <#mf-dry>`__ 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. .. _mf-output: 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 now defined in `geoclaw.visclaw.plottools.pcolorcells` that better plots finite volume grid cell data with proper alignment and boundary cells. See :ref:`pcolorcells`. 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 :ref:`force_dry`. 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 :ref:`ruled_rectangles`. This would give a minimal ruled rectangle covering all the chosen points. An example is given `below <#amr-flag>`__. .. _mf-mask: 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 the `rr.mask_outside` function to define the mask. See :ref:`ruled_rectangles`. .. _mf-prev: 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`. .. _mf-examples: Examples -------- First import some needed modules and set up color maps. .. code:: ipython3 import os,sys from numpy import ma # masked arrays from clawpack.visclaw import colormaps from clawpack.geoclaw import topotools, marching_front from clawpack.amrclaw import region_tools from clawpack.visclaw.plottools import pcolorcells zmin = -60. zmax = 40. cmap_land = 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]}) cmap_sea = colormaps.make_colormap({ 0.0:[0,0,1], 1.:[.8,.8,1]}) cmap_topo, norm_topo = colormaps.add_colormaps((cmap_land, cmap_sea), data_limits=(zmin,zmax), data_break=0.) cmap_sea_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((cmap_land, cmap_sea_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: .. code:: ipython3 region1_png = imread('region1.png') extent = [-122.46, -122.38, 47.93, 47.96] figure(figsize=(12,6)) imshow(region1_png, extent=extent) gca().set_aspect(1./cos(48*pi/180.)) .. image:: MarchingFront/output_13_0.png We read this small portion of the 1/3 arcsecond Puget Sound DEM, available from the NCEI thredds server: .. code:: ipython3 path = 'https://www.ngdc.noaa.gov/thredds/dodsC/regional/puget_sound_13_mhw_2014.nc' topo = topotools.read_netcdf(path, extent=extent) .. code:: ipython3 figure(figsize=(12,6)) pcolorcells(topo.X, topo.Y, topo.Z, cmap=cmap, norm=norm) colorbar(extend='both') gca().set_aspect(1./cos(48*pi/180.)) .. image:: MarchingFront/output_16_0.png 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. This problem is discussed in :ref:`force_dry`. Here we show how to choose only the DEM points that are close to the shore and/or below a given elevation. .. _mf-find-Zlow: 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). .. code:: ipython3 pts_chosen = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=15., max_iters=None) .. parsed-literal:: Selecting points with Z1 = 0, Z2 = 15, max_iters=279936 Done after 183 iterations with 89871 points chosen .. code:: ipython3 Zmasked = ma.masked_array(topo.Z, logical_not(pts_chosen)) figure(figsize=(12,6)) pcolorcells(topo.X, topo.Y, Zmasked, cmap=cmap, norm=norm) colorbar(extend='both') gca().set_aspect(1./cos(48*pi/180.)) .. image:: MarchingFront/output_22_0.png 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: .. code:: ipython3 pts_chosen = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=1e6, max_iters=20) .. parsed-literal:: 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: .. code:: ipython3 Zmasked = ma.masked_array(topo.Z, logical_not(pts_chosen)) figure(figsize=(12,6)) pcolorcells(topo.X, topo.Y, Zmasked, cmap=cmap, norm=norm) colorbar(extend='both') gca().set_aspect(1./cos(48*pi/180.)) .. image:: MarchingFront/output_27_0.png Then we augment the points already chosen with any points below 15 m and connected to the coast: .. code:: ipython3 pts_chosen = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=15., prev_pts_chosen=pts_chosen, max_iters=None) .. parsed-literal:: Selecting points with Z1 = 0, Z2 = 15, max_iters=279936 Done after 163 iterations with 94297 points chosen .. code:: ipython3 Zmasked = ma.masked_array(topo.Z, logical_not(pts_chosen)) figure(figsize=(12,6)) pcolorcells(topo.X, topo.Y, Zmasked, cmap=cmap, norm=norm) colorbar(extend='both') gca().set_aspect(1./cos(48*pi/180.)) .. image:: MarchingFront/output_30_0.png .. _mf-find-nearshore: Choose points only near shore ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ We can set `Z1 = 0` and `Z2 = -15` to select points that have elevation *greater than* -15 m and are connected to the coast: .. code:: ipython3 pts_chosen_shallow = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=-15., max_iters=None) .. parsed-literal:: Selecting points with Z1 = 0, Z2 = -15, max_iters=279936 Done after 177 iterations with 249577 points chosen .. code:: ipython3 Zshallow = ma.masked_array(topo.Z, logical_not(pts_chosen_shallow)) figure(figsize=(12,6)) pcolorcells(topo.X, topo.Y, Zshallow, cmap=cmap, norm=norm) colorbar(extend='both') gca().set_aspect(1./cos(48*pi/180.)) .. image:: MarchingFront/output_34_0.png Note that this chooses *all* onshore points in addition to offshore points with elevation greater than -15 m. We can take the intersection of this set of points with the onshore points previously chosen to get only the points that lie near the coast: .. code:: ipython3 pts_chosen_nearshore = logical_and(pts_chosen, pts_chosen_shallow) Znearshore = ma.masked_array(topo.Z, logical_not(pts_chosen_nearshore)) figure(figsize=(12,6)) pcolorcells(topo.X, topo.Y, Znearshore, cmap=cmap, norm=norm) colorbar(extend='both') gca().set_aspect(1./cos(48*pi/180.)) .. image:: MarchingFront/output_37_0.png .. _mf-fgmax-file: Write out the masked array indicating fgmax points -------------------------------------------------- One we have selected the desired fgmax points, these can be output in the style of a `topo_type=3` topography file, with a header followed point values at all points on a uniform grid. The values are simply the integer 1 for points that should be used as fgmax points and 0 for other points. Note that format `%1i` is used for compactness. .. code:: ipython3 fname_fgmax_mask = 'fgmax_pts_topostyle.data' topo_fgmax_mask = topotools.Topography() topo_fgmax_mask._x = topo.x topo_fgmax_mask._y = topo.y topo_fgmax_mask._Z = where(pts_chosen_nearshore, 1, 0) # change boolean to 1/0 topo_fgmax_mask.generate_2d_coordinates() topo_fgmax_mask.write(fname_fgmax_mask, topo_type=3, Z_format='%1i') print('Created %s' % fname_fgmax_mask) .. parsed-literal:: Created fgmax_pts_topostyle.data This file `fgmax_pts_topostyle.data` can then be read into GeoClaw, using the new capability of specifying fgmax grids in this manner, using `point_style == 4` as described in :ref:`fgmax`. .. _mf-amr-flag: Creating an AMR flag region --------------------------- Once a set of points has been selected, it can be used to define a ruled rectangle that might be used as an adaptive mesh refinement flag region, for example. See :ref:`ruled_rectangles` and :ref:`flagregions`. You might want to try both `ixy = 'x'` and `ixy = 'y'` in creating the ruled rectangle to see which one covers fewer non-selected points. In this case `ixy = 'x'` is better: .. code:: ipython3 figure(figsize=(12,5)) subplot(121) rr = region_tools.ruledrectangle_covering_selected_points(topo.X, topo.Y, pts_chosen_nearshore, ixy='x', method=0, padding=0, verbose=True) xv,yv = rr.vertices() pcolorcells(topo.X, topo.Y, Znearshore, cmap=cmap, norm=norm) axis([-122.47, -122.40, 47.925, 47.965]) gca().set_aspect(1./cos(48*pi/180.)) plot(xv, yv, 'r') title("With ixy = 'x'") subplot(122) rr = region_tools.ruledrectangle_covering_selected_points(topo.X, topo.Y, pts_chosen_nearshore, ixy='y', method=0, padding=0, verbose=True) xv,yv = rr.vertices() pcolorcells(topo.X, topo.Y, Znearshore, cmap=cmap, norm=norm) axis([-122.47, -122.40, 47.925, 47.965]) gca().set_aspect(1./cos(48*pi/180.)) plot(xv, yv, 'r') title("With ixy = 'y'") .. parsed-literal:: Extending rectangles to cover grid cells RuledRectangle covers 69788 grid points Extending rectangles to cover grid cells RuledRectangle covers 76005 grid points .. image:: MarchingFront/output_42_2.png .. _mf-dry: Determining dry areas below MHW ------------------------------- Note the blue region in the above plot that is inland from the coast and behind a green barrier. Examining Google Earth shows that this is wetland area that probably should not be initialized as a lake. We can identify this by using the marching front algorithm to start at `Z1 = -5` m and fill in up to `Z2 = 0` (MHW). .. code:: ipython3 wet_points = marching_front.select_by_flooding(topo.Z, Z1=-5., Z2=0., max_iters=None) .. parsed-literal:: Selecting points with Z1 = -5, Z2 = 0, max_iters=279936 Done after 112 iterations with 59775 points chosen .. code:: ipython3 Zdry = ma.masked_array(topo.Z, wet_points) figure(figsize=(12,6)) pcolorcells(topo.X, topo.Y, Zdry, cmap=cmap, norm=norm) colorbar(extend='both') gca().set_aspect(1./cos(48*pi/180.)) title('Dry land'); .. image:: MarchingFront/output_46_1.png Now the blue region above is properly identified as being dry land. See :ref:`force_dry` for more discussion of this example, showing how to create an array for GeoClaw in order to indicate that this region should be initialized as dry in spite of being below MHW.