topotools module for working with topography data

Warning

This describes new tools added in Clawpack 5.2.1

geoclaw/topotools_examples.ipynb illustrates how to use some of the tools.

The file $CLAW/geoclaw/tests/test_topotools.py contains some tests of these tools. Looking at these test routines may also give some ideas on how to use them.

Documentation auto-generated from the module docstrings

GeoClaw topotools Module $CLAW/geoclaw/src/python/geoclaw/topotools.py

Module provides several functions for reading, writing and manipulating topography (bathymetry) files.

Classes:
  • Topography
Functions:
  • determine_topo_type
  • create_topo_func
  • topo1writer
  • topo2writer
  • topo3writer
  • swapheader
TODO:
  • Add sub and super sampling capababilities
  • Add functions for creating topography based off a topo function, incorporate the create_topo_func into Topography class, maybe allow more broad initialization ability to the class to handle this?
  • Fix in_poly function
  • Add remove/fill no data value
  • Add more robust plotting capabilities
class clawpack.geoclaw.topotools.Topography(path=None, topo_func=None, topo_type=None, unstructured=False)

Base topography class.

A class representing a single topography file.

Properties:

Note: Modified to check the grid_registration when reading or writing topo files and properly deal with llcorner registration in which case the x,y data should be offset by dx/2, dy/2 from the lower left corner specified in the header of a DEM file.

Initialization:
Examples:
>>> import clawpack.geoclaw.topotools as topo
>>> topo_file = topo.Topography()
>>> topo_file.read('./topo.tt3', topo_type=3)
>>> topo_file.plot()
X

Two dimensional coordinate array in x direction.

Y

Two dimensional coordinate array in y direction.

Z

A representation of the data as a 2d array.

crop(filter_region=None, coarsen=1)

Crop region to filter_region

Create a new Topography object that is identical to this one but cropped to the region specified by filter_region

TODO:
  • Currently this does not work for unstructured data, could in principle
  • This could be a special case of in_poly although that routine could leave the resulting topography as unstructured effectively.
delta

Spacing of data points.

extent

Extent of the topography.

generate_2d_coordinates(mask=False)

Generate 2d coordinate arrays.

generate_2d_topo(mask=False)

Generate a 2d array of the topo.

in_poly(polygon)

Mask points (x,y) that are not in the specified polygon.

Uses simple ray casting algorithm for speed so beware of corner cases!

Input:
  • polygon (list) List of points that comprise the polygon. Note that order of the points will effect if this works (positive versus negative winding order). Points should be in counter-clockwise arrangement.
Returns:
  • X_mask (numpy.ma.MaskedArray) Masked array of X coordinates where those points outside of the polygon have been masked.
  • Y (numpy.ndarray) Coordinates in y direction in a meshgrid type of configuration.
interp_unstructured(fill_topo, extent=None, method='nearest', delta=None, delta_limit=20.0, no_data_value=-9999, buffer_length=100.0, proximity_radius=100.0, resolution_limit=2000)

Interpolate unstructured data on to regular grid.

Function to interpolate the unstructured data in the topo object onto a structured grid. Utilizes a bounding box plus a buffer of size buffer_length (meters) containing all data unless extent is not None is True. Then uses the fill topography fill_topo to fill in the gaps in the unstructured data. By default this is done by masking the fill data with the extents, the value no_data_value and if proximity_radius (meters) is not 0, by a radius of proximity_radius from all grid points in the object. Stores the result in the self.X, self.Y and self.Z object attributes. The resolution of the final grid is determined by calculating the minimum distance between all x and y data with a hard lower limit of delta_limit (meters).

Note that the function scipy.interpolate.griddata does not respect masks so a call to numpy.ma.MaskedArray.compressed() must be made to remove the masked data.

