This Jupyter notebook can be found in collection of Clawpack apps as the file
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.topotools module. These tools facilitate creating and manipulating the topo files that are required as GeoClaw input to specify topography and bathymetry (underwater topography).
See http://www.clawpack.org/topotools_module.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.
Populating the interactive namespace from numpy and matplotlib
from __future__ import print_function from clawpack.geoclaw import topotools import os CLAW = os.environ['CLAW'] datadir = os.path.join(CLAW,'geoclaw','scratch') # directory for some sample data files
Many GeoClaw examples are set up to use topography files that have already been created and archived on the web, e.g. the example found in
$CLAW/geoclaw/examples/tsunami/chile2010 uses a topo file that can be obtained by these commands (which are also found in the Python script
maketopo.py found in that directory):
from clawpack.clawutil.data import get_remote_file filename = 'etopo10min120W60W60S0S.asc' url = 'http://www.geoclaw.org/topo/etopo/' + filename get_remote_file(url=url, output_dir=datadir, force=True, verbose=True)
File already exists, not downloading
If all you want to do is use this topo file in a GeoClaw run, you do not need to use any further Python tools.
topo_path = os.path.join(datadir, filename) topo = topotools.Topography() topo.read(topo_path, topo_type=2)
We can now do various things with the data. First let's print out some basic attributes:
print("The extent of the data in longitude and latitude: ") print(topo.extent)
The extent of the data in longitude and latitude: [-120.0, -59.999999999879996, -60.0, 1.2000356264252332e-10]
print("The shapes of 1-dimensional arrays of longitude x and latitude y values:", topo.x.shape, topo.y.shape) print("The shapes of 2-dimensional arrays X,Y and the topography data Z:", topo.Z.shape)
The shapes of 1-dimensional arrays of longitude x and latitude y values: (361,) (361,) The shapes of 2-dimensional arrays X,Y and the topography data Z: (361, 361)
From the filename you might guess this is 10-arcminute data, we can check that it is:
print("topo.delta = ",topo.delta) print("10 arcminutes is 1/6 degree = %8.6f degree" % (1./6.))
topo.delta = (0.166666666667, 0.166666666667) 10 arcminutes is 1/6 degree = 0.166667 degree
A simple plot with various defaults used can be obtained simply by:
<matplotlib.axes._subplots.AxesSubplot at 0x106157c50>
You can specify the limits of the color range:
<matplotlib.axes._subplots.AxesSubplot at 0x103901898>
Of course you can plot the data any way you want using the X, Y, and Z attributes.
Here are a couple examples:
figure(figsize=(12,6)) subplot(121) contourf(topo.X, topo.Y, topo.Z, [-7000, -5000, -4000, -3000, -2000, 0]) colorbar() # rescale aspect ratio based on mean latitude so distances more correct: gca().set_aspect(1.0 / cos(pi / 180.0 * topo.y.mean())) subplot(122) contourf(topo.X, topo.Y, topo.Z, [0,10000], colors = '#55ff55') contour(topo.X, topo.Y, topo.Z, , linestyles='-', colors='k') gca().set_aspect(1.0 / cos(pi / 180.0 * topo.y.mean()))
Sometimes it is useful to crop a large topo file to create a smaller one that contains a sub-region, e.g.
filter_region = (-90, -60, -60, -40) topo2 = topo.crop(filter_region) topo2.Z.shape
<matplotlib.axes._subplots.AxesSubplot at 0x10395e198>
The cropped topography can be written out to a new file, which will be smaller than the original:
topo2_path = os.path.join(datadir,'tip_of_south_america.tt3') topo2.write(topo2_path, topo_type=3)
print("Size of the new file: %s bytes" % os.stat(topo2_path).st_size) print("Size of the original file: %s bytes" % os.stat(topo_path).st_size)
Size of the new file: 345984 bytes Size of the original file: 985529 bytes
To obtain topography for other regions, see the links at http://clawpack.github.io/doc/tsunamidata.html
crop function also has an optional argument
coarsen with default value 1. If the value is larger integer than it coarsens the grid by that factor in each direction.
In the example below
filter_region=None specifies that the original domain size should be retained, but one can also crop it at the same time as coarsening.
topo3 = topo.crop(filter_region=None, coarsen=4) print("The original grid was of shape %s, topo3 has shape %s" % (topo.Z.shape,topo3.Z.shape)) figure(figsize=(12,6)) ax1 = subplot(121) topo.plot(axes=ax1) title('Original') ax2 = subplot(122) topo3.plot(axes=ax2) title('Coarsened')
The original grid was of shape (361, 361), topo3 has shape (91, 91)
<matplotlib.text.Text at 0x1071e5b38>
Topography class also have facilities for handling unstructured topography data. To start we will randomly sample the data we were working with above to get an unstructured data file. Here we generate random tuples of longitude-latitude and interpolate the original data so we can get a field at the random points.
import scipy.interpolate as interpolate N = 10000 x = numpy.random.rand(N) * (topo.x[-1] - topo.x) + topo.x y = numpy.random.rand(N) * (topo.y[-1] - topo.y) + topo.y # Here we need to transpose the Z array so that the coordinate evaluation is correct bathy = interpolate.RectBivariateSpline(topo.x, topo.y, topo.Z.transpose()) # strange hack to evaluate the function `bathy` produced by RectBivariateSpline since bathy(x, y) does not work z = numpy.empty(x.shape) for i in range(x.shape): z[i] = bathy(x[i], y[i])
Now lets construct the unstructured topography object based on the random locations we just generated.
unstructured_topo = topotools.Topography(unstructured=True) unstructured_topo.x = x unstructured_topo.y = y unstructured_topo.z = z
We can plot the data as is to see what we randomly created so far
<matplotlib.axes._subplots.AxesSubplot at 0x107273908>
One of the primary things we want to be able to do with unstructured topography data is to turn it into structured data which GeoClaw can use. In order to do this, there is a method in the
Topography class called
interp_unstructured which has a number of options. Here we will just go with the defaults (please see the documentation for the full assortment of options). In this case we will provide two things to the routine. The first is topography to use to fill in possible gaps in the unstructured topography. This "fill topography" can be structured or unstructured. The second argument we provide is a limit on the calculated grid spacing. The routine attempts to calculate the minimum distance between two points and uses this as the basis of the structured grid. Placing a limit on this insures that the resulting grid is not too fine. Addtional options for filtering the data in the routine, smoothing and interpolation strategies are also available.
#unstructured_topo.interp_unstructured(topo, delta_limit=10000) # fills in with original topography, if available unstructured_topo.interp_unstructured(, delta_limit=10000) # does not fill in with other topo, more realistic
Now that the topography has been interpolated onto a new grid, we can check to see if the
unstructured attribute is in fact now set to
A new set of attributes
X, Y, Z have been created, interpolated to a regular grid:
print("unstructured_topo.z has shape %s, unstructured_topo.Z has shape %s" \ % (unstructured_topo.z.shape, unstructured_topo.Z.shape))
unstructured_topo.z has shape (10000,), unstructured_topo.Z has shape (578, 578)
Usually at this point you will want to save out a copy of the interpolated topography for use with GeoClaw, the usual functions from before work of course. More interesting is to plot the new topography which we can again do by using the
<matplotlib.axes._subplots.AxesSubplot at 0x108d207b8>