Create dtopo file for AKmaxWA source

Create the GeoClaw dtopo file AKmaxWA.tt3 for use in tsunami modeling projects for Washington State.

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [2]:
import os,sys
from clawpack.geoclaw import topotools,dtopotools
In [3]:
# topography for plots:
etopo1 = topotools.read_netcdf('etopo1', extent=[-160,-140,53,63], coarsen=3)

Using Clawpack 5.7.1 with with an updated version of https://sift.pmel.noaa.gov/ComMIT/compressed/info_sz.dat that will be added to v5.8.0, replacing the file clawpack/geoclaw/src/python/geoclaw/data/info_sz.dat.txt.

Check that some subfaults have been corrected:

In [4]:
AKmaxWA = dtopotools.SiftFault(longitude_shift=-360.)
SIFT = AKmaxWA.sift_subfaults

k = 'acsza38'
s = SIFT[k]
print('%s, dip = %i, depth = %.3f km' % (k,s.dip, s.depth/1e3))
assert s.dip==15, '*** wrong dip'
assert s.depth/1e3==17.94, '*** wrong depth'

k = 'acszz38'
s = SIFT[k]
print('%s, dip = %i, depth = %.3f km' % (k,s.dip, s.depth/1e3))
assert s.dip==15, '*** wrong dip'
assert s.depth/1e3==30.88, '*** wrong depth'

acsza38, dip = 15, depth = 17.940 km
acszz38, dip = 15, depth = 30.880 km

Note that we have also shifted longitudes by 360 to put into coordinates W used in WA state projects.

In [5]:
sift_slipAKmaxWA = {}
slip_uniform = 20.
for usrc in range(29,39):
    k = 'acsza%s' % usrc
    sift_slipAKmaxWA[k] = slip_uniform        
    k = 'acszz%s' % usrc
    sift_slipAKmaxWA[k] = slip_uniform
    
AKmaxWA.set_subfaults(sift_slipAKmaxWA)

print('With uniform slip %.2f meters, Magnitude = %.2f' % (slip_uniform,AKmaxWA.Mw()))

With uniform slip 20.00 meters, Magnitude = 9.24
In [6]:
print(' src       long      lat      dip   strike   rake  depth   slip')
for usrc in range(29,39):
    
    k = 'acsza%s' % usrc
    s = AKmaxWA.sift_subfaults[k]
    print('%s  %.4f  %.4f  %4.1f  %.4f  %.1f  %5.2f  %5.2f' \
        % (k, s.longitude, s.latitude, s.dip, s.strike, s.rake, (s.depth/1e3), s.slip))
    
    k = 'acszz%s' % usrc
    s = AKmaxWA.sift_subfaults[k]
    print('%s  %.4f  %.4f  %4.1f  %.4f  %.1f  %5.2f  %5.2f' \
        % (k, s.longitude, s.latitude, s.dip, s.strike, s.rake, (s.depth/1e3), s.slip))

 src       long      lat      dip   strike   rake  depth   slip
acsza29  -157.7390  55.1330  15.0  247.0000  90.0  17.94  20.00
acszz29  -158.1203  55.4908  15.0  246.2137  90.0  30.88  20.00
acsza30  -156.3960  55.5090  15.0  240.0000  90.0  17.94  20.00
acszz30  -156.8479  55.8534  15.0  240.4869  90.0  30.88  20.00
acsza31  -155.1050  55.9700  15.0  236.0000  90.0  17.94  20.00
acszz31  -155.5685  56.3016  15.0  235.6690  90.0  30.88  20.00
acsza32  -153.7920  56.4730  15.0  236.0000  90.0  17.94  20.00
acszz32  -154.2120  56.8210  15.0  235.4756  90.0  30.88  20.00
acsza33  -152.4630  56.9750  15.0  236.0000  90.0  17.94  20.00
acszz33  -152.8909  57.3227  15.0  235.4119  90.0  30.88  20.00
acsza34  -151.0629  57.5124  15.0  236.0000  90.0  17.94  20.00
acszz34  -151.5802  57.8213  15.0  234.6891  90.0  30.88  20.00
acsza35  -149.7403  58.0441  15.0  230.0000  90.0  17.94  20.00
acszz35  -150.3575  58.3252  15.0  230.1971  90.0  30.88  20.00
acsza36  -148.6751  58.6565  15.0  218.0000  90.0  17.94  20.00
acszz36  -149.4588  58.8129  15.0  217.3327  90.0  30.88  20.00
acsza37  -147.7495  59.2720  15.0  213.7100  90.0  17.94  20.00
acszz37  -148.3921  59.5820  15.0  214.2669  90.0  30.88  20.00
acsza38  -145.3445  60.1351  15.0  260.0800  90.0  17.94  20.00
acszz38  -145.4638  60.5429  15.0  259.0313  90.0  30.88  20.00
In [7]:
x1 = -160.783294677734
x2 = -139.783294677734
mx = 316
dx = (x2-x1) / (mx-1)

y1 = 50.64400100708008
y2 = 61.653900146484375
my = 301
dy = (y2-y1) / (my-1)

print('Will create %i by %i dtopo with spacing' % (mx,my))
print('  dx = %.14f,  dy = %.14f' % (dx,dy))

x_unif = linspace(x1, x2, mx)
y_unif = linspace(y1, y2, my)
dtopo_AKmaxWA = AKmaxWA.create_dtopography(x_unif, y_unif, times=[1.],verbose=True)

fname = 'AKmaxWA.tt3'
dtopo_AKmaxWA.write(fname, 3)
print('Created ', fname)

Will create 316 by 301 dtopo with spacing
  dx = 0.06666666666667,  dy = 0.03669966379801
Making Okada dz for each of 20 subfaults
0..1..2..3..4..5..6..7..8..9..10..11..12..13..14..15..16..17..18..19..
Done
Created  AKmaxWA.tt3
In [8]:
figure(figsize=(13,13))
contourf(etopo1.X, etopo1.Y, etopo1.Z, [0,10000], colors=[[0.6,1,0.6]])
contour(etopo1.X, etopo1.Y, etopo1.Z, [0], colors='g', linewidths=0.5)
plot([-149.9],[61.2],'ko')
text(-150.2,61.4,'Anchorage',color='k',fontsize=14)
axis([-158, -143, 54, 62])
gca().set_aspect(1/cos(57*pi/180))
grid(True)

clines_neg = arange(-10,0.,1.)
clines_pos = arange(1,16,1.)

dd = dtopo_AKmaxWA
dZ = dd.dZ[0,:,:]
contour(dd.X, dd.Y, dZ, clines_neg, colors='b', linestyles='-')
contour(dd.X, dd.Y, dZ, clines_pos, colors='r', linestyles='-')

dZmin = dZ.min()
dZmax = dZ.max()
axis([-160,-143,54,62]);
title('AKmaxWA Seafloor Deformation')
fname = 'AKmaxWA.png'
savefig(fname)
print('Created ',fname);

Created  AKmaxWA.png
In [ ]: