<jupyter_start><jupyter_text>Make Input Files for Whidbey1 Test Problem<jupyter_code>%matplotlib inline
from pylab import *
from IPython.display import Image
import os,sys
from clawpack.visclaw import colormaps
sys.path.insert(0,'../../new_python')
import region_tools, topotools, marching_front, dtopotools
import fgmax_tools
from plottools import pcolorcells
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.)<jupyter_output><empty_output><jupyter_text>Download topoTopo data for this test problem is taken from the 1/3 arcsecond Puget Sound DEM referenced to MHW, available on the NCEI thredds server.First we download data for the full domain at 2" resolution by coarsening the 1/3" data:<jupyter_code>try:
    # recommended approach: read topo directly from NCEI thredds server:
    extent = [-123, -122.2, 47.2, 48.2]
    path = 'https://www.ngdc.noaa.gov/thredds/dodsC/regional/puget_sound_13_mhw_2014.nc'
    topo_2sec = topotools.read_netcdf(path, extent=extent, coarsen=6, verbose=True)
except:
    # in case that fails:
    from clawpack.clawutil.data import get_remote_file
    print('Reading from thredds server failed, instead read cached version from geoclaw...')
    remote_topo_dir = 'http://depts.washington.edu/clawpack/geoclaw/topo/new_features_for_v5.7.0'
    fname = 'topo_2sec.tt3'
    path = os.path.join(remote_topo_dir, fname)
    get_remote_file(path, output_dir='input_files', file_name=fname, 
                    verbose=True)
    topo_2sec = topotools.Topography()
    topo_2sec.read(os.path.join('input_files',fname))
figure(figsize=(12,8))
ax = axes()
topo_2sec.plot(axes=ax)<jupyter_output><empty_output><jupyter_text>Also download topo data over the small region of interest at full 1/3" resolution:<jupyter_code>try:
    # recommended approach: read topo directly from NCEI thredds server:
    extent = [-122.46, -122.38, 47.93, 47.96]
    path = 'https://www.ngdc.noaa.gov/thredds/dodsC/regional/puget_sound_13_mhw_2014.nc'
    topo_1_3sec = topotools.read_netcdf(path, extent=extent, verbose=True)
except:
    # in case that fails:
    from clawpack.clawutil.data import get_remote_file
    print('Reading from thredds server failed, instead read cached version from geoclaw...')
    remote_topo_dir = 'http://depts.washington.edu/clawpack/geoclaw/topo/new_features_for_v5.7.0'
    fname = 'topo_1_3sec.tt3'
    path = os.path.join(remote_topo_dir, fname)
    get_remote_file(path, output_dir='input_files', file_name=fname, 
                    verbose=True)
    topo_1_3sec = topotools.Topography()
    topo_1_3sec.read(os.path.join('input_files',fname))
figure(figsize=(12,8))
ax = axes()
topo_1_3sec.plot(axes=ax)<jupyter_output><empty_output><jupyter_text>Save the topofiles for use in GeoClaw:<jupyter_code>fname = 'input_files/topo_2sec.tt3'
topo_2sec.write(fname, topo_type=3)

fname = 'input_files/topo_1_3sec.tt3'
topo_1_3sec.write(fname, topo_type=3)<jupyter_output><empty_output><jupyter_text>Determine fgmax pointsWe will choose fgmax points by selecting points that are<jupyter_code>topo = topo_1_3sec
pts_chosen = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=1e6, max_iters=20) 
pts_chosen = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=15., 
                                               prev_pts_chosen=pts_chosen,
                                               max_iters=None) 
pts_chosen_shallow = marching_front.select_by_flooding(topo.Z, Z1=0, Z2=-15., max_iters=None) 
pts_chosen_nearshore = logical_and(pts_chosen, pts_chosen_shallow)
Znearshore = ma.masked_array(topo.Z, logical_not(pts_chosen_nearshore))
print('Total number of fgmax points = %i' % Znearshore.count())
figure(figsize=(12,6))
pcolorcells(topo.X, topo.Y, Znearshore, cmap=cmap, norm=norm)
colorbar(extend='both')
gca().set_aspect(1./cos(48*pi/180.))<jupyter_output><empty_output><jupyter_text>Write out fgmax points as an array to be read into GeoClaw:<jupyter_code>fname_fgmax_mask = 'input_files/fgmax_pts_whidbey1.data'
topo_fgmax_mask = topotools.Topography()
topo_fgmax_mask._x = topo.x
topo_fgmax_mask._y = topo.y     
topo_fgmax_mask._Z = where(pts_chosen_nearshore, 1, 0)  # change boolean to 1/0
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)<jupyter_output><empty_output><jupyter_text>Create an AMR flagregion covering the fgmax points<jupyter_code>rr_fgmax = region_tools.ruledrectangle_covering_selected_points(topo.X, topo.Y, 
                                                                pts_chosen_nearshore, 
                                                                ixy='y', method=0,
                                                                padding=0, verbose=True)
