Skip to content

Commit

Permalink
added MakeFlagregionsCoast.ipynb
Browse files Browse the repository at this point in the history
  • Loading branch information
@rjleveque
rjleveque committed May 21, 2020
1 parent f7234b3 commit 1a1d343
Showing 1 changed file with 356 additions and 0 deletions.
356 changes: 356 additions & 0 deletions notebooks/geoclaw/MakeFlagregionsCoast.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,356 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Make Flagregions along the Coast\n",
"\n",
"This example shows how to make a [Ruled Rectangle](http:https://www.clawpack.org/ruled_rectangles.html) for [AMR flagging](http:https://www.clawpack.org/flagregions.html) 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](http:https://www.clawpack.org/marching_front.html) to a topography DEM at a suitable resolution.\n",
"\n",
"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. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pylab import *"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import os,sys\n",
"from imp import reload\n",
"from clawpack.amrclaw import region_tools\n",
"from clawpack.geoclaw import topotools, dtopotools, kmltools, fgmax_tools, marching_front\n",
"from clawpack.visclaw import colormaps\n",
"from IPython.display import Image"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"zmin = -60.\n",
"zmax = 40.\n",
"\n",
"land_cmap = colormaps.make_colormap({ 0.0:[0.1,0.4,0.0],\n",
" 0.25:[0.0,1.0,0.0],\n",
" 0.5:[0.8,1.0,0.5],\n",
" 1.0:[0.8,0.5,0.2]})\n",
"\n",
"sea_cmap = colormaps.make_colormap({ 0.0:[0,0,1], 1.:[.8,.8,1]})\n",
"\n",
"cmap, norm = colormaps.add_colormaps((land_cmap, sea_cmap),\n",
" data_limits=(zmin,zmax),\n",
" data_break=0.)\n",
" \n",
"sea_cmap_dry = colormaps.make_colormap({ 0.0:[1.0,0.8,0.8], 1.:[1.0,0.8,0.8]})\n",
"cmap_dry, norm_dry = colormaps.add_colormaps((land_cmap, sea_cmap_dry),\n",
" data_limits=(zmin,zmax),\n",
" data_break=0.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Download etopo1 data coarsened to 4 minute resolution\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"extent = [-132, -122, 44, 52]\n",
"topo = topotools.read_netcdf('etopo1', extent=extent, coarsen=4)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"figure(figsize=(13,10))\n",
"ax = axes()\n",
"topo.plot(axes=ax)\n",
"title('etopo1 coarsend to 4-minute');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Restrict region of interest:\n",
"\n",
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Specify a RuledRectangle the our flagregion should lie in:\n",
"rrect = region_tools.RuledRectangle()\n",
"rrect.ixy = 'y' # s = latitude\n",
"rrect.s = array([46,48,48.8,50,51.])\n",
"rrect.lower = -131*ones(rrect.s.shape)\n",
"rrect.upper = array([-123.,-124,-124,-126,-128])\n",
"rrect.method = 1\n",
"\n",
"xr,yr = rrect.vertices()\n",
"\n",
"figure(figsize=(10,10))\n",
"ax = axes()\n",
"topo.plot(axes=ax)\n",
"plot(xr,yr,'r');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Start with a mask defined by the ruled rectangle `rrect` defined above:\n",
"mask_out = rrect.mask_outside(topo.X, topo.Y)\n",
"\n",
"# select onshore points within 2 grip points of shore:\n",
"pts_chosen_Zabove0 = marching_front.select_by_flooding(topo.Z, mask=mask_out, \n",
" prev_pts_chosen=None, \n",
" Z1=0, Z2=1e6, max_iters=2)\n",
"# select offshore points down to 1000 m depth:\n",
"pts_chosen_Zbelow0 = marching_front.select_by_flooding(topo.Z, mask=None, \n",
" prev_pts_chosen=None, \n",
" Z1=0, Z2=-1000., max_iters=None)\n",
"# buffer offshore points with another 10 grid cells:\n",
"pts_chosen_Zbelow0 = marching_front.select_by_flooding(topo.Z, mask=None, \n",
" prev_pts_chosen=pts_chosen_Zbelow0, \n",
" Z1=0, Z2=-5000., max_iters=10)\n",
"\n",
"# Take the intersection of the two sets of points selected above:\n",
"nearshore_pts = where(pts_chosen_Zabove0+pts_chosen_Zbelow0 == 2, 1, 0)\n",
"print('Number of nearshore points: %i' % nearshore_pts.sum())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now create the minimal ruled rectangle `rr2` that covers the grid cells selected above:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rr2 = region_tools.ruledrectangle_covering_selected_points(topo.X, topo.Y,\n",
" nearshore_pts, ixy='y', method=0,\n",
" verbose=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the topography, masked to only show the points selected in `nearshore_pts`, along with the minimal ruled rectangle that covers these points."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Z = ma.masked_where(logical_not(nearshore_pts), topo.Z)\n",
"masked_topo = topotools.Topography()\n",
"masked_topo.set_xyZ(topo.x, topo.y, Z)\n",
"\n",
"figure(figsize=(10,10))\n",
"ax = axes()\n",
"masked_topo.plot(axes=ax)\n",
"\n",
"x2,y2 = rr2.vertices()\n",
"plot(x2,y2,'r')\n",
"ylim(45,52);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the ruled rectangle in the bigger context:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"figure(figsize=(10,10))\n",
"ax = axes()\n",
"topo.plot(axes=ax)\n",
"x2,y2 = rr2.vertices()\n",
"plot(x2,y2,'r')\n",
"ylim(45,52)\n",
"title('RuledRectangle_Coast_46_51')\n",
"savefig('RuledRectangle_Coast_46_51.png')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From this we can make a `.data` file that can be used as a flagregion `spatial_extent_file`, see [FlagRegions.ipynb](FlagRegions.ipynb) for discussion."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rr_name = 'RuledRectangle_Coast_46_51'\n",
"rr2.write(rr_name + '.data')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rr2.make_kml(fname=rr_name+'.kml', name=rr_name)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's a screenshot:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Image('http:https://www.clawpack.org/gallery/_static/figures/RuledRectangle_Coast_46_51_GE.png', width=400)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Use as a flagregion in GeoClaw\n",
"\n",
"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`:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" from clawpack.amrclaw.data import FlagRegion\n",
"\n",
" # append as many flagregions as desired to this list:\n",
" flagregions = rundata.flagregiondata.flagregions \n",
"\n",
" flagregion = FlagRegion(num_dim=2)\n",
" flagregion.name = 'Region_Coast_46_51'\n",
" flagregion.minlevel = 3\n",
" flagregion.maxlevel = 5\n",
" flagregion.t1 = 0.\n",
" flagregion.t2 = 1e9\n",
" flagregion.spatial_region_type = 2 # Ruled Rectangle\n",
" flagregion.spatial_region_file = 'RuledRectangle_Coast_46_51.data'\n",
" flagregions.append(flagregion)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:\n",
"\n",
"- [AMR refinement criteria](http:https://www.clawpack.org/refinement.html)\n",
"- [flagregions](http:https://www.clawpack.org/flagregions.html)\n",
"- [making flagregions using ruled rectangles](http:https://www.clawpack.org/marching_front.html#mf-amr-flag)\n",
"- [ruled rectangles](http:https://www.clawpack.org/ruled_rectangles.html)\n",
"- [marching front algorithm](http:https://www.clawpack.org/marching_front.html)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.10"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

0 comments on commit 1a1d343

Please sign in to comment.