Input:
  • fill_topo (list) - List of Topography objects to use as fill data in the projection.
  • extent (tuple) - A tuple defining the rectangle of the sub-section. Must be in the form (x lower,x upper,y lower, y upper).
  • method (string) - Method used for interpolation, valid methods are found in scipy.interpolate.griddata. Default is nearest.
  • delta (tuple) - Directly set the grid spacing of the interpolation rather than determining it from the data itself. Defaults to None which causes the method to determine this value itself. Should be a 2-tuple of floats (delta_x, delta_y).
  • delta_limit (float) - Limit of finest horizontal resolution, default is 20 meters.
  • no_data_value (float) - Value to use if no data was found to fill in a missing value, ignored if method = ‘nearest’. Default is -9999.
  • buffer_length (float) - Buffer around bounding box, only applicable when extent is None. Default is 100.0 meters.
  • proximity_radius (float) - Radius every unstructured data point used to mask the fill data with. Default is 100.0 meters.
  • resolution_limit (int) - Limit the number of grid points in a single dimension. Raises a ValueError if the limit is violated. Default value is 2000.

Sets this object’s unstructured attribute to False if successful.

make_shoreline_xy(sea_level=0)

Returns an array shoreline_xy with 2 columns containing x and y values for all segements of the shoreline (defined to be the contour where self.z = sea_level) separated by [nan,nan] pairs. This allows all shorelines to be quickly plotted via:

>>> plot(shoreline_xy[:,0], shoreline_xy[:,1])

The shoreline can be saved as a binary .npy file via:

>>> numpy.save(filename, shoreline_xy)

which is much smaller than the original topography file. Reload via:

>>> shoreline_xy = numpy.load(filename)
plot(axes=None, contour_levels=None, contour_kwargs={}, limits=None, cmap=None, add_colorbar=True, plot_box=False, fig_kwargs={}, data_break=0.0)

Plot the topography.

Input:
  • axes (matplotlib.pyplot.axes) - If passed in, plot will be added to this axes. Otherwise a new plot figure will be created (using fig_kwargs) and a new axes object created and returned.
  • contour_levels (list) - levels for contour lines if these are to be added (default None). Set to [0.] to plot shoreline.
  • contour_kwargs (dict) - keyword arguments to be passed to contour command, e.g. {‘colors’:’r’, ‘linestyles’: ‘-‘}. Default is empty dict.
  • limits (list) - (min, max) of topo values for color map. Defaults to None, in which case (self.Z.min(), self.Z.max()) used.
  • cmap (matplotlib.colors.Colormap) - colormap, defaults to specialized map with blues for bathymetry and green/browns for topo.
  • fig_kwargs (dict) - keyword arguments to be passed to figure.
  • plot_box (bool or color specifier) - If evaluates to True, plot a box around limits of this topo.
  • data_break (float) - when default cmap is used, the value to use to break between water and land colormaps. Defaults to 0., but for some topo files may need to use e.g. 0.01 Or may want to show plots at different tide stage.
Output:
  • axes (matplotlib.pyplot.axes) - the axes on which plot created.
Note that:
  • if type(self.Z) is numpy.ma.MaskedArray then pcolor is used,
  • if type(self.Z) is numpy.ndarray then imshow is used. (This is faster for large files)
read(path=None, topo_type=None, unstructured=False, mask=False, filter_region=None, force=False, stride=[1, 1], nc_params={})

Read in the data from the object’s path attribute.

Stores the resulting data in one of the sets of x, y, and z or X, Y, and Z.

Input:
  • path (str) file to read
  • topo_type (int) - GeoClaw format topo_type
  • unstructured (bool) - default is False for lat-long grids.
  • mask (bool) - whether to store as masked array for missing values (default if False)
  • filter_region (tuple)
  • stride (list) - List of strides for the x and y dimensions respectively. Default is [1, 1]. Note that this is only implemented for NetCDF reading currently.
  • nc_params (dict) -

The first three might have already been set when instatiating object.

read_header()

Read in header of topography file at path.

If a value returns numpy.nan then the value was not retrievable. Note that this routine can read in headers whose values and labels are swapped.

replace_no_data_values(method='fill')

Replace no_data_value with other values as specified by method.

self.no_data_value

Input:

  • method can be one of:

    • fill - Fill in all no_data_value locations with value
    • nearest - Fill in no_data_value locations with average of nearest neighbors.
replace_values(indices, value=nan, method='fill')

Replace the values at indices by the specified method

Methods:
  • “fill”
  • “nearest”
smooth_data(indices, r=1)

Filter topo data at indices by averaging surrounding data.

Surrounding data is considered within the ball of radius r in the inf-norm. Acts as a low-band pass filter and removes oscillatory data.

Input:
  • indices (list)
  • r (int)
Output:

None
write(path, topo_type=None, no_data_value=None, masked=True, header_style='geoclaw', Z_format='%15.7e', grid_registration=None)

Write out a topography file to path of type topo_type.

Writes out a topography file of topo type specified with topo_type or inferred from the output file’s extension, defaulting to 3, to path from data in Z. The rest of the arguments are used to write the header data.

Input:
  • path (str) - file to write
  • topo_type (int) - GeoClaw format topo_type Note: this is second positional argument, agreeing with the read function in this class. It was the third argument in GeoClaw version 5.3.1 and earlier.
  • no_data_value - values used to indicate missing data
  • masked (bool) - unused??
  • header_style (str) - indicates format of header lines
    ‘geoclaw’ or ‘default’ ==> write value then label
    with grid_registration == ‘lower’ as default
    ‘arcgis’ or ‘asc’ ==> write label then value
    with grid_registration == ‘llcorner’ as default (needed for .asc files in ArcGIS)
  • Z_format (str) - string format to use for Z values The default format “%15.7e” gives at least millimeter precision for topography with abs(Z) < 10000 and results in smaller files than the previous default of “%22.15e” used in GeoClaw version 5.3.1 and earlier. A shorter format can be used if the user knows there are fewer significant digits, e.g. etopo1 data is integers and so has a resolution of 1 meter. In this case a cropped or coarsened version might be written with Z_format = “%7i”, for example.
  • grid_registration (str) - ‘lower’, ‘llcorner’, ‘llcenter’
    or None for defaults described above.
x

One dimensional coorindate array in x direction.

y

One dimensional coordinate array in y direction.

z

A representation of the data as an 1d array.

clawpack.geoclaw.topotools.create_topo_func(loc, verbose=False)

Given a 1-dimensional topography profile specfied by a set of (x,z) values, create a lambda function that when evaluated will give the topgraphy at the point (x,y). (The resulting function is constant in y.)

Example:
>>> f = create_topo_func(loc)
>>> b = f(x, y)
Input:

  • loc (list) - Create a topography file with the profile denoted by the tuples inside of loc. A sample set of points are shown below. Note that the first value of the list is the x location and the second is the height of the topography.

    z (m)
^                                                  o loc[5]  o
|                                                    
|                                          loc[4]   
|--------------------------------------------o-----> x (m) (sea level)
|                                            
|                                o loc[2] o loc[3]
|                         
|                         
|                           o loc[1]
|           
|                               
|__________________o loc[0]
0.0
clawpack.geoclaw.topotools.determine_topo_type(path, default=None)

Using the file suffix of path, attempt to deterimine the topo type.

Input:
  • path (string) - Path to the file. Can include archive extensions (they will be stripped off).
  • default (object) - Value returned if no suitable topo type was determined. Default is None.

returns integer between 1-3 or default if nothing matches.

clawpack.geoclaw.topotools.fetch_topo_url(url, local_fname=None, force=None, verbose=False, ask_user=False)

DEPRECATED: Use clawpack.clawutil.data.get_remote_file instead (see note below).

Replaces get_topo function.

Download a topo file from the web, provided the file does not already exist locally.

Input:
  • url (str) URL including file name
  • local_fname (str) name of local file to create. If local_fname == None, take file name from URL
  • force (bool) If False, prompt user before downloading.

