make_input_files -- SW_Whidbey

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.

In [1]:
%matplotlib inline
In [2]:
from pylab import *
In [3]:
import os,sys
from imp import reload
from clawpack.geoclaw import topotools, dtopotools
from clawpack.visclaw import colormaps

Set some values needed below:

In [4]:
loc = 'SW_Whidbey'   ## SET THIS to name to be used for this study region

save_figs = True        # make png files for figures?
In [5]:
try:
    rootdir = os.environ['WA_EMD_2019']
except:
    print('*** Need to set environment variable WA_EMD_2019')
print('Set rootdir to ', rootdir)

Set rootdir to  /Users/rjl/git/WA_EMD_2019
In [6]:
# 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')
In [7]:
import fgmax_tools, fgmax_routines, kmltools  # from new_code
import region_tools, marching_front # from topo
#reload(marching_front)  # reload if debugging module
In [8]:
def savefigp(fname):
    global save_figs
    if save_figs:
        savefig(fname, bbox_inches='tight')
        print('Created ', fname)
    else:
        print('save_figs = False')

Load Google Earth image, if available

In [9]:
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);

GE extent:  [-122.531389, -122.409537, 47.89963, 48.040463]
In [10]:
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.)
In [11]:
PTPSm = topotools.read_netcdf(topodir + '/topofiles/PTPS_merged_13s.nc')

Ruled rectangle used as boundary polygon for fgmax points

We will select fgmax points that lie within this polygon, with other conditions also imposed later.

In [12]:
fg_poly = region_tools.RuledRectangle()
fg_poly.s = array([-122.531, -122.46 , -122.41 ])
fg_poly.lower = array([47.95, 47.91, 47.9])
fg_poly.upper = array([48, 48.04, 47.95])
fg_poly.ixy = 1
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)

Cropping PTPSm topo to:  [-122.53146296296296, -122.40953703703704, 47.899537037037035, 48.04046296296296]
New extent of topo:  (-122.53138890705502, -122.40953705861502, 47.89962962123999, 48.040462962949995)
In [13]:
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)
In [14]:
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)

Created  fg_poly_SW_Whidbey.png

Make kml files

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.

In [15]:
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)

Polygon:     -122.531000   47.950000
             -122.460000   47.910000
             -122.410000   47.900000
             -122.410000   47.950000
             -122.460000   48.040000
             -122.531000   48.000000
             -122.531000   47.950000
Created  SW_Whidbey_fgmax_poly.kml
Box:   -122.531389  -122.409537   47.899630   48.040463
Created  SW_Whidbey_fgmax_bbox.kml

Select fgmax points within this polygon

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:

  • the point is offshore at any depth of water, i.e., Z < 0,
  • the point lies within onshore_buffer grid points of the shoreline, or
  • the topography satisfies B < 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'.

In [16]:
# mask the points not in fg_poly:
mask_out = fg_poly.mask_outside(topo.X, topo.Y)
In [17]:
fgmax_pts_chosen = marching_front.select_by_flooding(topo.Z, mask=mask_out, 
                                                     prev_pts_chosen=None, 
                                                     Z1=0, Z2=1e6, max_iters=10)

Selecting points with Z1 = 0, Z2 = 1e+06, max_iters=10
Done after 10 iterations with 707137 points chosen
In [18]:
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) 

Selecting points with Z1 = 0, Z2 = 15, max_iters=2004474
Done after 287 iterations with 758485 points chosen
In [19]:
# 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())

Selected 758485 fgmax points
In [20]:
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)

Created  fgmax_points_SW_Whidbey.png

Write out the masked array indicating fgmax points

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.

In [21]:
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)

Created fgmax_pts_topostyle_SW_Whidbey.data

Write out topo file with DEM Z values

Don't need this since we'll put topo in netCDF file.

In [22]:
if 0:
    fname_fgmax_topo = 'fgmax_topo_%s.nc' % loc
    topo.write(fname_fgmax_topo, topo_type=4)
    print('Created %s' % fname_fgmax_topo)

Make RuledRectangle for AMR flagging to the finest level

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.

In [23]:
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)

RuledRectangle rr1 covers 880728 grid points
RuledRectangle rr2 covers 785267 grid points
Created  RuledRectangle_SW_Whidbey_ixy1.data
Created  RuledRectangle_SW_Whidbey_ixy2.data
In [24]:
# 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)

Created  RuledRectangles_SW_Whidbey.png

Deal with dry land below MHW

If this region does not have any land below MHW, can simply set:

In [25]:
#allow_wet_init = None

Otherwise, determine which points should be intially wet. This is only necessary if there are regions near the coast where the topography value Z is below 0 but the initial water depth should be 0 because the region is behind a dike or levy.

The basic approach is to process the DEM to determine what points are both below MHW and connected to deep water by a path of grid points that is entirely below MHW.

Normally the approach below is sufficient: Start by selecting points with Z < -5 and then iterate to flood all points up to Z = 0 (MHW):

In [26]:
wet_points = marching_front.select_by_flooding(topo.Z, mask=None, 
                                                     prev_pts_chosen=None, 
                                                     Z1=-5., Z2=0., max_iters=None)

Selecting points with Z1 = -5, Z2 = 0, max_iters=2004474
Done after 466 iterations with 971163 points chosen
In [27]:
mask_dry = logical_not(wet_points)
mm = logical_or(mask_fg, wet_points)
Z_dry = ma.masked_array(topo.Z, mm)
mask_dry_onshore = logical_and(mask_dry, topo.Z<0.)
mm = logical_or(mask_fg, mask_dry_onshore)
Z_fg2= ma.masked_array(topo.Z, mm)

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_fg2, cmap=cmap, norm=norm)
#cb = colorbar(pc, shrink=0.5, extend='both')
#cb.set_label('meters')
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);


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)

The plots above shows areas marked below MHW that were determined to be dry as pink. Note that there is a region east of Deer Harbor in Useless Bay (near $(-122.48, 48)$) where the Google Earth image shows dry land (mostly a golf course) while the algorithm above declares this region wet. This is due to a gap in the dike and/or the fact that MHW is not referenced exactly correctly in this location (it varies a lot in Puget Sound, so a single DEM can't get it right everywhere).

To address this, we first flood only up to $Z = -0.5$ meters, which was determined to be low enough that water does not overtop the dike in question. Then we flood everywhere except in this problem region up to $Z=0$ so that GeoClaw initializes properly elsewhere. (In the problem region GeoClaw may initialize some cells near the shoreline as dry even though they are on the wet side of the dike, generating a bit of flow will immediately. Water will also slowly leak through the low point in the dike as well, but this is the best we can do.)

In [28]:
wet_points = marching_front.select_by_flooding(topo.Z, mask=None, 
                                                     prev_pts_chosen=None, 
                                                     Z1=-5., Z2=-0.5, max_iters=None)

# generate mask to avoid changing wet/dry values in the problem region:
x1,x2,y1,y2 = [-122.485, -122.48, 47.985, 48.015]
xy_in_bad_region = logical_and(logical_and(topo.X>=x1, topo.X<=x2),\
                               logical_and(topo.Y>=y1, topo.Y<=y2))
mask_out_bad_region = logical_or(mask_out, xy_in_bad_region)

wet_points = marching_front.select_by_flooding(topo.Z, mask=mask_out_bad_region, 
                                                     prev_pts_chosen=wet_points, 
                                                     Z1=-5., Z2=0., max_iters=None)

Selecting points with Z1 = -5, Z2 = -0.5, max_iters=2004474
Done after 195 iterations with 945754 points chosen
Selecting points with Z1 = -5, Z2 = 0, max_iters=2004474
Done after 58 iterations with 948052 points chosen
In [29]:
mask_dry = logical_not(wet_points)
mm = logical_or(mask_fg, wet_points)
Z_dry = ma.masked_array(topo.Z, mm)
mask_dry_onshore = logical_and(mask_dry, topo.Z<0.)
mm = logical_or(mask_fg, mask_dry_onshore)
Z_fg2= ma.masked_array(topo.Z, mm)

figure(figsize=(8,8))
pc = pcolormesh(topo.X, topo.Y, Z_dry, cmap=cmap_dry, norm=norm_dry)
pc = pcolormesh(topo.X, topo.Y, Z_fg2, 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(dry regions below MHW pink)' % (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')

Now we observe dry land where expected in the plot above.

Make figure for viewing as kml

In [30]:
reload(kmltools)
Z_dry_above0 = ma.masked_where(Z_dry<0., 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)

Using kml_dpi = 203,figure size 12.975369 by 14.995074 inches
Figure has 1317 by 1522 grid cells of uniform color
Dots per cell in x: 2.000000,  in y: 2.000000
       These should be integers
Created  dry_above0_for_kml.png
In [31]:
Z_dry_below0 = ma.masked_where(Z_dry>0., Z_dry)
fig,ax,png_extent,kml_dpi = kmltools.pcolorcells_for_kml(topo.x, topo.y, Z_dry_below0,
                                                 png_filename='dry_below0_for_kml.png',
                                                 dpc=2, cmap=cmap_dry, norm=norm_dry)

Using kml_dpi = 203,figure size 12.975369 by 14.995074 inches
Figure has 1317 by 1522 grid cells of uniform color
Dots per cell in x: 2.000000,  in y: 2.000000
       These should be integers
Created  dry_below0_for_kml.png
In [32]:
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)

Using kml_dpi = 203,figure size 12.975369 by 14.995074 inches
Figure has 1317 by 1522 grid cells of uniform color
Dots per cell in x: 2.000000,  in y: 2.000000
       These should be integers
Created  wet_for_kml.png
In [33]:
kmltools.png2kml(topo.extent, png_files=['dry_above0_for_kml.png', 'dry_below0_for_kml.png',
                                         'wet_for_kml.png'],
                 png_names=['dry above MHW','dry below MHW','allow wet'], 
                 name='fgmax_points_%s' % loc, 
                 fname='fgmax_points_%s.kmz' % loc)

Extent:   -122.531389  -122.409537   47.899630   48.040463
Created  fgmax_points_SW_Whidbey.kmz

Make allow_wet_init array after padding

The wet_points identified above are points that the DEM indicates should be initially wet.

In GeoClaw we don't use the DEM topography directly, instead the B values are cell averages over finite volume cells, that ideally are centered at the DEM points with the same resolution.

In GeoClaw we want to force a point to be dry only if wet_points[i,j] == 0 and also the 8 neighboring points have this property. This picks up all the locations well separated from the coast. We want to avoid forcing a GeoClaw cell to be dry if it is at the coastline, since the B value might be slightly negative even if the DEM Z value is positive. Forcing it to be dry while the neighbor offshore is wet will cause immediate flow near the coast.

We accomplish this by setting a new array allow_wet_init We will allow GeoClaw to initialize as wet or dry depending on the value of B in cells where allow_wet_init[i,j] = 1 and only force cells to be dry if allow_wet_init[i,j] = 0.

Summing wet_points over 3x3 blocks of point centered at each [i,j] allows quickly determining if any of the 8 neighboring points are potentially wet.

In [34]:
wet_points_sum = wet_points[1:-1,1:-1] + wet_points[0:-2,1:-1] + wet_points[2:,1:-1] + \
                 wet_points[1:-1,0:-2] + wet_points[0:-2,0:-2] + wet_points[2:,0:-2] + \
                 wet_points[1:-1,2:] + wet_points[0:-2,2:] + wet_points[2:,2:]
allow_wet_init = ones(wet_points.shape)
allow_wet_init[1:-1,1:-1] = wet_points_sum
allow_wet_init = where(allow_wet_init > 0, 1, 0)  # 1 if point or any neighbor is wet

Confirm that allow_wet_init contains more points near the shore than wet_points:

In [35]:
print('wet_points has %i nonzeros, allow_wet_init has %i nonzeros'\
       % (wet_points.sum(), allow_wet_init.sum()))

wet_points has 948052 nonzeros, allow_wet_init has 955419 nonzeros

In the figure below, wet_points==1 at blue points and red points are dry according to the results of the marching algorithm applied to the DEM. After the padding above, only red points inside the black contour are forced to be dry and other points are initialized wet or dry according to the GeoClaw B value.

In [36]:
figure(figsize=(8,8))
contourf(topo.X, topo.Y, wet_points, [-1,0.5,2], colors=[[1,.5,.5],[.5,.5,1]])
contour(topo.X, topo.Y, allow_wet_init, [0.5], colors='k')
axis([-122.495,-122.475,47.986,47.998])
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
In [37]:
mask_dry = logical_not(allow_wet_init)
Z_wet = ma.masked_array(topo.Z, mask_dry)

# create masked array Z_dry to plot dry fgmax points as pink:
mm = logical_or(mask_fg, logical_not(mask_dry))
Z_dry = ma.masked_array(topo.Z, mm)

# create Z_fg2 with onshore dry points masked
# only for plotting purpose, so pink shows through:
mask_dry_onshore = logical_and(mask_dry, topo.Z<0.)
mm = logical_or(mask_fg, mask_dry_onshore)
Z_fg2= ma.masked_array(topo.Z, mm)

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_fg2, 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(dry regions below MHW pink)' % (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)

Created  allow_wet_init_SW_Whidbey.png

Create the allow_wet_init file used by GeoClaw

First plot the full array to make sure it looks reasonable:

In [38]:
figure(figsize=(10,8))
subplot(121)
contourf(topo.X, topo.Y, allow_wet_init, [-1,0.5,2], colors=[[1,.5,.5],[.5,.5,1]])
#axis([-122.495,-122.475,47.986,47.998])
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);

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')
In [39]:
allow_wet_init_topo = topotools.Topography()
allow_wet_init_topo._x = topo.x
allow_wet_init_topo._y = topo.y     
allow_wet_init_topo._Z = allow_wet_init
allow_wet_init_topo.generate_2d_coordinates()

fname_allow_wet_init = 'allow_wet_init_%s.data' % loc
allow_wet_init_topo.write(fname_allow_wet_init, topo_type=3, Z_format='%1i')
print('Created %s' % fname_allow_wet_init)

Created allow_wet_init_SW_Whidbey.data

Write fgmax file as xyz.data

This format is no longer used, but for previous projects we needed this list of all fgmax points.

In [40]:
if 0:
    fg = fgmax_tools.FGmaxGrid()
    Z1d = ravel(Z_fg)
    X1d = ravel(topo.X)
    Y1d = ravel(topo.Y)
    km, = where(logical_not(Z1d.mask))
    print('Found %i unmasked fgmax points' % len(km))
    fg.X = X1d[km]
    fg.Y = Y1d[km]
    fg.Z = Z1d[km]
    fg.point_style = 0  # points
    fg.npts = len(fg.X)
    fg.xy_fname = 'fgmax_xyz_%s.data' % loc

    xydata = np.vstack([fg.X,fg.Y,fg.Z]).T
    np.savetxt(fg.xy_fname, xydata, header='%8i' % len(fg.X), 
                  comments='', fmt='%24.14e')
    print('Created ', fg.xy_fname)

Create netCDF file with topo, fgmax points, wet mask:

First crop the arrays down to the bounding box of the fgmax points.

In [41]:
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]

original shape of X,Y,Z:  (1522, 1317)
First and last iindex:       5,    1362
First and last jindex:       5,    1311
new shape of X,Y,Z:  (1358, 1307)

Plot to make sure it looks right

In [42]:
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');
In [43]:
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]

Created fgm object
Check that this has the expected number of fgmax points:  758485
In [44]:
fname_nc = '%s_input.nc' % loc
fgmax_routines.make_nc_input(fname_nc, fgm, force=True)
print('Created ',fname_nc)

Created SW_Whidbey_input.nc
History:   Created with input data Thu Sep 26 18:04:06 2019 in /Users/rjl/git/WA_EMD_2019/Runs/SW_Whidbey;  
Created  SW_Whidbey_input.nc

Save resulting input files

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.