figure(figsize=(12,6))
xv_whidbey1,yv_whidbey1 = rr_fgmax.vertices()
pcolorcells(topo.X, topo.Y, Znearshore, cmap=cmap, norm=norm)
axis([-122.47, -122.40, 47.925, 47.965])
gca().set_aspect(1./cos(48*pi/180.))
plot(xv_whidbey1, yv_whidbey1, 'r')
title("With ixy = 'y'")<jupyter_output><empty_output><jupyter_text>Write this ruled rectangle out as a file that can be read into GeoClaw:<jupyter_code>rr_name = 'input_files/RuledRectangle_whidbey1'
rr_fgmax.write(rr_name + '.data', verbose=True)
rr_fgmax.make_kml(rr_name + '.kml', verbose=True)<jupyter_output><empty_output><jupyter_text>Create an AMR flagregion covering Admiralty InletWe also create a ruled rectangle covering Admiralty Inlet, using code from [RuledRectangles.ipynb](RuledRectangles.ipynb):<jupyter_code>slu = \
array([[  47.851, -122.75 , -122.300],
       [  47.955, -122.75 , -122.300],
       [  48.   , -122.8  , -122.529],
       [  48.036, -122.8  , -122.578],
       [  48.12 , -122.9  , -122.577],
       [  48.187, -122.9  , -122.623],
       [  48.191, -122.9  , -122.684],
       [  48.221, -122.9  , -122.755]])
       
rr_admiralty = region_tools.RuledRectangle(slu=slu)
rr_admiralty.ixy = 'y'
rr_admiralty.method = 1

rr_name = 'input_files/RuledRectangle_AdmiraltyInlet'
rr_admiralty.write(rr_name + '.data', verbose=True)
rr_admiralty.make_kml(rr_name + '.kml', verbose=True)

xv_admiralty, yv_admiralty = rr_admiralty.vertices()<jupyter_output><empty_output><jupyter_text>Plot these two flagregions:<jupyter_code>figure(figsize=(12,8))
ax = axes()
topo_2sec.plot(axes=ax)
plot(xv_whidbey1, yv_whidbey1, 'r')
plot(xv_admiralty, yv_admiralty, 'yellow')
ylim(47.8, 48.2)<jupyter_output><empty_output><jupyter_text>Make `force_dry_init` arrayThe code below is taken from the notebook [ForceDry.ipynb](ForceDry.ipynb) found in the notebooks directory.  See that notebook for more description of what is being done in each step, and some additional plots.<jupyter_code>wet_points = marching_front.select_by_flooding(topo.Z, Z1=-5., Z2=0., max_iters=None)

dry_points = 1 - wet_points
dry_points_sum = dry_points[1:-1,1:-1] + dry_points[0:-2,1:-1] + dry_points[2:,1:-1] + \
                 dry_points[1:-1,0:-2] + dry_points[0:-2,0:-2] + dry_points[2:,0:-2] + \
                 dry_points[1:-1,2:] + dry_points[0:-2,2:] + dry_points[2:,2:]
        
# initialize array to 0 everywhere:
force_dry_init = zeros(dry_points.shape)
# reset in interior to 1 if all points in the 3x3 block around it are dry:
force_dry_init[1:-1,1:-1] = where(dry_points_sum == 9, 1, 0)

force_dry_init_topo = topotools.Topography()
force_dry_init_topo.set_xyZ(topo.x, topo.y, force_dry_init)

fname_force_dry_init = 'input_files/force_dry_init.data'
force_dry_init_topo.write(fname_force_dry_init, topo_type=3, Z_format='%1i')
print('Created %s' % fname_force_dry_init)
figure(figsize=(10,12))

# first plot all dry points as pink:

# Create a version of topo.Z with wet points masked out:
mask_dry = logical_not(wet_points)
Z_dry = ma.masked_array(topo.Z, wet_points) 
pcolorcells(topo.X, topo.Y, Z_dry, cmap=cmap_dry, norm=norm_dry)

# then plot colored by topography, except at dry points below MHW:

# Create a version of topo.Z with only dry points below MHW masked out:
mask_dry_onshore = logical_and(mask_dry, topo.Z<0.)
Z_allow_wet= ma.masked_array(topo.Z, mask_dry_onshore)
pcolorcells(topo.X, topo.Y, Z_allow_wet, cmap=cmap, norm=norm)

gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);<jupyter_output><empty_output><jupyter_text>Create dtopo file for earthquake sourceThis is the Seattle Fault "large" earthquake denoted by SF-L in our recent modeling work. The fault cuts across Puget Sound between Bainbridge Island and Seattle.<jupyter_code>subfault_fname = 'input_files/SFL.csv'
dtopo_fname = 'input_files/SFL.tt3'
print("Using Okada model to create %s " % dtopo_fname)

defaults = {}
coordinate_specification = 'top center'
input_units = {'slip': 'm', 'depth': 'km', 
           'length': 'km', 'width': 'km'}
rupture_type = 'static'

fault = dtopotools.CSVFault()
fault.read(subfault_fname, input_units, coordinate_specification,
       rupture_type)
print("There are %s subfaults" % len(fault.subfaults))
print("Mw = %5.2f" % fault.Mw())

x = arange(-123, -122.2, 2./3600.)
y = arange(47., 48.6, 2./3600.)
times = [1.]

fault.create_dtopography(x,y,times)
dtopo = fault.dtopo
dtopo.write(dtopo_fname, dtopo_type=3)

from importlib import reload
reload(dtopotools)
figure(figsize=(12,8))
ax = axes()
contourf(topo_2sec.X, topo_2sec.Y, topo_2sec.Z, [0,1000],colors=[[.4,1,.4]])
dtopo.plot_dZ_colors(2, axes=ax, dZ_interval=0.5)
ax.set_ylim(47.2,48.2)
ax.set_aspect(1./cos(48*pi/180.))<jupyter_output><empty_output><jupyter_text>Make netCDF file with inputCopy arrays into a new `fgmax_tools.FGmaxMaskedGrid` object and then write this out.First crop arrays down to bounding box of fgmax points:<jupyter_code>Z_fg = Znearshore 
print('original shape of X,Y,Z: ', Z_fg.shape)
print('original number of fgmax points: %i' % Z_fg.count())
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]
Xcrop = topo.X[i1:i2+1, j1:j2+1]
Ycrop = topo.Y[i1:i2+1, j1:j2+1]
force_dry_init_crop = force_dry_init[i1:i2+1, j1:j2+1]
print('new shape of X,Y,Z: ', Zcrop.shape)
print('new number of fgmax points: %i' % Zcrop.count())<jupyter_output><empty_output><jupyter_text>Make a final check that fgmax points look correct:<jupyter_code>figure(figsize=(10,12))
pcolorcells(Xcrop, Ycrop, Zcrop, cmap=cmap, norm=norm)
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('%i fgmax points on %i by %i grid' \
      % (Znearshore.count(), Zcrop.shape[0], Zcrop.shape[1]));
fgm = fgmax_tools.FGmaxMaskedGrid()
fgm.id = 'Whidbey1'
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 force_dry_init is not None:
    fgm.force_dry_init = force_dry_init_crop
print('new bounding box: ', fgm.bounding_box())
fname_nc = 'input.nc'
fgmax_tools.make_nc_input(fname_nc, fgm, force=True)
print('Created ',fname_nc)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13220338983050847, 'fraction_numerical': 0.043819760231500623, 'mean_word_length': 4.09519797809604, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 57, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '576565', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5116, 'n_tokens_neox': 4780, 'n_words': 890, 'path': 'examples/geoclaw_whidbey1/MakeInputFiles_Whidbey1.ipynb', 'repo_name': 'clawpack/new_features_for_v5.7.0', 'size': '12095'}