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.

In [1]:
%matplotlib inline
In [2]:
from pylab import *
In [3]:
import os,sys
from imp import reload
from clawpack.geoclaw import topotools, dtopotools, kmltools
from clawpack.visclaw import colormaps
from IPython.display import Image
In [4]:
sys.path.insert(0,'../new_python')
In [5]:
import fgmax_tools, fgmax_routines 
import region_tools, marching_front
In [6]:
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 [7]:
extent = [-132, -122, 44, 52]
topo = topotools.read_netcdf('etopo1', extent=extent, coarsen=4)
In [8]:
topo.plot()
Out[8]:
<matplotlib.axes._subplots.AxesSubplot at 0x121cfe908>

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 [9]:
# 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')
Out[9]:
[<matplotlib.lines.Line2D at 0x1222f7a58>]

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 [10]:
# 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 6350 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: 2458

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

In [11]:
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 2558 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 [12]:
marching_front.plot_masked_Z(topo, nearshore_pts, zmin=-3000, zmax=1000)

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

And the ruled rectangle in the bigger context:

In [13]:
figure(figsize=(10,10))
ax = axes()
topo.plot(axes=ax)
x2,y2 = rr2.vertices()
plot(x2,y2,'r')
ylim(45,52);

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

In [14]:
rr_name = 'RuledRectangle_Coast_46_51'
rr2.write(rr_name + '.data')
rr2.make_kml(fname=rr_name+'.kml', name=rr_name)

The above cell also makes a kml file RuledRectangle_Coast_46_51.kml that can be opened on Google Earth to show this region. Here's a screenshot:

In [15]:
Image('figs/RuledRectangle_Coast_46_51_GE.png', width=500)
Out[15]:

Back to Index.html