Download and plot hmax results

Original Author: R.J. LeVeque, Modified by Christopher Liu to work also on Windows.

For data from the paper

"A Source Clustering Approach for Efficient Inundation Modeling and Regional Scale PTHA" by A.L. Williamson, D. Rim, L.M. Adams, R.J. LeVeque, D. Melgar, and F.I. González, Frontiers in Earth Science, 8 (2020) p. 442. [paper and more data].

This illustrates how to download the hmax (maximum inundation) results for Westport from some of the 2000 GeoClaw fine-grid runs archived at http://krakatoa.uoregon.edu/cascadia_ptha/ground_truth/Westport_fine_hfiles. These give the maximum inundation at grid points on a 1 arcsecond grid covering part of Westport, WA.

The corresponding coarse-grid run results (on a 9 arcsecond grid) are in the directory http://krakatoa.uoregon.edu/cascadia_ptha/coarse_model_runs/Westport_coarse_hfiles.

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [2]:
import os
from urllib.request import urlretrieve

Load the topography

This loads topo data on the same spatial grids as the hmax files:

In [3]:
X = load('B0/fgmax0001_1sec_X.npy')
Y = load('B0/fgmax0001_1sec_Y.npy')
B0 = load('B0/fgmax0001_1sec_B0.npy')
print('These arrays have shape ', B0.shape)

These arrays have shape  (216, 180)

Function to download hmax file for a set of runs

In [4]:
# where to put the files downloaded  
datadir = os.path.join(os.getcwd(), 'westport_data')

def get_hmax_urllib(rnums):

    # create directory for data if it does not exist
    if not os.path.isdir(datadir):
        os.mkdir(datadir)
    
    for rnum in rnums:
        fname = 'h_%s.npy' % str(rnum).zfill(6)
        url = 'http://krakatoa.uoregon.edu/cascadia_ptha/ground_truth/Westport_fine_hfiles/'
        rundir = os.path.join(datadir, 'run%s' % str(rnum).zfill(6))
        
        if not os.path.isdir(rundir):
            os.mkdir(rundir)
        
        print('Data for run %s will go into directory %s' % (rnum, rundir))
        
        local_fname = os.path.join(rundir, 'hmax.npy')
        if os.path.isfile(local_fname):
            print('%s has already been downloaded' % fname)
        else:
            remote_fname = url + '/' + fname
            print('Downloading %s from %s' % (fname,url))
            print('Saving as %s' % local_fname)
            a = urlretrieve(remote_fname, local_fname)

For illustration purposes we will just download a few samples, including the ones shown in the paper.

The script download_hmax.py can be used to download all 2000 or whatever subset you want.

In [5]:
get_hmax_urllib([1665, 1670, 1999])

Data for run 1665 will go into directory /Users/rjl/git/Forks/jdf-autoencoder/frontiers2020_westport/westport_data/run001665
h_001665.npy has already been downloaded
Data for run 1670 will go into directory /Users/rjl/git/Forks/jdf-autoencoder/frontiers2020_westport/westport_data/run001670
h_001670.npy has already been downloaded
Data for run 1999 will go into directory /Users/rjl/git/Forks/jdf-autoencoder/frontiers2020_westport/westport_data/run001999
Downloading h_001999.npy from http://krakatoa.uoregon.edu/cascadia_ptha/ground_truth/Westport_fine_hfiles/
Saving as /Users/rjl/git/Forks/jdf-autoencoder/frontiers2020_westport/westport_data/run001999/hmax.npy

Load an hmax file:

In [6]:
rnum = 1670
rundir = os.path.join(datadir, 'run%s' % str(rnum).zfill(6))
fname = os.path.join(rundir, 'hmax.npy')
hmax = load(fname)
print('hmax array has shape ', hmax.shape)

hmax array has shape  (216, 180)

Define array with nan values offshore

hmax includes the depth at offshore points. If you only want to plot onshore flooding depth, then you can replace offshore points with NaN:

In [7]:
hmax_onshore = where(B0>0, hmax, nan)

Sample plot with default colormap

In [8]:
pcolor(X,Y,hmax_onshore)
colorbar()
gca().set_aspect(1/cos(47*pi/180)) # fix aspect ratio based on latitude

Nicer plots

This uses some tools from visclaw to facilitate setting up color maps and plotting cell-based values, but once you have read in the data, you can manipulate and plot with whatever tools you want to use.

In [9]:
from clawpack.visclaw import colormaps, plottools
In [10]:
# Color map for topography / bathymetry:

zmin = -15.
zmax = 15.

cmap_land = 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]})

cmap_water = colormaps.make_colormap({ 0.0:[0,0,1], 1.:[.8,.8,1]})

cmap_topo, norm_topo = colormaps.add_colormaps((cmap_land, cmap_water),
                                     data_limits=(zmin,zmax),
                                     data_break=0.)

# Color map for hmax:

clines = [0,.5,1.0] + list(linspace(2.0,16.0,8))
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)])
hmax_colors = list(zip(Red,Green,Blue))
In [11]:
figure(figsize=(10,10))
plottools.pcolorcells(X,Y,B0,cmap=cmap_topo,norm=norm_topo)
colorbar(extend='both')
gca().set_aspect(1/cos(47*pi/180))
title('Westport, WA -- pre-seismic topography');
In [12]:
def plot_hfile(rnum):
    rundir = os.path.join(datadir, 'run%s' % str(rnum).zfill(6))
    fname = os.path.join(rundir, 'hmax.npy')
    hmax = load(fname)
    
    figure(figsize=(10,10))
    contour(X,Y,B0,[0,5,10], colors='k')
    hmax_onshore = where(B0>0, hmax, nan)
    cf = contourf(X,Y,hmax_onshore,clines,colors=hmax_colors)
    cbar = colorbar(label='meters')
    cbar.set_ticks(clines)
    title('Inundation depth for run %s' % str(rnum).zfill(6));
In [13]:
plot_hfile(1665)
In [14]:
plot_hfile(1670)
In [15]:
plot_hfile(1999)