<jupyter_start><jupyter_text>process fgmax outputAfter running the GeoClaw code, use this notebook to read in and plot the fgmax results.Change the next cell to `%matplotlib notebook` if you want the plots to be live in the notebook (so you can zoom in on them, for example), but then remember to close figures when you are done with them.<jupyter_code>%matplotlib inline
from pylab import *
from clawpack.geoclaw import fgmax_tools, topotools, dtopotools
from clawpack.visclaw import colormaps
from scipy.interpolate import RegularGridInterpolator
import numpy
GEmap = imread('EagleHarborGE.jpg')
GEextent = [-122.55,-122.48,47.61,47.64]
figure(figsize=(8,6))
imshow(GEmap,extent=GEextent)
gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
fg = fgmax_tools.FGmaxGrid()
fg.read_input_data('fgmax_eagle_harbor.txt')
fg.read_output(outdir='_output')
fg.h.shape<jupyter_output><empty_output><jupyter_text>Define colormaps:see also https://matplotlib.org/tutorials/colors/colormaps.html<jupyter_code>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.)<jupyter_output><empty_output><jupyter_text>Compute subsidence during earthquakeThe values stored in `fg.B` are the cell-averaged topography values used in GeoClaw *after* uplift or subsidence due to the earthquake.  Since Eagle Harbor lies right on top of the Seattle Fault, there is significant subsidence in this region.  It is important to correct for this before plotting the "maximum flow depth on shore", for example, since it affects what cells are considered "on shore", as we illustrate below.  We generally want to show the flooding depth in regions that are on shore before the earthquake, even if they end up below sea level afterward.Read in the sea floor deformation file and interpolate to the locations of the fgmax grid.  **Note this is not exactly the right correction** in cases like this where the uplift varies greatly over a small area, but is pretty good and fine for most cases.<jupyter_code>dtopo_path = '../dtopo/dtopofiles/seattlefault_uniform.tt3'
dtopo = dtopotools.DTopography(dtopo_path, dtopo_type=3)
x1d = dtopo.X[0,:]
y1d = dtopo.Y[:,0]
dtopo_func = RegularGridInterpolator((x1d,y1d), dtopo.dZ[-1,:,:].T, 
                method='linear', bounds_error=False, fill_value=0.)
        
dz = dtopo_func(list(zip(numpy.ravel(fg.X), numpy.ravel(fg.Y))))
dz = numpy.reshape(dz, fg.X.shape)

# Estimate of topography before earthquake:
B0 = fg.B - dz

print('The maximum subsidence over the fgmax region is %.2f meters' % dz.min())
print('The maximum uplift over the fgmax region is %.2f meters' % dz.max())<jupyter_output><empty_output><jupyter_text>There is massive uplift just south of Eagle Harbor but subsidence on the harbor shore, as seen in this plot of the deformation (with the original coast in green):<jupyter_code>figure(figsize=(8,6))
ax = axes()
contour(fg.X, fg.Y, B0, [0], colors='g')  # original shoreline
dtopo.plot_dZ_colors(t=2, axes=ax)
axis(GEextent)
gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);<jupyter_output><empty_output><jupyter_text>Compare topography before and after event<jupyter_code>figure(figsize=(8,10))
subplot(2,1,1)
pc = pcolormesh(fg.X, fg.Y, B0, cmap=cmap, norm=norm)  
cb = colorbar(pc, extend='both', shrink=0.5)
cb.set_label('meters')
gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
title('GeoClaw topography B, before earthquake')

subplot(2,1,2)
pc = pcolormesh(fg.X, fg.Y, fg.B, cmap=cmap, norm=norm)  
cb = colorbar(pc, extend='both', shrink=0.5)
cb.set_label('meters')
gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
title('GeoClaw topography B, after earthquake')<jupyter_output><empty_output><jupyter_text>Compare the shorelines<jupyter_code>figure(figsize=(8,6))
imshow(GEmap,extent=GEextent)
contour(fg.X, fg.Y, B0, [0], colors='b')
contour(fg.X, fg.Y, fg.B, [0], colors='r')

gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
title('Shore before (blue) and after (red) earthquake')<jupyter_output><empty_output><jupyter_text>Define onshore depth based on B0Note that you get misleading results if you replace `B0` by `fg.B` below!<jupyter_code>h_onshore = ma.masked_where(B0 < 0., fg.h)
print('Maximum flow depth onshore: %.2f m' % h_onshore.max())
#bounds_depth = array([1e-6,0.25,0.5,0.75,1,1.25,1.5])
bounds_depth = array([1e-6,0.5,1.,2.,3.,4.,5.])

cmap_depth = mpl.colors.ListedColormap([[.7,.7,1],\
                 [.5,.5,1],[0,0,1],\
                 [1,.7,.7], [1,.4,.4], [1,0,0]])

# Set color for value exceeding top of range to purple:
cmap_depth.set_over(color=[1,0,1])

# Set color for land points without inundation to light green:
cmap_depth.set_under(color=[.7,1,.7])

norm_depth = mpl.colors.BoundaryNorm(bounds_depth, cmap_depth.N)
    
figure(figsize=(8,6))
pc = pcolormesh(fg.X, fg.Y, h_onshore, cmap=cmap_depth, norm=norm_depth)
cb = colorbar(pc, extend='max', shrink=0.5)
cb.set_label('meters')
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('Maximum flow depth')
#savefigp('h_onshore.png')<jupyter_output><empty_output><jupyter_text>Plot on map imageTo plot on a map we want to mask out the onshore dry points so they don't all get colored green and obscure the map.(Again, observe what happens is you replace `B0` by `fg.B`).<jupyter_code>h_onshore = ma.masked_where(B0 < 0, fg.h)

# define h_wet_onshore to suppress plotting dry onshore fgmax points:
h_wet_onshore = ma.masked_where(fg.h < 0.01, h_onshore)

figure(figsize=(8,6))
imshow(GEmap,extent=GEextent)

pc = pcolormesh(fg.X, fg.Y, h_wet_onshore, cmap=cmap_depth, norm=norm_depth)
cb = colorbar(pc, extend='max', shrink=0.5)
cb.set_label('meters')
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('Maximum flow depth')<jupyter_output><empty_output><jupyter_text>Plot the maximum speed<jupyter_code>print('Maximum flow speed: %.2f m/s' % fg.s.max())
bounds_speed = np.array([1e-6,0.5,1.5,2,2.5,3,4.5,6])
cmap_speed = mpl.colors.ListedColormap([[.9,.9,1],[.6,.6,1],\
                 [.3,.3,1],[0,0,1],\
                 [1,.7,.7], [1,.4,.4], [1,0,0]])

# Set color for value exceeding top of range to purple:
cmap_speed.set_over(color=[1,0,1])

# Set color for land points without inundation to light green:
cmap_speed.set_under(color=[.7,1,.7])

#imshow(GEmap,extent=GEextent)
norm_speed = mpl.colors.BoundaryNorm(bounds_speed, cmap_speed.N)

figure(figsize=(8,6))
pc = pcolormesh(fg.X, fg.Y, fg.s, cmap=cmap_speed, norm=norm_speed)
cb = colorbar(pc, extend='max')
cb.set_label('m/s')

contour(fg.X, fg.Y, B0, [0], colors='k')  # plot original shoreline

gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('Maximum speed')<jupyter_output><empty_output><jupyter_text>Plot on the mapMask out the speed where the depth is 0 before plotting over the map.We also include a transparency `alpha` in the call to `pcolormesh` so you can see the image below.<jupyter_code>s_wet = ma.masked_where(fg.h < 0.01, fg.s)

figure(figsize=(8,6))
imshow(GEmap,extent=GEextent)

pc = pcolormesh(fg.X, fg.Y, s_wet, alpha = 0.5, edgecolors='',
                cmap=cmap_speed, norm=norm_speed)
cb = colorbar(pc, extend='max')
cb.set_label('m/s')
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('Maximum speed')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14627127127127126, 'fraction_numerical': 0.03778778778778779, 'mean_word_length': 4.973841554559043, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 41, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '408612', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3199, 'n_tokens_neox': 3063, 'n_words': 689, 'path': 'eagle_harbor_test/process_fgmax.ipynb', 'repo_name': 'rjleveque/WA_DNR_GeoClaw', 'size': '7992'}