Make Input Files -- Test 1

For this example simple artificial topography is generated in order to illustrate various things.

Contents:

In [1]:
%matplotlib inline
In [2]:
from pylab import *
from scipy.interpolate import interp1d
In [3]:
sys.path.insert(0,'../../new_python')
import topotools, marching_front, plottools, dtopotools

Define ocean topography

This simple topography is piecewise linear in $x$ (longitude) with a continental shelf and beach, and constant in the $y$ (latitude) direction. It is placed at the equator so distances are roughly equal in $x$ and $y$, and also placed at longitude 0.

In [4]:
# Define piecewise linear function (unequally spaced):
xocean = array([-2,-1,-0.5,-0.1,0.1])
zocean = array([-3000,-3000,-100,-100,100])

# Interpolate to equally spaced grid for topofile:
xo = arange(-2,0.2,0.1)
yo = array([-2,2])
zfunc = interp1d(xocean,zocean,fill_value="extrapolate")
zo = zfunc(xo)

# Convert to 2d arrays:
Xo,Yo = meshgrid(xo,yo)
Zo = vstack((zo,zo))
In [5]:
figure(figsize=(12,5))
subplot(121)
contourf(Xo,Yo,Zo)
colorbar()
title('Ocean Topography')
subplot(122)
plot(xo,zo,'k-')
fill_between(xo,zo,maximum(zo,0),color=[.5,.5,1])
title('Topography on transect')
Out[5]:
Text(0.5, 1.0, 'Topography on transect')

Save as a topofile:

In [6]:
topo = topotools.Topography()
topo.set_xyZ(xo,yo,Zo)

topo.write('input_files/topo_ocean.tt3', topo_type=3, Z_format="%11.3e")

Define topo for small coastal region

We define some more complicated topography on a finer grid over a small coastal region with 1/3 arcsecond resolution, chosen to be aligned with integer multiples of degrees (e.g. a grid point at longitude x=0 and latitude y=0) as typical of real DEMs from NCEI. This is important when aligning computational grids and fgmax grids (if used) in setrun.py.

We will use a cutoff function so that this fine-scale topo matches the linear beach profile of the ocean topography along the edges of this rectangle. The cutoff is 1 in the center of the rectangle and decays to 0 at the edges:

In [7]:
# choose DEM grid points:
arcsec13 = 1./(3*3600.)  # 1/3 arcsecond
print('arcsec13 = %.6f degrees = %.2f meters' % (arcsec13,arcsec13*111e3))
x = arange(-100*arcsec13, 150*arcsec13, arcsec13)
y = arange(-55*arcsec13, 55*arcsec13, arcsec13)
X,Y = meshgrid(x,y)
print('X.shape = ', X.shape)

x1,x2 = x.min(), x.max()
y1,y2 = y.min(), y.max()
print('Extent of coastal topo: (%.6f, %.6f, %.6f, %.6f)' % (x1,x2,y1,y2))

# define and plot the cutoff function:

w = 0.001 # width of cutoff layer
cutoff = 1. / (1. + exp(1e4*(X-(x2-w))) + exp(1e4*((x1+w)-X)) \
                  + exp(1e4*(Y-(y2-w))) + exp(1e4*((y1+w)-Y)))

figure(figsize=(10,6))
contourf(X,Y,cutoff)
colorbar(shrink=0.5)
gca().set_aspect(1)
title('Cutoff function');

arcsec13 = 0.000093 degrees = 10.28 meters
X.shape =  (110, 250)
Extent of coastal topo: (-0.009259, 0.013796, -0.005093, 0.005000)

The topography in this region is the linearly sloping beach augmented by a Gaussian dip. The beach slope is chosen to agree with the ocean topography offshore (1 km / degree, about 1/100), while onshore there is a smaller slope in this region for illustration.

In [8]:
Z0 = 1e3*X  # sloping beach matching ocean topography
Z1 = where(X<0, 1e3*X, 0.2e3*X)  # smaller slope on shore
R1 = (X-0.004)**2 + (Y-0.002)**2
Z1 += -4*exp(-500000*R1)         # Gaussian dip
Z = (1-cutoff)*Z0 + cutoff*Z1

Plot the coastal topography:

In [9]:
# colors:
c = [[.2,.2,1],[.5,.5,1],[.8,.8,1],[.7,1,.7],[.2,.8,0],[.9,.8,.2]]

figure(figsize=(12,7))
subplot(211)
contourf(X,Y,Z,[-2,-1,0,1,2],colors=c,extend='both')
cb = colorbar(shrink=0.9)
cb.set_label('meters')
contour(X,Y,Z,[-2,-1,0,1,2],colors=['b','b','r','g','g'])
gca().set_aspect(1.)
xticks(rotation=20)
xlabel('Longitude')
ylabel('Latitude')

subplot(212)
contourf(X*111e3,Y*111e3,Z,[-2,-1,0,1,2],colors=c,extend='both')
cb = colorbar(shrink=0.9)
cb.set_label('meters')
contour(X*111e3,Y*111e3,Z,[-2,-1,0,1,2],colors=['b','b','r','g','g'])
gca().set_aspect(1.)
xticks(rotation=20)
xlabel('meters')
ylabel('meters')

tight_layout();

The lower plot in the figure above shows the same topography as on the top, but with x,y units of meters to better show the scale. Recall that 1 degree is about 111 km and 1/3 arcsec is about 10 meters.

In the plots above, the red contour is at $Z = 0$, and hence is the "shoreline". However, the isolated "lake" with elevation $Z < 0$ could be dry land below sea level. Normally with GeoClaw this region would be filled with water initially up to $Z = 0$ everywhere. Below in the Force_Dry section, we discuss how to force this region to be initialized as dry if it is in fact dry land.

Save this as a topofile:

In [10]:
topo = topotools.Topography()
topo.set_xyZ(x,y,Z)

topo.write('input_files/topo_shore.tt3', topo_type=3, Z_format="%11.3e")
#topo.plot()

Plot both topo sets together

The coastal region above is very small compared to the ocean region defined above. Here we plot both together:

In [11]:
def plot_topo(add_colorbar=False):
    contourf(Xo,Yo,Zo,[-2,-1,0,1,2],colors=c,extend='both')
    contourf(X,Y,Z,[-2,-1,0,1,2],colors=c,extend='both')
    if add_colorbar: 
        cb = colorbar()
        cb.set_label('meters')
    #contour(X,Y,Z,[-2,-1,0,1,2],colors=['b','b','r','g','g'])
    plot([x1,x1,x2,x2,x1],[y1,y2,y2,y1,y1],'k')
    gca().set_aspect(1.)
    xticks(rotation=20)
    xlabel('Longitude')
    ylabel('Latitude')

figure(figsize=(12,6))
subplot(121)
plot_topo()
title('Ocean Topography')
subplot(122)
plot_topo(add_colorbar=True)
axis([-0.005, 0.015, -0.01, 0.01])
title('Zoom around shore')
tight_layout()

In the plot on the left above, the black rectangle showing the extent of the coastal DEM is barely visible. Zooming in shows that the topography does match up near the edges of this rectangle. In GeoClaw the finest available topography is used when computing cell-averaged topo values, so the coastal DEM will be used for any cell that overlaps this region.

Create dtopo for an earthquake source:

We define a simple earthquake in which there is uniform slip on a single subfault. The parameters are chosen to be somewhat reasonable for a subduction zone event offshore, but the shape is a bit odd (width 100 km and length 50 km) in order to give a smallish event with the desired onshore subsidence, for illustration purposes.

In [12]:
subfault = dtopotools.SubFault()
subfault.strike = 0.
subfault.length = 50.e3
subfault.width = 100.e3
subfault.depth = 10.e3
subfault.slip = 5.
subfault.rake = 90.
subfault.dip = 10.
subfault.longitude = -1.
subfault.latitude = 0.
subfault.coordinate_specification = "top center"

fault = dtopotools.Fault()
fault.subfaults = [subfault]

print("Earthquake magnitude: Mw = %.2f" % fault.Mw())
dtopo_fname = 'input_files/dtopo_test.tt3'
print("Using Okada model to create dtopo file", dtopo_fname)

x_deform = linspace(-2, 1, 100)
y_deform = linspace(-1, 1, 100)
times = [1.]

fault.create_dtopography(x_deform,y_deform,times)
dtopo = fault.dtopo

dtopo.write(dtopo_fname, dtopo_type=3)

figure(figsize=(12,6))
ax = subplot(121)
dtopo.plot_dZ_colors(2.,axes=ax,dZ_interval=0.5)
contour(Xo,Yo,Zo,[-110,-90,0],colors=['b','b','r'],linestyles='--')
ax.set_aspect(1.)
axis([-2,0.5,-2,2])
xlabel('Longitude')
ylabel('Latitude')

ax = subplot(122)
ylat = 0.
jlat = where(dtopo.y<=ylat)[0].max()
plot(dtopo.x, dtopo.dZ[0,jlat,:],'g')
plot(dtopo.x, 0*dtopo.x, 'k')
xlabel('Longitude')
title('Vertical displacement on transect at latitude %.2f' % ylat);

Earthquake magnitude: Mw = 7.97
Using Okada model to create dtopo file input_files/dtopo_test.tt3

The left plot above shows the sea floor deformation as contours and colors, along with the extent of the continental shelf as blue dashed lines and the shoreline as a red dashed line. The plot on the right shows the vertical deformation along a transect at latitude 0 going through the coastal region of interest.

We can compute the subsidence at the location on the shoreline where our fine scale topography is defined as:

