GeoClaw dtopotools examples¶
This Jupyter notebook can be found in collection of Clawpack apps as the file $CLAW/apps/notebooks/geoclaw/dtopotools_examples.ipynb.
To run this notebook, install Clawpack, and clone the apps repository.
A static view of this and other notebooks can be found in the Clawpack Gallery of Jupyter notebooks.
This notebook contains some examples of working with the clawpack.geoclaw.dtopotools module. These tools facilitate creating and manipulating the dtopo files that are required as GeoClaw input to specify moving topography, in particular the seafloor motion due to a subduction zone earthquake. All the examples in this notebook are of this nature: slip on a fault surface (or in geneneral a collection of rectangular planar subfaults making up the fault) is converted to seafloor deformation by applying the Okada model.
See http://clawpack.github.io/doc/dtopotools.html for general documentation of these tools and their use in the context of GeoClaw, and http://www.clawpack.org for more information on the Clawpack suite of software.
Setup notebook and environment:¶
%pylab inline
from clawpack.geoclaw import dtopotools, topotools
import os
CLAW = os.environ['CLAW']
datadir = os.path.join(CLAW,'geoclaw','scratch') # directory for some sample data files
Reading a csv file specifying the subfaults¶
Read in a csv file where the first line labels the columns and the remaining lines give the parameters for each subfault. This is a model of the 1964 Alaska earthquake using 12 subfaults from the NOAA SIFT database of unit sources.
See http://nctr.pmel.noaa.gov/propagation-database.html and gica2937.pdf.
The full database is available at http://sift.pmel.noaa.gov/ComMIT/compressed/info_sz.dat and a later example in this notebook illustrates using the dtopotools.SiftFault class to access this directly. But this example shows how a csv file can be specified. The first line labels the columns with standard attributes needed to specify each subfault.
Here is the raw file:
subfault_fname = 'data/alaska1964.csv'
print(open(subfault_fname).read())
Note that the last column has the name of the unit source from the NOAA SIFT database and will be ignored. The other columns have labels that are expected. Specify the units used for each and the coordinate system for the (longitude,latitude) and depth when reading the file in:
input_units = {"length":"km", "width":"km", "depth":"km", "slip":"m"}
fault = dtopotools.CSVFault()
fault.read(subfault_fname, input_units=input_units, coordinate_specification="noaa sift")
You could skip directly to Create the dtopo file needed by GeoClaw, but first we'll explore this fault model a bit...
The Seismic moment Mo and Moment Magnitude Mw can be computed using the Fault class methods:
print("The seismic moment is %g N-m" % fault.Mo())
print("The Moment magnitude is %g" % fault.Mw())
print(" (Assuming the rigidity mu of all subfaults is the default value %g Pa)"\
% fault.subfaults[0].mu)
Plot the subfault geometry and slip¶
First with a vector showing the rake (direction of slip):
fault.plot_subfaults(plot_rake=True)
Plot again showing the magnitude of slip on each subfault:
fault.plot_subfaults(slip_color=True)
Plot the coastline by downloading a file containing the x-y coordinates of shorelines based on etopo 4 minute data:
from clawpack.clawutil.data import get_remote_file
filename = 'pacific_shorelines_east_4min.npy'
url = 'http://www.geoclaw.org/topo/' + filename
get_remote_file(url=url, output_dir=datadir, force=True, verbose=True)
shorelines_file = os.path.join(datadir, filename)
fault.plot_subfaults()
shore = load(shorelines_file)
plot(shore[:,0], shore[:,1], 'g')
axis([170,240,40,65])
gca().set_aspect(1./cos(55*pi/180.))
Create the dtopo file needed by GeoClaw:¶
Generate a dtopo grid by applying the Okada model to each subfault. First we need to specify the region over which the seafloor deformation should be specified and the resolution of the grid:
xlower = 203.5
xupper = 214. # approximate - adjusted below
ylower = 54.5
yupper = 60. # approximate - adjusted below
# dtopo parameters:
points_per_degree = 60 # 1 minute resolution
dx = 1./points_per_degree
mx = int((xupper - xlower)/dx + 1)
xupper = xlower + (mx-1)*dx
my = int((yupper - ylower)/dx + 1)
yupper = ylower + (my-1)*dx
x = linspace(xlower,xupper,mx)
y = linspace(ylower,yupper,my)
Now apply Okada to create the static deformation at a single time $t = 1$ second:
dtopo = fault.create_dtopography(x,y,times=[1.], verbose=True)
We can save this deformation as a dtopo file with various dtopo_types recognized by GeoClaw. The most compact is dtopo_type==3, which specifies a header followed by all the dZ data:
dtopo_fname = os.path.join(datadir, 'alaska1964.tt3')
dtopo.write(dtopo_fname, dtopo_type=3)
Read the file in and print just the header. (The remaining lines contain all the data.)
lines = open(dtopo_fname).readlines()
for k in range(9):
print(lines[k][:-1])
Plot the seafloor displacement¶
dtopo.plot_dZ_colors(t=1)
# add shoreline and Anchorage for orientation:
plot(shore[:,0], shore[:,1], 'g')
plot([210.1],[61.2],'ro')
text(210.3,61.5,'Anchorage')
axis([xlower,xupper,ylower,62])
Read in an existing dtopo file:¶
An existing dtopo file can be read in for plotting purposes or to further manipulate it. To illustrate, we read in the file we just created...
dtopo2 = dtopotools.DTopography()
dtopo2.read(dtopo_fname, dtopo_type=3)
# Check that this data looks right:
assert len(dtopo2.x) == 631, "*** length of x is wrong"
assert len(dtopo2.y) == 331, "*** length of y is wrong"
dz_max = abs(dtopo2.dZ).max()
assert abs(dz_max - 15.368266081250006) < 1e-5, "*** dz_max is wrong: %g" % dz_max
print("Looks ok")
Create a fault using the NOAA SIFT database directly¶
The example above showed how to read a csv file with arbitrary columns. Since this particular fault is actually specified in terms of the NOAA SIFT unit sources, another option to create the same fault is to use the dtopotools.SiftFault class, which takes as an argument a dictionary sift_slip specifying the unit sources to be used and the slip on each.
To illustrate, here we specify only the northern two and southern two unit sources from the example above, but with a dictionary of 12 unit sources we could recreate the full fault.
sift_slips = {'acsza31':19.5, 'acszb31':19.5, 'acsza36':34.5, 'acszb36':34.5}
f2 = dtopotools.SiftFault(sift_slips)
f2.plot_subfaults(slip_color=True)
Time-dependent fault rupture¶
This is synthetic data -- the same set of faults as before but with rupture occuring from south to north. Note that new columns rupture_time and rise_time have been added.
subfault_fname = 'data/timedep.csv'
print(open(subfault_fname).read())
Read in this csv file and set the parameters for the desired dtopo file:
input_units = {"length":"km", "width":"km", "depth":"km",
"slip":"m", "mu":"dyne/cm^2"}
fault = dtopotools.CSVFault()
fault.read(subfault_fname, input_units=input_units,
coordinate_specification="noaa sift")
fault.rupture_type = 'dynamic'
print("%s subfaults read in " % len(fault.subfaults))
xlower = 203.5
xupper = 214.
ylower = 54.5
yupper = 60.
xylim = [xlower,xupper,ylower,yupper]
# dtopo parameters:
points_per_degree = 60 # 1 minute resolution
mx = int((xupper - xlower)*points_per_degree + 1)
my = int((yupper - ylower)*points_per_degree + 1)
x = linspace(xlower,xupper,mx)
y = linspace(ylower,yupper,my)
print("dZ arrays will have shape %s by %s" % (len(y),len(x)))
Specify the desired output times, in this case 6 times between $t=0$ and $t=13$, when the rupture motion is complete. Then apply the Okada model at each time (to the set of all subfaults with slips computed at each time based on the rupture_time and rise_time for that subfault). Then dtopo.dZ will be a 3 dimensional array where each field at 6 different time points will be represented:
times = linspace(0,13,6)
dtopo = fault.create_dtopography(x,y,times=times)
print("Created array of dtopo deformations with shape", dtopo.dZ.shape)
Create an animation of the fault rupture¶
First define a function that will show time-dependent slip on the fault planes and resulting time-dependent seafloor motion in two side-by-side axes in the same figure:
# for setting color scale:
print("maximum abs(dz) over the full rupture time:", \
abs(dtopo.dZ).max())
# Read in shoreline segments to give context in plots:
#shoreline_fname = os.path.join(datadir,'tohoku_shoreline_1min.npy')
#shoreline_xy = load(shoreline_fname)
shoreline_xy = shore
# Incorporate this function in dtopotools to replace animate_dz_colors?
def plot_subfaults_dZ(t, fig):
fig.clf()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
fault.plot_subfaults(axes=ax1, slip_color=True,
slip_time=t, xylim=xylim)
dtopo.plot_dZ_colors(axes=ax2, t=t, cmax_dZ=dz_max)
ax1.plot(shoreline_xy[:,0],shoreline_xy[:,1],'g')
ax2.plot(shoreline_xy[:,0],shoreline_xy[:,1],'g')
ax2.set_xlim(xlower,xupper)
ax2.set_ylim(ylower,yupper)
return fig
The loop below makes each plot and then saves it as a png file:
from clawpack.visclaw.JSAnimation import IPython_display
import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J
plotdir = '_plots'
J.make_plotdir(plotdir, clobber=True)
fig = plt.figure(figsize=(12,5))
for k,t in enumerate(times):
plot_subfaults_dZ(t,fig)
J.save_frame(k, verbose=True)
print("Final frame will be displayed below:")
Combine the png files into an animation and display it:
anim = J.make_anim(plotdir)
anim