Notebook to make an array of fgmax points based on a RuledRectangle and parameters set in params.py
, and a RuledRectangle covering these points for AMR flagging.
Also make the allow_initially_wet
array for this area, if needed. This part of the notebook may need to be customized.
In the next cell use %matplotlib notebook
for live interactive plots or %matplotlib inline
for non-interactive plots.
%matplotlib inline
from pylab import *
import os,sys
from imp import reload
from clawpack.geoclaw import topotools, dtopotools
from clawpack.visclaw import colormaps
loc = 'SE_Whidbey' ## SET THIS to name to be used for this study region
save_figs = True # make png files for figures?
try:
rootdir = os.environ['WA_EMD_2019']
except:
print('*** Need to set environment variable WA_EMD_2019')
print('Set rootdir to ', rootdir)
# directory for topo files:
topodir = rootdir + '/topo'
# directories where some python modules can be found:
sys.path.insert(0, rootdir + '/new_code')
sys.path.insert(0, rootdir + '/topo')
import fgmax_tools, fgmax_routines, kmltools # from new_code
import region_tools, marching_front # from topo
#reload(marching_front) # reload if debugging module
def savefigp(fname):
global save_figs
if save_figs:
savefig(fname, bbox_inches='tight')
print('Created ', fname)
else:
print('save_figs = False')
GEdir = rootdir + '/info/fgmax_regions'
try:
GEfile = GEdir + '/%s_bbox_GE.png' % loc
GEmap = imread(GEfile)
GEextent_file = GEdir + '/%s_bbox_GE_extent.txt' % loc
GEextent_tokens = open(GEextent_file,'r').read().strip().split(',')
GEextent = [float(s) for s in GEextent_tokens]
print('GE extent: ', GEextent)
except:
print('Problem reading Google Earth image file or extent')
GEmap = None
if GEmap is not None:
figure(figsize=(8,8))
imshow(GEmap, extent=GEextent)
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
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.)
PTPSm = topotools.read_netcdf(topodir + '/topofiles/PTPS_merged_13s.nc')
We will select fgmax points that lie within this polygon, with other conditions also imposed later.
fg_poly = region_tools.RuledRectangle()
slu = \
array([[ 47.890, -122.415, -122.303],
[ 48.020, -122.415, -122.303],
[ 48.050, -122.450, -122.345],
[ 48.12 , -122.450, -122.435]])
fg_poly.s = slu[:,0]
fg_poly.lower = slu[:,1]
fg_poly.upper = slu[:,2]
fg_poly.ixy = 2
fg_poly.method = 1
# Cropped topo file:
#extent = [-122.54, -122.41, 47.9, 48.04] # custom extent if desired
x1,x2,y1,y2 = fg_poly.bounding_box()
padding = 5./(3*3600.) # add a few grid cells around bounding box
extent = [x1-padding, x2+padding, y1-padding, y2+padding]
print('Cropping PTPSm topo to: ', extent)
topo = PTPSm.crop(extent)
extent = topo.extent
print('New extent of topo: ', extent)
fg_poly.write('fg_poly_%s.slu' % loc)
def plot_topo(Z, show_cb=True):
#figure(figsize=(8,8))
pc = pcolormesh(topo.X, topo.Y, Z, cmap=cmap, norm=norm)
if show_cb:
cb = colorbar(pc)
cb.set_label('meters')
axis(extent)
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
figure(figsize=(8,8))
plot_topo(topo.Z)
xpoly,ypoly = fg_poly.vertices()
plot(xpoly,ypoly,'w')
title('Topography and fg_poly')
savefigp('fg_poly_%s.png' % loc)
Make a kml file showing the polygon used above, and one showing the bounding box as a rectangle. The latter is useful for capturing a Google Earth image of the fgmax region in order to plot or topography on top of the image.
kmltools.poly2kml((xpoly,ypoly), name='%s_fgmax_poly' % loc, fname='%s_fgmax_poly.kml' % loc)
kmltools.box2kml(extent, name='%s_fgmax_bbox' % loc, fname='%s_fgmax_bbox.kml' % loc)
First use the criterion that points are selected everywhere in the bounding box rectangle, as specified by the extent of the topo
array, and if one of these conditions holds:
Z < 0
,onshore_buffer
grid points of the shoreline, orB < onshore_Z
, and is connected to the shore by a chain of grid points also satisfying this condition.The latter requirement avoids selecting points in a local depression that are separated from the coast by high ground that won't be overtopped.
This is accomplished by a marching front algorithm defined in rootdir + '/topo/marching_front.py'
.
# mask the points not in fg_poly:
mask_out = fg_poly.mask_outside(topo.X, topo.Y)
#mask_fg = logical_or(mask_out, mask_near)
fgmax_pts_chosen = marching_front.select_by_flooding(topo.Z, mask=mask_out,
prev_pts_chosen=None,
Z1=0, Z2=1e6, max_iters=10)
fgmax_pts_chosen = marching_front.select_by_flooding(topo.Z, mask=mask_out,
prev_pts_chosen=fgmax_pts_chosen,
Z1=0, Z2=15., max_iters=None)
# mask the points not selected above:
mask_fg = logical_not(fgmax_pts_chosen)
# topography as masked array:
Z_fg = ma.masked_array(topo.Z, mask_fg)
print('Selected %i fgmax points' % Z_fg.count())
figure(figsize=(8,8))
plot_topo(Z_fg)
plot(xpoly,ypoly,'r')
title('%i fgmax points in specified polygon' % Z_fg.count())
savefigp('fgmax_points_%s.png' % loc)
This is written 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.
fname_fgmax_mask = 'fgmax_pts_topostyle_%s.data' % loc
topo_fgmax_mask = topotools.Topography()
topo_fgmax_mask._x = topo.x
topo_fgmax_mask._y = topo.y
topo_fgmax_mask._Z = where(Z_fg.mask, 0, 1)
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)
Don't need this since we'll put topo in netCDF file.
if 0:
fname_fgmax_topo = 'fgmax_topo_%s.nc' % loc
topo.write(fname_fgmax_topo, topo_type=4)
print('Created %s' % fname_fgmax_topo)
We wish to force GeoClaw to refine the region around the fgmax points to the finest level (generally 1/3" to agree with the resolution of the fgmax points). Rather than trying to specify rectangular regions that cover the fgmax points efficiently, we can use a RuledRectangle
.
We can make a RuledRectangle
that includes all the fgmax points in two different ways, either with s = x
so the rules go from ylower
to yupper
on each column of grid points, or with s = y
so that the rules go from xlower
to xupper
along each row of grid points.
Which approach is best depends on the geometry of the fgmax points. We would like the RuledRectangle
to cover as few points as possible that are not fgmax points.
Here we compute and display both options.
x = topo.x
y = topo.y
# Ruled rectangle with s = x:
rr1 = region_tools.RuledRectangle()
rr1.ixy = 1
s = []
lower = []
upper = []
for i in range(len(x)):
if fgmax_pts_chosen[:,i].sum() > 0:
j = where(fgmax_pts_chosen[:,i]==1)[0]
j1 = j.min()
j2 = j.max()
s.append(x[i])
lower.append(y[j1])
upper.append(y[j2])
rr1.s = array(s)
rr1.lower = array(lower)
rr1.upper = array(upper)
rr1.method = 0
rr1.ds = s[1] - s[0]
rr1_npts = int(ceil(sum(rr1.upper - rr1.lower) *3*3600))
print('RuledRectangle rr1 covers %s grid points' % rr1_npts)
# Ruled rectangle with s = y:
rr2 = region_tools.RuledRectangle()
rr2.ixy = 2
s = []
lower = []
upper = []
for j in range(len(y)):
if fgmax_pts_chosen[j,:].sum() > 0:
i = where(fgmax_pts_chosen[j,:]==1)[0]
i1 = i.min()
i2 = i.max()
s.append(y[j])
lower.append(x[i1])
upper.append(x[i2])
rr2.s = array(s)
rr2.lower = array(lower)
rr2.upper = array(upper)
rr2.method = 0
rr2.ds = s[1] - s[0]
rr2_npts = int(ceil(sum(rr2.upper - rr2.lower)*3*3600))
print('RuledRectangle rr2 covers %s grid points' % rr2_npts)
# Create files for AMR flag regions:
fname1 = 'RuledRectangle_%s_ixy1.data' % loc
rr1.write(fname1)
print('Created ',fname1)
fname2 = 'RuledRectangle_%s_ixy2.data' % loc
rr2.write(fname2)
print('Created ',fname2)
# Plot comparison of two approaches:
figure(figsize=(10,8))
subplot(121)
plot_topo(Z_fg, show_cb=False)
x2,y2 = rr1.vertices()
plot(x2,y2,'r')
title(fname1 + '\n%i points covered' % rr1_npts)
subplot(122)
plot_topo(Z_fg, show_cb=False)
x2,y2 = rr2.vertices()
plot(x2,y2,'r')
title(fname2 + '\n%i points covered' % rr2_npts)
savefigp('RuledRectangles_%s.png' % loc)
If this region does not have any land below MHW, can simply set:
allow_wet_init = None
figure(figsize=(10,8))
subplot(121)
#pc = pcolormesh(topo.X, topo.Y, Z_dry, cmap=cmap_dry, norm=norm_dry)
pc = pcolormesh(topo.X, topo.Y, Z_fg, cmap=cmap, norm=norm)
#cb = colorbar(pc)
#cb.set_label('meters')
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('%s: %i fgmax points\n(no dry regions below MHW)' % (loc, Z_fg.count()))
subplot(122)
if GEmap is not None:
imshow(GEmap, extent=GEextent)
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
else:
print('Could not make plot with image')
savefigp('allow_wet_init_%s.png' % loc)
Z_dry = ma.masked_where(logical_or(topo.Z<0, mask_fg), topo.Z)
Z_dry_above0 = Z_dry
fig,ax,png_extent,kml_dpi = kmltools.pcolorcells_for_kml(topo.x, topo.y, Z_dry_above0,
png_filename='dry_above0_for_kml.png',
dpc=2, cmap=cmap_dry, norm=norm_dry)
mask_dry = logical_not(Z_dry.mask)
mm = logical_or(mask_fg, mask_dry)
Z_wet = ma.masked_array(topo.Z, mm)
fig,ax,png_extent,kml_dpi = kmltools.pcolorcells_for_kml(topo.x, topo.y, Z_wet,
png_filename='wet_for_kml.png',
dpc=2, cmap=cmap, norm=norm)
kmltools.png2kml(topo.extent, png_files=['dry_above0_for_kml.png',
'wet_for_kml.png'],
png_names=['dry above MHW','allow wet'],
name='fgmax_points_%s' % loc,
fname='fgmax_points_%s.kmz' % loc)
print('original shape of X,Y,Z: ', Z_fg.shape)
iindex = array([i for i in range(Z_fg.shape[0]) \
if ma.notmasked_edges(Z_fg[i,:]) is not None])
jindex = array([j for j in range(Z_fg.shape[1]) \
if ma.notmasked_edges(Z_fg[:,j]) is not None])
i1 = iindex.min(); i2 = iindex.max()
j1 = jindex.min(); j2 = jindex.max()
print('First and last iindex: %7i, %7i' % (i1,i2))
print('First and last jindex: %7i, %7i' % (j1,j2))
Zcrop = Z_fg[i1:i2+1, j1:j2+1]
print('new shape of X,Y,Z: ', Zcrop.shape)
Xcrop = topo.X[i1:i2+1, j1:j2+1]
Ycrop = topo.Y[i1:i2+1, j1:j2+1]
figure(figsize=(8,8))
pc = pcolormesh(Xcrop, Ycrop, Zcrop, cmap=cmap, norm=norm)
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
#axis(extent)
#plot(xpoly,ypoly,'r');
fgm = fgmax_tools.FGmaxMaskedGrid()
fgm.id = loc
fgm.X = Xcrop
fgm.Y = Ycrop
fgm.Z = Zcrop
fgm.fgmax_point = where(Zcrop.mask, 0, 1)
print('Created fgm object')
print('Check that this has the expected number of fgmax points: ',\
fgm.fgmax_point.sum())
if allow_wet_init is not None:
fgm.allow_wet_init = allow_wet_init_topo.Z[i1:i2+1, j1:j2+1]
fname_nc = '%s_input.nc' % loc
fgmax_routines.make_nc_input(fname_nc, fgm, force=True)
print('Created ',fname_nc)
Save the *.data, *.nc
(and optionally *.png
and *.kml
) files created by this notebook to a (new) directory input_files
if you don't already have them saved in this directory, since this is where GeoClaw will look for them. Back up an old version first if you want to, for safety.