Make Flagregions along the Coast

This example shows how to make a Ruled Rectangle for AMR flagging along the coast, based on first selecting points within some region where the water depth satisfies a costraint. This is done by applyting the marching front algorithm to a topography DEM at a suitable resolution.

Here we make a ruled rectangle that covers the continental shelf from the Columbia River up to Vancouver Island, a region where we might want to require finer AMR grids for tsunami modeling in Washington State in order to capture edge waves that are trapped on the shelf.

In [1]:
%matplotlib inline
In [2]:
from pylab import *
In [3]:
import os,sys
from imp import reload
from clawpack.amrclaw import region_tools
from clawpack.geoclaw import topotools, dtopotools, kmltools, fgmax_tools, marching_front
from clawpack.visclaw import colormaps
from IPython.display import Image
In [4]:
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.8,0.8], 1.:[1.0,0.8,0.8]})
cmap_dry, norm_dry = colormaps.add_colormaps((land_cmap, sea_cmap_dry),
                                     data_limits=(zmin,zmax),
                                     data_break=0.)

Download etopo1 data coarsened to 4 minute resolution

Even if the computational grid will be finer than 4 arcminutes, this is sufficient resolution for creating a ruled rectangle to use as a flagregion.

In [5]:
extent = [-132, -122, 44, 52]
topo = topotools.read_netcdf('etopo1', extent=extent, coarsen=4)
In [6]:
figure(figsize=(13,10))
ax = axes()
topo.plot(axes=ax)
title('etopo1 coarsend to 4-minute');

Restrict region of interest:

We first define a simple ruled region with the region of interest. In this cae we are interested in making a flagregion that goes along the continental shelf from the Columbia River at about 46N to the northern tip of Vancouver Island at about 51N, and extending only slightly into the Strait of Juan de Fuca (since we used other flagregions in the Strait and Puget Sound).

In [7]:
# Specify a RuledRectangle the our flagregion should lie in:
rrect = region_tools.RuledRectangle()
rrect.ixy = 'y'  # s = latitude
rrect.s = array([46,48,48.8,50,51.])
rrect.lower = -131*ones(rrect.s.shape)
rrect.upper = array([-123.,-124,-124,-126,-128])
rrect.method = 1

xr,yr = rrect.vertices()

figure(figsize=(10,10))
ax = axes()
topo.plot(axes=ax)
plot(xr,yr,'r');

We could use the red ruled rectangle plotted above directly as our flagregion, but instead we will trim this down to only cover the continental shelf a shore region:

In [8]:
# Start with a mask defined by the ruled rectangle `rrect` defined above:
mask_out = rrect.mask_outside(topo.X, topo.Y)

# select onshore points within 2 grip points of shore:
pts_chosen_Zabove0 = marching_front.select_by_flooding(topo.Z, mask=mask_out, 
                                                       prev_pts_chosen=None, 
                                                       Z1=0, Z2=1e6, max_iters=2)
# select offshore points down to 1000 m depth:
pts_chosen_Zbelow0 = marching_front.select_by_flooding(topo.Z, mask=None, 
                                                       prev_pts_chosen=None, 
                                                       Z1=0, Z2=-1000., max_iters=None)
# buffer offshore points with another 10 grid cells:
pts_chosen_Zbelow0 = marching_front.select_by_flooding(topo.Z, mask=None, 
                                                       prev_pts_chosen=pts_chosen_Zbelow0, 
                                                       Z1=0, Z2=-5000., max_iters=10)

# Take the intersection of the two sets of points selected above:
nearshore_pts = where(pts_chosen_Zabove0+pts_chosen_Zbelow0 == 2, 1, 0)
print('Number of nearshore points: %i' % nearshore_pts.sum())

Selecting points with Z1 = 0, Z2 = 1e+06, max_iters=2
Done after 2 iterations with 6587 points chosen
Selecting points with Z1 = 0, Z2 = -1000, max_iters=18271
Done after 28 iterations with 8552 points chosen
Selecting points with Z1 = 0, Z2 = -5000, max_iters=10
Done after 10 iterations with 10085 points chosen
Number of nearshore points: 2537

Now create the minimal ruled rectangle rr2 that covers the grid cells selected above:

In [9]:
rr2 = region_tools.ruledrectangle_covering_selected_points(topo.X, topo.Y,
                                                          nearshore_pts, ixy='y', method=0,
                                                          verbose=True)

Extending rectangles to cover grid cells
RuledRectangle covers 2640 grid points

Here's the topography, masked to only show the points selected in nearshore_pts, along with the minimal ruled rectangle that covers these points.

In [10]:
Z = ma.masked_where(logical_not(nearshore_pts), topo.Z)
masked_topo = topotools.Topography()
masked_topo.set_xyZ(topo.x, topo.y, Z)

figure(figsize=(10,10))
ax = axes()
masked_topo.plot(axes=ax)

x2,y2 = rr2.vertices()
plot(x2,y2,'r')
ylim(45,52);

And the ruled rectangle in the bigger context:

In [11]:
figure(figsize=(10,10))
ax = axes()
topo.plot(axes=ax)
x2,y2 = rr2.vertices()
plot(x2,y2,'r')
ylim(45,52)
title('RuledRectangle_Coast_46_51')
savefig('RuledRectangle_Coast_46_51.png')

From this we can make a .data file that can be used as a flagregion spatial_extent_file, see FlagRegions.ipynb for discussion.

In [12]:
rr_name = 'RuledRectangle_Coast_46_51'
rr2.write(rr_name + '.data')

We can also write out the ruled rectangle as a kml file RuledRectangle_Coast_46_51.kml that can be opened on Google Earth to show this region.

In [13]:
rr2.make_kml(fname=rr_name+'.kml', name=rr_name)

Here's a screenshot:

In [14]:
Image('http://www.clawpack.org/gallery/_static/figures/RuledRectangle_Coast_46_51_GE.png', width=400)
Out[14]:

Use as a flagregion in GeoClaw

To use this ruled rectangle as a flagregion, specifying for example that within this region at least 3 levels of AMR patches and at most 5 levels should be used, lines similar to those shown below should be added to setrun.py:

from clawpack.amrclaw.data import FlagRegion

# append as many flagregions as desired to this list:
flagregions = rundata.flagregiondata.flagregions 

flagregion = FlagRegion(num_dim=2)
flagregion.name = 'Region_Coast_46_51'
flagregion.minlevel = 3
flagregion.maxlevel = 5
flagregion.t1 = 0.
flagregion.t2 = 1e9
flagregion.spatial_region_type = 2  # Ruled Rectangle
flagregion.spatial_region_file = 'RuledRectangle_Coast_46_51.data'
flagregions.append(flagregion)

Of course other flagregions can also be specified if there is a coastal location where much higher resolution is required. See the documentation for more details:

In [ ]: