Flatten nearshore topography to minimize reflections

%matplotlib inline

from pylab import *
from clawpack.geoclaw import topotools, kmltools

topo = topotools.Topography('../etopo22_15s_-137_-121_37_55.asc',3)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[3], line 1
----> 1 topo = topotools.Topography('../etopo22_15s_-137_-121_37_55.asc',3)

File ~/git/clawpack/geoclaw/src/python/geoclaw/topotools.py:446, in Topography.__init__(self, path, topo_type, topo_func, unstructured, **kwargs)
    443 self.coordinate_transform = lambda x,y: (x,y)
    445 if path:
--> 446     self.read(path=path, topo_type=topo_type,
    447               unstructured=unstructured, **kwargs)

File ~/git/clawpack/geoclaw/src/python/geoclaw/topotools.py:706, in Topography.read(self, path, topo_type, unstructured, mask, filter_region, force, stride, nc_params)
    702     self._delta = (dx,dy)
    704 elif abs(self.topo_type) in [2,3]:
    705     # Get header information
--> 706     N = self.read_header()  # note this also sets self._extent
    707                             # self._x, self._y, self._delta,
    708                             # and  self.grid_registration
    710     if abs(self.topo_type) == 2:
    711         # Data is read in as a single column, reshape it

File ~/git/clawpack/geoclaw/src/python/geoclaw/topotools.py:838, in Topography.read_header(self)
    835 self._extent = [numpy.nan,numpy.nan,numpy.nan,numpy.nan]
    836 self._delta = numpy.nan
--> 838 with open(self.path, 'r') as topo_file:
    839     # Check to see if we need to flip the header values
    840     first_line = topo_file.readline()
    841     try:

FileNotFoundError: [Errno 2] No such file or directory: '../etopo22_15s_-137_-121_37_55.asc'

#topo.crop(coarsen=4).plot(limits=(-1000,1000))

extent = [-124.75, -123., 44.25, 47.25]
topo_coast = topo.crop(extent)

coast_xy = topo_coast.make_shoreline_xy()
save('etopo22_crop_coast.npy', coast_xy)

Replace topo above some cutoff with flat bottom¶

Bcutoff = -60.

fig,ax = subplots(figsize=(8,12))
topo_coast.plot(axes=ax,limits=(-150,150))
contour(topo_coast.X, topo_coast.Y, topo_coast.Z, 
        [-100, Bcutoff-1], colors='yellow', linestyles='-')
title(f'Topography with contours \nat -100m and {Bcutoff-1}m');
fname = 'topo_coast.png'
savefig(fname, bbox_inches='tight')
print('Created ',fname)

Z_nocoast = where(topo_coast.Z > Bcutoff, Bcutoff, topo_coast.Z)

topo_nocoast = topotools.Topography()
topo_nocoast.set_xyZ(topo_coast.x, topo_coast.y, Z_nocoast)

fig,ax = subplots(figsize=(8,12))
topo_nocoast.plot(axes=ax,limits=(-150,150))
contour(topo_nocoast.X, topo_nocoast.Y, topo_nocoast.Z, 
        [-100, Bcutoff-1], colors='yellow', linestyles='-')
# show original coastline for reference:
contour(topo_coast.X, topo_coast.Y, topo_coast.Z, 
        [0], colors='g', linestyles='-', linewidths=0.5)
title(f'Flattened Topography with contours \nat -100m and {Bcutoff-1}m');
fname = f'topo_flat{int(Bcutoff)}.png'
savefig(fname, bbox_inches='tight')
print('Created ',fname)

fname = f'etopo22_flat{int(Bcutoff)}.asc'
topo_nocoast.write(fname, topo_type=3, header_style='arcgis', Z_format='%.3f')
print('Created ',fname)

Version with vertical wall at depth Bcutoff:¶

Bcutoff = -10.
Bwall = 100.
Z_wall = where(topo_coast.Z > Bcutoff, Bwall, topo_coast.Z)

topo_wall = topotools.Topography()
topo_wall.set_xyZ(topo_coast.x, topo_coast.y, Z_wall)

fig,ax = subplots(figsize=(8,12))
topo_wall.plot(axes=ax,limits=(-150,150))
contour(topo_wall.X, topo_wall.Y, topo_wall.Z, 
        [-100, Bcutoff-1], colors='yellow', linestyles='-')
# show original coastline for reference:
contour(topo_coast.X, topo_coast.Y, topo_coast.Z, 
        [0], colors='g', linestyles='-', linewidths=0.5)
title(f'Coastal wall Topography with contours \nat -100m and {Bcutoff-1}m');

fname = f'etopo22_wall{int(Bcutoff)}.asc'
topo_wall.write(fname, topo_type=3, header_style='arcgis', Z_format='%.3f')
print('Created ',fname)

topo_nocoast