Plot fgmax results

Plot fgmax results#

This Jupyter notebook is available in $GTT/CopalisBeach/example3/plot_fgmax_grid_Copalis.ipynb within the GeoClaw Tsunami Tutorial.

The notebook plot_fgmax_grid_Copalis.ipynb creates an fgmax grid for use with point_style==4 in GeoClaw, covering the points around Copalis Beach that are below some specified elevation.

In this notebook we read these results in and make some plots. This notebook is based on ../example2/plot_fgmax_results and produces similar plots.

from pylab import *
import os
from clawpack.geoclaw import fgmax_tools
from clawpack.visclaw import geoplot, gridtools, plottools

if 0:
    # use this to fetch sample_results from the online data repository:
    if not os.path.isdir('sample_results'):
        import fetch_sample_results  # fetches directory
    outdir = 'sample_results/_output'
else:
    # use this if you have run the code locally to create '_output'
    outdir = '_output'

fgmax = fgmax_tools.FGmaxGrid()

fgmax.outdir = outdir  # as set at top of this notebook

# read the input data used for this run:
data_file = os.path.join(fgmax.outdir, 'fgmax_grids.data')
fgmax.read_fgmax_grids_data(fgno=1, data_file=data_file)

# read fgmax results:
fgmax.read_output()

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[3], line 7
      5 # read the input data used for this run:
      6 data_file = os.path.join(fgmax.outdir, 'fgmax_grids.data')
----> 7 fgmax.read_fgmax_grids_data(fgno=1, data_file=data_file)
      9 # read fgmax results:
     10 fgmax.read_output()

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/clawpack/geoclaw/fgmax_tools.py:93, in FGmaxGrid.read_fgmax_grids_data(self, fgno, data_file)
     87 def read_fgmax_grids_data(self, fgno, data_file='fgmax_grids.data'):
     88     """
     89     Read input info for fgmax grid number fgno from the data file
     90     fgmax_grids.data, which should have been created by setrun.py.
     91     This file now contains info about all fgmax grids.
     92     """
---> 93     with open(data_file) as filep:
     94         lines = filep.readlines()
     95     #print('+++ opened fgmax_new.data with %i lines' % len(lines))

FileNotFoundError: [Errno 2] No such file or directory: '_output/fgmax_grids.data'

Note that the cell above is identical to the one that reads the fgmax data from a uniform 2D grid in ../example2/plot_fgmax_results. The code in fgmax_tools converts the data that is output only at the points indicated as fgmax points in the file fgmax_pts_CopalisBeach.data into 2D arrays that can be worked with as usual.

Compute pre-earthquake topography#

dtopodir = '../../dtopo/dtopofiles'
dtopofile = os.path.join(dtopodir, 'ASCE_SIFT_Region2.dtt3')
fgmax.interp_dz(dtopofile, dtopo_type=3)
fgmax.B0 = fgmax.B - fgmax.dz

Plot onshore maximum water depth#

In this notebook we concentrate on illustrating some features of the fgmax solution computed using point_style==4 and so we just plot the solution without any underlying image. See ../example2/plot_fgmax_results for such plots.

# colormap:
clines = [0.01] + list(arange(3,31,3))
nlines = len(clines)
n1 = int(floor((nlines-1)/2.))
n2 = nlines - 1 - n1
Green = hstack([linspace(1,1,n1),linspace(1,0,n2)])
Red = hstack([linspace(0,0.8,n1), ones(n2)])
Blue = hstack([linspace(1,0.2,n1), zeros(n2)])
Alpha = 1.0*ones(nlines)  # transparency
colors = list(zip(Red,Green,Blue,Alpha))

onshore = fgmax.B0 > 0
h_onshore = where(onshore, fgmax.h, nan)
h_masked = where(fgmax.h.mask, 1, 0)
h_dry = where(fgmax.h < 0.01, 1, 0)


fig,ax = subplots(figsize=(6,10))
#contourf(fgmax.X,fgmax.Y, h_dry, [.5, 1.5], colors='g')
contourf(fgmax.X,fgmax.Y, h_masked, [.5, 1.5], colors=[.8,.8,.8])

cs = contourf(fgmax.X,fgmax.Y,h_onshore,clines,colors=colors,
              extend='both')
cs.cmap.set_over('m')
cs.cmap.set_under('g')
cs.changed()

fig.colorbar(cs, shrink=0.7, label='meters')

contourf(fgmax.X, fgmax.Y, fgmax.B0, [-1000,0], colors=['b'])
contour(fgmax.X, fgmax.Y, fgmax.B0, [0], colors='b', linewidths=0.7)
title('Max water depth onshore')

# fix axes:
ticklabel_format(style='plain',useOffset=False)
xticks(rotation=20)
gca().set_aspect(1./cos(fgmax.Y.mean()*pi/180.));

#axis([-124.195, -124.155, 47.11, 47.145]) # to compare with fgmax results from example2

In the plot above:

Blue regions show regions where B0 < 0 (the ocean and river) without any indication of water depth
Dark green shows regions that were captured in the fgmax results but where fgmax.h==0 (stayed dry)
Gray shows onshore regions that were not included in the fgmax points being monitored (presumably stayed dry due to the high elevation)
Other colors show the maximum water depth as recorded in the fgmax results.

Things to note:

Onshore regions that flooded are mostly surrounded by a dark green border, which indicates that the fgmax points selected went up to a sufficiently high elevation to capture the inundation.
The maximum water depth goes to zero (cyan color) at the margins of the inundated region, as expected.

If such plots show regions where the water depth does not seem to go to zero near the margins, or where there is no surrounding dark green cells, this may be an indication that this event reached higher elevations than expected and the run should be re-done after selecting more fgmax points in the Make fgmax_grid data for Copalis Beach notebook.

In the plot above, we might be concerned about the region around latitude 47.15, where there are some steep cliffs onshore, and also in the steep valley near the northern edge of the domain. These concerns are reinforced by the next plot…

Maximum water elevation (runup) onshore#

For inundation maps one typically wants to plot the maximum water depth as shown above. However, for some purposes it is also interesting to plot the maximum water surface elevation, as in the plot below.

fig,ax = subplots(figsize=(6,10))
#imshow(fg_image, extent=fg_extent)

h_masked = where(fgmax.h.mask, 1, 0)
h_dry = where(fgmax.h < 0.01, 1, 0)
#contourf(fgmax.X,fgmax.Y, h_dry, [.5, 1.5], colors='g')
contourf(fgmax.X,fgmax.Y, h_masked, [.5, 1.5], colors=[.8,.8,.8])

zeta_wet_onshore = where(logical_and(fgmax.h > 0.01, logical_not(fgmax.h.mask)), h_onshore+fgmax.B0, 0.)
cs = contourf(fgmax.X, fgmax.Y, zeta_wet_onshore,
              clines, colors=colors, extend='max')
cs.cmap.set_over('m')
cs.changed()

fig.colorbar(cs, shrink=0.7, label='meters')

contourf(fgmax.X, fgmax.Y, fgmax.B0, [-1000,0], colors=['b'])
contour(fgmax.X, fgmax.Y, fgmax.B0, [0], colors='b', linewidths=0.7)
title('Max water surface elevation onshore')

# fix axes:
ticklabel_format(style='plain',useOffset=False)
xticks(rotation=20)
gca().set_aspect(1./cos(fgmax.Y.mean()*pi/180.));