For GeoClaw examples, some topo files can be found in `http://www.geoclaw.org/topo`_ See that website for a list of archived topo datasets.

If force==False then prompt the user to make sure it’s ok to download,

If force==None then check for environment variable CLAW_TOPO_DOWNLOAD and if this exists use its value. This is useful for the script python/run_examples.py that runs all examples so it won’t stop to prompt.

This routine has been deprecated in favor of clawpack.clawutil.data.get_remote_file. All the functionality should be the same but calls the other routine internally.

clawpack.geoclaw.topotools.get_topo(topo_fname, remote_directory, force=None)

DEPRECATED: Use clawpack.geoclaw.util.get_remote_file instead

Download a topo file from the web, provided the file does not already exist locally.

remote_directory should be a URL. For GeoClaw data it may be a subdirectory of http://www.clawpack.org/geoclaw/topo See that website for a list of archived topo datasets.

If force==False then prompt the user to make sure it’s ok to download, with option to first get small file of metadata.

If force==None then check for environment variable CLAW_TOPO_DOWNLOAD and if this exists use its value. This is useful for the script python/run_examples.py that runs all examples so it won’t stop to prompt.

clawpack.geoclaw.topotools.read_netcdf(path, zvar=None, extent='all', coarsen=1, return_topo=True, return_xarray=False, verbose=False)
Input:
  • path (str) - Path to the file to read, or url to remote file, or a key into the topotools.remote_topo_urls dictionary.
  • zvar (str) - variable to read as Z=elevation. if None, will try ‘Band1’, ‘z’, ‘elevation’.
  • extent - [x1,x2,y1,y2] for desired subset, or ‘all’ for entire file
  • coarsen (int) - factor to coarsen by, 1 by default.
  • return_topo (bool) - if True, return a topotools.Topography object. default is True
  • return_xarray (bool) - if True, return an xarray.Dataset object. default is False
Output:
  • topo and/or xarray_ds depending on what was requested. (either a single object or a tuple of two objects.)

If return_xarray == True then xarray is used to read the data, otherwise netCDF4 is used directly.

Sample usage:

from clawpack.geoclaw import topotools extent = [-126,-122,46,49] path = ‘etopo1’ topo = topotools.read_netcdf(path, extent=extent, coarsen=2, verbose=True)

# to plot: topo.plot()

# to save topofile for input to GeoClaw: topo.write(‘etopo_sample_2min.tt3’, topo_type=3, Z_format=’%.0f’)

This should give a 2-minute resolution DEM of the Western Washington coast. Note that etopo1 Z values are integers (vertical resolution is 1 meter) and using Z_format=’%.0f’ will save as integers to minimize file size.

clawpack.geoclaw.topotools.swapheader(inputfile, outputfile)

Swap the order of key and value in header to value first.

Note that this is a wrapper around functionality in the Topography class.

clawpack.geoclaw.topotools.topo1writer(outfile, topo, xlower, xupper, ylower, yupper, nxpoints, nypoints)

Function topo1writer will write out the topofiles by evaluating the function topo on the grid specified by the other parameters.

Assumes topo can be called on arrays X,Y produced by numpy.meshgrid.

Output file is of “topotype1,” which we use to refer to a file with (x,y,z) values on each line, progressing from upper left corner across rows, then down.

clawpack.geoclaw.topotools.topo2writer(outfile, topo, xlower, xupper, ylower, yupper, nxpoints, nypoints, nodata_value=-99999)

Write out a topo type 2 file by evaluating the function topo.

This routine is here for backwards compatibility and simply creates a new topography object and writes it out.

clawpack.geoclaw.topotools.topo3writer(outfile, topo, xlower, xupper, ylower, yupper, nxpoints, nypoints, nodata_value=-99999)

Write out a topo type 3 file by evaluating the function topo.

This routine is here for backwards compatibility and simply creates a new topography object and writes it out.