In [13]:
xlon = 0.
ilon = where(dtopo.x<=xlon)[0].max()
ylat = 0.
jlat = where(dtopo.y<=ylat)[0].max()
#print(ilon,jlat)
dz0 = dtopo.dZ[0,jlat,ilon]
print('Surface deformation at x=%.2f, y=%.2f is dz = %.2f meters' \
      % (xlon,ylat,dz0))

Surface deformation at x=0.00, y=0.00 is dz = -0.68 meters

This subsidence is enough to significantly change the shoreline location, as seen below:

In [14]:
figure(figsize=(12,6))
subplot(211)
contourf(X,Y,Z,[-2,-1,0,1,2],colors=c,extend='both')
cb = colorbar(shrink=0.9)
cb.set_label('meters')
contour(X,Y,Z,[-2,-1,0,1,2],colors=['b','b','r','g','g'])
gca().set_aspect(1.)
xticks(rotation=20)
#xlim(-0.002,0.008)
xlabel('Longitude')
ylabel('Latitude')
title('Original topo')

subplot(212)
Z_postquake = Z + dz0
contourf(X,Y,Z_postquake,[-2,-1,0,1,2],colors=c,extend='both')
cb = colorbar(shrink=0.9)
cb.set_label('meters')
contour(X,Y,Z_postquake,[-2,-1,0,1,2],colors=['b','b','r','g','g'])
gca().set_aspect(1.)
xticks(rotation=20)
#xlim(-0.002,0.008)
xlabel('Longitude')
ylabel('Latitude')
title('Subsided topo, dz = %.2f m' % dz0);

tight_layout()
savefig('topo_with_dz.png')

Force Dry array

Now suppose that the onshore lake shown in the plots above is really a depression that should be dry land in spite of being below sea level. We can use the marching front algorithm described in the notebook MarchingFront.html to identify points that are below sea level but disconnected from the coast.

We use the marching front algorithm starting by assuming any point with Z < Z1 = -5 meters should be wet and marching to find all connected points with elevation up to Z = Z2 = 0:

In [15]:
wet_points = marching_front.select_by_flooding(topo.Z, Z1=-5., Z2=0., max_iters=None)

Selecting points with Z1 = -5, Z2 = 0, max_iters=27500
Done after 55 iterations with 11110 points chosen

See the notebook ForceDry.html for more discussion of the cells below...

In [16]:
Zdry = ma.masked_array(topo.Z, wet_points)
Zwet = ma.masked_array(topo.Z, logical_not(wet_points))

figure(figsize=(12,6))
subplot(211)
contourf(X,Y,Zdry,[-2,-1,0,1,2],colors=c,extend='both')
cb = colorbar(shrink=0.9)
cb.set_label('meters')
contour(X,Y,Z,[-2,-1,0,1,2],colors='k',linewidths=0.8)
gca().set_aspect(1.)
xticks(rotation=20)
#xlim(-0.002,0.008)
xlabel('Longitude')
ylabel('Latitude')
title('Colored points are identified as initially dry');

subplot(212)
contourf(X,Y,Zwet,[-2,-1,0,1,2],colors=c,extend='both')
cb = colorbar(shrink=0.9)
cb.set_label('meters')
contour(X,Y,Z,[-2,-1,0,1,2],colors='k',linewidths=0.8)
gca().set_aspect(1.)
xticks(rotation=20)
#xlim(-0.002,0.008)
xlabel('Longitude')
ylabel('Latitude')
title('Colored points are identified as initially wet');
tight_layout()

Create force_dry_init array for GeoClaw

First we buffer the points identified above as discussed in the ForceDry.html notebook.

In [17]:
dry_points = 1 - wet_points
dry_points_sum = dry_points[1:-1,1:-1] + dry_points[0:-2,1:-1] + dry_points[2:,1:-1] + \
                 dry_points[1:-1,0:-2] + dry_points[0:-2,0:-2] + dry_points[2:,0:-2] + \
                 dry_points[1:-1,2:] + dry_points[0:-2,2:] + dry_points[2:,2:]
        
# initialize array to 0 everywhere:
force_dry_init = zeros(dry_points.shape)
# reset in interior to 1 if all points in the 3x3 block around it are dry:
force_dry_init[1:-1,1:-1] = where(dry_points_sum == 9, 1, 0)

And finally create the input file needed for GeoClaw:

In [18]:
force_dry_init_topo = topotools.Topography()
force_dry_init_topo.set_xyZ(topo.x,topo.y,force_dry_init)

fname_force_dry_init = 'input_files/force_dry_init.data'
force_dry_init_topo.write(fname_force_dry_init, topo_type=3, Z_format='%1i')
print('Created %s' % fname_force_dry_init)

Created input_files/force_dry_init.data

See RunGeoclaw.html for more discussion and sample GeoClaw results.