Illustrate topo function¶
The make_function method will be added to clawpack.geoclaw.Topography in v5.14.0.
In [1]:
%matplotlib inline
In [2]:
from pylab import *
from clawpack.geoclaw import topotools
Read etopo 2022 30 arcsec topo, coarsened by 2 so resolution is 1 arcminute
In [3]:
etopo22 = topotools.read_netcdf('etopo22_30sec', extent=[-125,-122,47,49], coarsen=2, buffer=1)
Specify points on a transect:
In [4]:
xtrans = linspace(-124.75, -122.3, 1000)
ytrans = linspace(47.25, 48.5, 1000)
In [5]:
etopo22.plot()
plot(xtrans, ytrans, 'k')
title('topography and transect')
Out[5]:
Text(0.5, 1.0, 'topography and transect')
Extract topo values on the transect and plot¶
First use the default linear interpolation:
In [6]:
topofunc = etopo22.make_function()
etopo22_ztrans = topofunc(xtrans, ytrans)
figure(figsize=(12,6))
plot(xtrans, etopo22_ztrans, 'g')
grid(True)
xlabel('longitude')
ylabel('topo elevation (m)')
title('topography on transect');
Now use method='nearest' in interpolator to better show the resolution of the topography.
Note that the transect may cut across corners of some topo DEM cells, centers of others, accounting for spikiness.
In [7]:
topofunc = etopo22.make_function(interp_kwargs={'method':'nearest'})
ztrans = topofunc(xtrans, ytrans)
figure(figsize=(12,6))
plot(xtrans, ztrans, 'g')
grid(True)
xlabel('longitude')
ylabel('topo elevation (m)')
title('topography on transect');
Compare to older etopo1 topography¶
The make_topo function is handy for comparing two different topo files.
Compare the etopo 2022 at 1 arcminute obtained above with the older etopo1 topography at this resolution:
In [8]:
etopo1 = topotools.read_netcdf('etopo1', extent=[-125,-122,47,49], coarsen=1, buffer=1)
In [9]:
etopo1func = etopo1.make_function()
etopo1_ztrans = etopo1func(xtrans, ytrans)
figure(figsize=(12,6))
plot(xtrans, etopo22_ztrans, 'g', label='etopo 2022 at 1 arcmin')
plot(xtrans, etopo1_ztrans, 'r', label='older etopo1')
legend()
grid(True)
xlabel('longitude')
ylabel('topo elevation (m)')
title('topography on transect');
Compare the nearshore topography on the outer coast:
In [10]:
figure(figsize=(12,6))
plot(xtrans, etopo22_ztrans, 'g', label='etopo 2022 at 1 arcmin')
plot(xtrans, etopo1_ztrans, 'r', label='older etopo1')
xlim(-124.75, -124)
ylim(-150,150)
legend()
grid(True)
xlabel('longitude')
ylabel('topo elevation (m)')
title('topography on transect');
In [ ]: