topotools module for working with topography data

Warning

Many changes are being implemented in the way topo and dtopo files are handled, in both the Python tools and the Fortran code. See Changes to topo and dtopo handling (planned for v5.15.0) for a summary.

See also

Preprocessing attributes

Topography objects support seven preprocessing attributes that are applied automatically by read() in this order:

  1. negate_z — flip sign of Z (independent of topo_type < 0).

  2. z_shift — add a constant to all non-missing Z values.

  3. x_shift — add a constant to all x coordinates.

  4. crop_extent, buffer, align, coarsen — crop and subsample via crop().

Set them before calling read():

from clawpack.geoclaw.topotools import Topography

topo = Topography()
topo.crop_extent = [-100., -60., 10., 50.]
topo.coarsen = 2
topo.z_shift = 10.0
topo.read('bathymetry.nc', topo_type=4)

See Topography preprocessing attributes for a full attribute table with types and defaults, and non-obvious behavior notes.

Lazy-load pattern for NetCDF (read_header)

For NetCDF files (topo_type=4), read_header() uses TopoInspector to detect coordinate variable names via CF conventions (standard_name, axis, or common fallback names such as lon/lat, x/y). After calling read_header(), extent and delta are populated without loading the elevation array. Accessing Z triggers a deferred read():

topo = Topography()
topo.path = 'large_dem.nc'
topo.topo_type = 4
topo.read_header()     # fast: reads coordinate arrays only
print(topo.extent)     # available immediately
print(topo.delta)      # available immediately
z = topo.Z             # triggers full read on first access

This pattern is also used internally by TopographyData._compute_priority_order to determine file resolution without loading elevation data.

Vertical datum metadata

Topography and DTopography carry an optional datum attribute recording the vertical reference of the elevation/deformation data (e.g. 'MSL' or 'NAVD88'). It defaults to None, is populated from a NetCDF vertical_datum (or datum) attribute on read, and is written back out on NetCDF write (topo_type=4 / dtopo_type=4). ASCII formats have no datum field, so the value is not persisted for types 1/2/3.

The datum is informational only – GeoClaw performs no vertical-datum transformation and always uses the Z values as given. As a guard, both write() (when producing topo.data) and write() (dtopo.data) issue a warning if the files they are given carry more than one distinct datum, since mixing vertical references without converting between them is a likely error.

Deprecated since version ``topo_type=1``: (three-column x y z ASCII, one point per line) is deprecated. Reading emits a DeprecationWarning; writing also emits a DeprecationWarning; setting any preprocessing attribute before reading raises NotImplementedError.

To convert a type-1 file:

topo = Topography()
topo.read('old.tt1', topo_type=1)   # DeprecationWarning
topo.write('new.tt2', topo_type=2)  # save as type 2

Genuinely unstructured (scattered) point data cannot be converted this way. Either grid it externally (e.g. scipy.interpolate, GMT) or use interp_unstructured() (see Gridding unstructured (scattered) data below).

Gridding unstructured (scattered) data

GeoClaw’s solvers require topography on a regular, logically rectangular grid (Cartesian x-y in meters, or lon-lat – see coordinate_system in Specifying GeoClaw parameters in setrun.py), but survey or sounding data often comes as scattered (x, y, z) points. A Topography constructed with unstructured=True holds such points in its x, y, and z arrays, and interp_unstructured() interpolates them onto a regular grid, optionally filling gaps from one or more coarser (structured or unstructured) “fill” topographies:

from clawpack.geoclaw.topotools import Topography

# Scattered survey points (x, y, z)
survey = Topography(unstructured=True)
survey.x = x_points
survey.y = y_points
survey.z = z_values

# A coarser regional DEM used to fill gaps between survey points
regional = Topography(path='regional.tt3', topo_type=3)
regional.read()

survey.interp_unstructured(regional, extent=[x1, x2, y1, y2],
                           proximity_radius=100.0)
survey.unstructured        # now False -- it holds a regular grid
survey.write('combined.tt3', topo_type=3)

The grid spacing is taken from the minimum spacing of the scattered points (bounded below by delta_limit meters) unless set explicitly with delta. A fill point is used only where it lies inside extent, carries valid data, and is more than proximity_radius meters from every scattered point, so the higher-resolution survey data is preferred wherever it exists. Missing fill values (NaN in memory, or a numeric no_data_value) are dropped. See interp_unstructured() for the full set of parameters.

The notebook topotools_examples illustrates how to use some of the tools. Needs to be updated.

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. Needs to be updated.

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?

  • Add more robust plotting capabilities

class clawpack.geoclaw.topotools.Topography(path=None, topo_type=None, topo_func=None, unstructured=False, **kwargs)

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()
property X

Two dimensional coordinate array in x direction.

property Y

Two dimensional coordinate array in y direction.

property Z

A representation of the data as a 2d array.

crop(filter_region=None, coarsen=1, buffer=0, align=None)

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

Input:

  • filter_region (tuple): (x1,x2,y1,y2) desired new extent

  • coarsen (int): coarsening factor (by subsampling)

  • buffer (int): when possible, have at least this many points

    outside the filter_region on each side

  • align (tuple): (xalign,yalign) = desired alignment if coarsening

Setting buffer > 0 may be useful to insure that the computational domain lies entirely inside a cropped topo file (In GeoClaw, cell-centered topo values B are computed by integrating topo file values that are viewed as pointwise values, so topo point values are needed out to the domain edges.)

When subsampling with coarsen > 1, the align parameter may be useful to insure that the subsampling starts at an appropriate index. For example, if the original topo has

topo.x = [0, 0.5, 1, 1.5, 2, 2.5]

then coarsening by 2 would result in

newtopo.x = [0, 1, 2] # if align[0] is an integer

or

newtopo.x = [0.5, 1.5, 2.5] # if align[0] is an integer + 0.5

Often in GeoClaw, if the original topofile is aligned with integer longitudes and latitudes, for example, then we want any subsampled topo to have the same property.

In general, it tries to choose a starting index so that

(newtopo.x[0] - align[0]) / dx_new

is an integer, where dx_new is the spacing of points in the new topo after coarsening. This may not be possible, since it depends on the alignment of the original topography, in which case it will choose the index for which the misalignment is minimized.

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.

property delta

Spacing of data points.

property 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)

Return a boolean mask of grid points inside polygon.

Input:

  • polygon - sequence of (x, y) vertices describing a closed polygon (the closing edge from the last vertex back to the first is implied).

Output:

  • mask (numpy.ndarray of bool) - True at points lying inside polygon, False elsewhere. For a structured grid the mask has the same shape as self.X / self.Y; for unstructured data it is 1-D over the scattered self.x / self.y points.

Uses matplotlib.path.Path for a robust point-in-polygon test: it handles concave polygons and is independent of vertex winding order.

Example – keep only topography inside a region of interest:

mask = topo.in_poly(region_vertices)
topo.Z[~mask] = numpy.nan
interp_unstructured(fill_topo, extent=None, method='nearest', delta=None, delta_limit=20.0, no_data_value=-99999, 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 -99999.

  • 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_function(interp_kwargs={})

Create a function of (x,y) that returns the topo Z interpolated to a point (or to a 1D transect or 2D grid of points).

Inputs:
interp_kwargs: dictionary of parameter values to be passed to

RegularGridInterpolator. See defaults below.

Outputs:

topo_func: The function created

See the docstring in topo_func below for details on what shapes its arguments (x,y) can be.

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, long_lat=True, fig_kwargs={}, data_break=0.0, cb_kwargs={})

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.

  • long_lat (bool) - If this is a longitude-latitude plot then set the aspect of the plot to compensate for stretching. If not then the aspect is set to “equal”.

  • 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.

  • cb_kwargs (dict) - keyword arguments to be passed to colorbar e.g. ‘shrink’, ‘extend’, ‘label’. Can also set ‘title’ for cbar

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) - options for NetCDF (topo_type=4) reading:

    • z_var (str): name of the elevation variable, if it cannot be auto-detected by CF standard_name or common names.

    • assume_units (str): unit to assume for the elevation variable when the file has no units attribute (e.g. “m”). Units are otherwise required and never silently assumed: a file whose elevation variable lacks units, or whose units are not meters, raises ValueError (GeoClaw does not convert on read; pre- convert non-meter data to meters first).

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(value=nan, method='fill')

Replace missing (NaN) cells in Z using method.

Missing data is represented in memory as numpy.nan (the numeric no_data_value is only the on-file sentinel). This locates those cells and replaces them via replace_values().

Input:

  • value (float) - constant used when method == 'value'.

  • method (str) - one of 'value', 'nearest', 'linear', or 'fill'; see replace_values().

replace_values(indices, value=nan, method='fill')

Replace the Z values at indices in place using method.

Input:

  • indices - sequence of (i, j) index pairs identifying the cells to replace, e.g. the output of numpy.argwhere(condition).

  • value (float) - constant used when method == 'value'. Default numpy.nan.

  • method (str) - how to choose replacement values:

    • 'value' - set the cells to the constant value.

    • 'nearest' - nearest-neighbor value from the remaining (non-replaced) cells.

    • 'linear' - linear interpolation from the remaining cells; cells outside the convex hull of the remaining data are left as numpy.nan.

    • 'fill' - replace each cell with the average of the nearest surrounding non-replaced cells, growing the search box until at least one is found (the default).

'nearest' and 'linear' interpolate in index space, which is equivalent to physical space for a regularly spaced grid.

set_xyZ(X, Y, Z)

Set _x, _y, and _Z attributes and then generate X,Y,Z.

If X,Y are 1d arrays, then shape of Z should be (len(Y), len(X)).

Allow X,Y to be 2d arrays of shape Z.shape, in which case first extract x,y

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, fill_value=None, header_style='geoclaw', Z_format='%15.7e', grid_registration=None, z_dtype='float32', compression=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

  • fill_value (float) - value to use if filling a masked array

  • 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.

  • z_dtype (str) - on-disk dtype of the elevation variable when

    writing NetCDF (topo_type=4). Default “float32”, which still gives sub-millimeter precision for Earth topography (abs(Z) < 10000 m) while halving file size; pass “float64” for full double precision. Ignored for ASCII topo types. The elevation is always written with a CF units = “m” attribute (GeoClaw requires meters; see topo_netcdf).

  • compression - NetCDF (topo_type=4) zlib compression for the

    elevation variable. None/False (default) writes uncompressed; True uses zlib level 1 + byte shuffle; an int selects the zlib complevel; a dict is passed through verbatim. The compressed file stays randomly readable and needs no reader change. See netcdf_utils.compression_encoding.

property x

One dimensional coorindate array in x direction.

property y

One dimensional coordinate array in y direction.

property 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.extract_datum(*attr_dicts)

Return the first recognized vertical-datum attribute, or None.

Searches each mapping in attr_dicts (e.g. a NetCDF variable’s attrs then the dataset’s global attrs) for any of _DATUM_ATTR_NAMES and returns the first value found. Informational only; GeoClaw applies no vertical-datum transformation.

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, buffer=0, align=None, 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

  • buffer (int): when possible, have at least this many points

    outside the filter_region on each side

  • align (tuple): (xalign,yalign) = desired alignment if coarsening

    See the doc string for Topography.crop()

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, buffer=1, align=(-126,46), verbose=True)

# results in topo.x = array([-126.03333333, -126., …])

# 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.

Note that the newer etopo 2022 30 arcsecond DEM can be sampled using path = ‘etopo22_30sec’, but this topo is aligned differently with e.g. x = -126. falling half way between points. Also note that Z values in the newer dataset are no longer integers.

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.