dtopotools module for moving topography

See also Earthquake sources: Fault slip and the Okada model for a discussion of the subfault parameters required for the Okada model of seafloor deformation resulting from slip on a specified fault plane.

Warning

Starting in Version 5.5.0, the subfault parameter rise_time now refers to the total rise time of a subfault, while rise_time_starting is the rise time up to the break in the piecewise quadratic function defining the rise. By default rise_time_ending is set equal to rise_time_starting. See the module function rise_fraction below for more description. (In earlier versions, rise_time read in from csv files, for example, was erroneously interpreted as rise_time_starting is now.)

The notebook dtopotools_examples.ipynb illustrates how to use some of the tools.

The notebook Okada.ipynb illustrates the Okada model using some of these tools.

The file $CLAW/geoclaw/tests/test_dtopotools.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 dtopotools Module $CLAW/geoclaw/src/python/geoclaw/dtopotools.py

Module provides several functions for dealing with changes to topography (usually due to earthquakes) including reading sub-fault specifications, writing out dtopo files, and calculating Okada based deformations.

Classes:

  • DTopography

  • SubFault

  • Fault

  • UCSBFault

  • CSVFault

  • SiftFault

  • SegmentedPlaneFault

  • Fault1d

  • SubFault1d

  • DTopography1d

Functions:

  • plot_dz_contours

  • plot_dz_colors

  • Mw

  • strike_direction

  • rise_fraction

class clawpack.geoclaw.dtopotools.CSVFault(subfaults=None, input_units={}, coordinate_specification=None)

Fault subclass for reading in CSV formatted files

Assumes that the first row gives the column headings

read(path, input_units={}, coordinate_specification='top center', rupture_type='static', verbose=False)

Read in subfault specification at path.

Creates a list of subfaults from the subfault specification file at path.

class clawpack.geoclaw.dtopotools.DTopography(path=None, dtopo_type=None)

Basic object representing moving topography

dZ_at_t(t)

Interpolate dZ to specified time t and return deformation.

dZ_max()

Return max(abs(dZ)) over all dz in self.dZ, the maximum surface deformation for this dtopo. DEPRECATE? – it’s now a 1-liner

plot_dZ_colors(t, axes=None, cmax_dZ=None, dZ_interval=None, colorbar_ticksize=10, colorbar_labelsize=10, fig_kwargs={})

Interpolate self.dZ to specified time t and then call module function plot_dZ_colors.

plot_dZ_contours(t, dZ_interval=0.5, axes=None, fig_kwargs={})

Interpolate self.dZ to specified time t and then call module function plot_dZ_contours.

read(path=None, dtopo_type=None, verbose=False)

Read in a dtopo file and use to set attributes of this object.

Input:

  • path (path) - Path to existing dtopo file to read in.

  • dtopo_type (int) - Type of topography file to read. Default is 3

    if not specified or apparent from file extension.

write(path=None, dtopo_type=None, dZ_format='%.3f')

Write out subfault resulting dtopo to file at path.

Input:

  • path (path) - Path to the output file to written to.

  • dtopo_type (int) - Type of topography file to write out. Default 3.

  • dZ_format (str) - format for dZ values printed. Default ‘%.3f’

class clawpack.geoclaw.dtopotools.DTopography1d(path=None, dtopo_type=None)

Basic object representing moving topography in 1d

dZ_at_t(t)

Interpolate dZ to specified time t and return deformation.

dZ_cellave_at_t(t, celledges)

Interpolate dZ to specified time t and return deformation.

read(path=None, dtopo_type=None, verbose=False)

Read in a dtopo file and use to set attributes of this object.

Input:

  • path (path) - Path to existing dtopo file to read in.

  • dtopo_type (int) - Type of topography file to read. Default is 3

    if not specified or apparent from file extension.

write(path=None, dtopo_type=None)

Write out subfault resulting dtopo to file at path.

Input:

  • path (path) - Path to the output file to written to.

  • dtopo_type (int) - Type of topography file to write out. Default is 3.

class clawpack.geoclaw.dtopotools.Fault(subfaults=None, input_units={}, coordinate_specification=None)

Base Fault class

A class describing a fault possibly composed of subfaults.

Properties:

Initialization:

Examples:

Mo()

Calculate the seismic moment for a fault composed of subfaults, in units N-m.

Mw()

Calculate the moment magnitude for a fault composed of subfaults.

containing_rect()

Find containing rectangle of fault in x-y plane.

Returns tuple of x-limits and y-limits.

create_dtopo_xy(rect=None, dx=0.016666666666666666, buffer_size=0.5)

Create coordinate arrays containing fault with a buffer.

Input:
  • rect - if None, use self.containing_rect

    Otherwise a list [x1,x2,y1,y2]

  • dx (int) - Spatial resolution. Defaults to 1” resolution.

  • buffer_size (float) - Buffer distance around edge of fault in degrees, defaults to 0.5 degrees.

Output:

  • x,y 1-dimensional arrays that cover the desired rect. They start at (x1,y1) and may go a bit beyond (x2,y2) depending on dx

create_dtopography(x, y, times=[0.0, 1.0], verbose=False)

Compute change in topography and construct a dtopography object.

Use subfaults’ okada routine and add all deformations together.

Raises a ValueError exception if the rupture_type is an unknown type.

returns a :class`DTopography` object.

plot_subfaults(axes=None, plot_centerline=False, slip_color=False, cmap_slip=None, cmin_slip=None, cmax_slip=None, slip_time=None, plot_rake=False, xylim=None, plot_box=True, colorbar_shrink=1, verbose=False, colorbar_labelsize=10, colorbar_ticksize=10)

Plot each subfault projected onto the surface.

axes can be passed in to specify the matplotlib.axes.AxesSubplot on which to add this plot. If axes == None, a new figure window will be opened. The axes on which it is plotted is the return value of this call.

If plot_centerline == True, plot a line from the centroid to the top center of each subfault to show what direction is up-dip.

If slip_color == True then use the color map cmap_slip (which defaults to matplotlib.cm.jet) to color the subplots based on the magnitude of slip, scaled between cmin_slip and cmax_slip. (If these are None then scaled automatically based on range of slip.) If slip_time == None then colors are based on the final slip. For dynamic faults, slip_time can be set to a time and the dynamic timing of each subfault will be used to compute and plot the slip at this time.

If plot_rake == True, plot a line from the centroid pointing in the direction of the rake (the direction in which the top block is moving relative to the lower block. The distance it moves is given by the slip.)

xylim can be set to a list or tuple of length 4 of the form [x1,x2,y1,y2] to specify the x- and y-axis limits.

If plot_box == True, a box will be drawn around each subfault.

plot_subfaults_depth(axes=None)

Plot the depth of each subfault vs. x and vs. y in a second plot.

read(path, column_map, coordinate_specification='centroid', rupture_type='static', skiprows=0, delimiter=None, input_units={}, defaults=None)

Read in subfault specification at path.

Creates a list of subfaults from the subfault specification file at path.

Inputs:

  • path (str) file to read in, should contain subfaults, one per line

  • column_map (dict) specifies mapping from parameter to the column of the input file that contains values for this parameter, e.g.

    column_map = {“latitude”:0, “longitude”:1, “depth”:2, “slip”:3, “rake”:4, “strike”:5, “dip”:6}

  • coordinate_specification (str) specifies the location on each subfault that corresponds to the (longitude,latitude) and depth of the subfault. See the documentation for SubFault.calculate_geometry.

  • rupture_type (str) either “static” or “kinematic”

  • skiprows (int) number of header lines to skip before data

  • delimiter (str) e.g. ‘,’ for csv files

  • input_units (dict) indicating units for length, width, slip, depth,

    and for rigidity mu as specified in file. These will be converted to “standard units”.

  • defaults (dict) default values for all subfaults, for values not

    included in subfault file on each line.

set_dynamic_slip(t)

Set slip_at_dynamic_t attribute of all subfaults to slip at the requested time t.

Input:

  • t (float) -

Raises a ValueError exception if this object’s rupture_type attribute is set to static.

write(path, style=None, column_list=None, output_units={}, delimiter='  ')

Write subfault format file with one line for each subfault. Can either specify a style that determines the columns, or a column_list. Must specify one but not both. See below for details.

Inputs:

  • path (str) file to write to.

  • style (str) to write in a style that matches standard styles adopted by various groups. One of the following:

    • “usgs” (Not implemented)

    • “noaa sift” (Not implemented)

    • “ucsb” (Not implemented)

  • column_list (list) specifies what order the parameters should be written in the output file, e.g.

    column_list = [‘longitude’,’latitude’,’length’,’width’, ‘depth’,’strike’,’rake’,’dip’,’slip’]

  • output_units (dict) specifies units to convert to before writing. Defaults to “standard units”.

  • delimiter (str) specifies delimiter between columns, e.g. “,” to create a csv file. Defaults to “ “.

class clawpack.geoclaw.dtopotools.Fault1d(subfaults=None, input_units={}, coordinate_specification=None)

Base Fault1d class

A class describing a 1d fault possibly composed of subfaults.

clawpack.geoclaw.dtopotools.Mw(Mo, units='N-m')

Calculate moment magnitude based on seismic moment Mo. Follows USGS recommended definition from

The SubFault and Fault classes each have a function Mo to compute the seismic moment for a single subfault or collection respectively.

class clawpack.geoclaw.dtopotools.SiftFault(sift_slip=None, longitude_shift=0.0)

Define a fault by specifying the slip on a subset of the SIFT unit sources. The database is read in by load_sift_unit_sources. See http://www.pmel.noaa.gov/pubs/PDF/gica2937/gica2937.pdf for a discussion of these unit sources, although the database used is more recent than what is reported in that paper and uses different notation for the subfault names. The subfault database used was downloaded from

Example:

>>> sift_slip = {'acsza1':2, 'acszb1':3}
>>> fault = SiftFault(sift_slip)

results in a fault with two specified subfaults with slip of 2 and 3 meters.

set_subfaults(sift_slip)
sift_slip (dict) is a dictionary with key = name of unit source

and value = magnitude of slip to assign (in meters).

class clawpack.geoclaw.dtopotools.SubFault

Basic sub-fault specification.

Locate fault plane in 3D space Note that the coodinate specification is in reference to the fault

Coordinates of Fault Plane:

The attributes centers and corners are described by the figure below.

centers[0,1,2] refer to the points labeled 0,1,2 below.

In particular the centroid is given by centers[1]. Each will be a tuple (x, y, depth).

corners[0,1,2,3] refer to the points labeled a,b,c,d resp. below.

Each will be a tuple (x, y, depth).

Top edge    Bottom edge
  a ----------- b          ^ 
  |             |          |         ^
  |             |          |         |
  |             |          |         | along-strike direction
  |             |          |         |
  0------1------2          | length  |
  |             |          |
  |             |          |
  |             |          |
  |             |          |
  d ----------- c          v
  <------------->
       width

  <-- up dip direction
Mo()

Calculate the seismic moment for a single subfault

Returns in units of N-m and assumes mu is in Pascals.

calculate_geometry()

Calculate the fault geometry.

Routine calculates the class attributes corners and centers which are the corners of the fault plane and points along the centerline respecitvely in 3D space.

Note: self.coordinate_specification specifies the location on each subfault that corresponds to the (longitude,latitude) and depth of the subfault. Currently must be one of these strings:

  • “bottom center”: (longitude,latitude) and depth at bottom center

  • “top center”: (longitude,latitude) and depth at top center

  • “centroid”: (longitude,latitude) and depth at centroid of plane

  • “noaa sift”: (longitude,latitude) at bottom center, depth at top,

    This mixed convention is used by the NOAA SIFT database and “unit sources”, see: http://nctr.pmel.noaa.gov/propagation-database.html

  • “top upstrike corner”: (longitude,latitude) and depth at

    corner of fault that is both updip and upstrike.

The Okada model is expressed assuming (longitude,latitude) and depth are at the bottom center of the fault plane, so values must be shifted or other specifications.

calculate_geometry_triangles()

Calculate geometry for triangular subfaults

  • Uses corners to calculate centers, longitude, latitude, depth, strike, dip, length, width.

  • sets coordinate_specification as “triangular”

  • Note that calculate_geometry() computes long/lat/strike/dip/length/width to calculate centers/corners

property centers

Coordinates along the center-line of the fault plane.

convert_to_standard_units(input_units, verbose=False)

Convert parameters from the units used for input into the standard units used in this module.

coordinate_specification

Specifies where the latitude, longitude and depth are measured from.

property corners

Coordinates of the corners of the fault plane.

depth

Depth of subfault based on coordinate_specification in meters.

dip

Subfault’s angle of dip

dynamic_slip(t)

For a dynamic fault, compute the slip at time t. Assumes the following attributes are set:

  • rupture_time

  • rise_time

  • rise_time_starting: optional, defaults to rise_time/2

  • rise_shape: optional, defaults to quadratic

property gauss_pts

Coordinates along the center-line of the fault plane.

latitude

Latitutde of the subfault based on coordinate_specification.

length

Length of subfault in meters.

longitude

Longitude of the subfault based on coordinate_specification.

mu

Rigidity (== shear modulus) in Pascals.

okada(x, y)

Apply Okada to this subfault and return a DTopography object.

Input:

  • x,y are 1d arrays

Output:
  • DTopography object with dZ array of shape (1,len(x),len(y))

    with single static displacement and times = [0.].

Currently only calculates the vertical displacement.

Okada model is a mapping from several fault parameters to a surface deformation. See Okada 1985 [Okada85], or Okada 1992, Bull. Seism. Soc. Am.

okadamap function riginally written in Python by Dave George for Clawpack 4.6 okada.py routine, with some routines adapted from fortran routines written by Xiaoming Wang.

Rewritten and made more flexible by Randy LeVeque

Note: self.coordinate_specification (str) specifies the location on each subfault that corresponds to the (longitude,latitude) and depth subfault.

See the documentation for SubFault.calculate_geometry for dicussion of the possible values self.coordinate_specification can take.

rake

Rake of subfault movement in degrees.

rise_shape

Shape of rise, ‘linear’ or ‘quadratic’ (piecewise)

rise_time

Total rise time of subfault (sec)

rise_time_starting

Optional time of first part of rise, before break in pw smooth rise

rupture_time

Time subfault starts to rupture (sec)

rupture_type

Either ‘static’ or ‘kinematic’

set_corners(corners, projection_zone=None)

Set three corners for a triangular fault. :Inputs:

  • corners should be iterable of length 3.

  • rake should be between -180. and 180.

slip

Slip on subfault in strike direction in meters.

strike

Strike direction of subfault in degrees.

width

Width of subfault in meters.

class clawpack.geoclaw.dtopotools.SubFault1d(strike=0, length=1000000.0)
coordinate_specification

location of x0: ‘top center’, ‘centroid’, or ‘bottom center’

latitude

latitude where placed on sphere to use 2d routines

length

length (m) in strike direction orthogonal to x

longitude

longitude at location specified by coordinate_specification

rake

rake=90 ==> top at right if strike=0, at left if strike=180

strike

rake=90, so top at right if strike=0, at left if strike=180

class clawpack.geoclaw.dtopotools.SubdividedPlaneFault(base_subfault, nstrike=1, ndip=1, slip_function=None, Mo=None)

Define a fault by starting with a single fault plane (specified as base_subfault of class SubFault) and subdividing the fault plane into a rectangular array of nstrike by ndip equally sized subfaults.

By default,the slip on each subfault will be initialized to base_subfault.slip so that the slip is uniform over the original plane and the seismic moment is independent of the number of subdivisions.

Alternatively, the slip distribution can be specified by providing a function slip_distribution, which should be a function of (xi,eta) with each variable ranging from 0 to 1. xi varies from 0 at the top of the fault to 1 at the bottom in the down-dip direction. eta varies from one edge of the fault to the other moving in the strike direction. This function will be evaluated at the centroid of each subfault to set the slip.

Can also specify a desired seismic moment Mo in which case the slips will be rescaled at the end so the total seismic moment is Mo. In this case the slip_distribution function only indicates the relative slip between subfaults.

subdivide(nstrike=1, ndip=1, slip_function=None, Mo=None)

Subdivide the fault plane into nstrike * ndip subfaults.

class clawpack.geoclaw.dtopotools.TensorProductFault(fault_plane, slip_along_strike=None, slip_down_dip=None, nstrike=1, ndip=1, Mo=None)

Define a fault by starting with a single fault plane (specified as fault_plane of class SubFault) and subdividing the fault plane into a rectangular array of nstrike by ndip equally sized subfaults.

Then define the slip on each subfault via two one-dimensional functions slip_along_strike and slip_down_dip that specify the slip as a function of fractional distance in the along-strike and down-dip direction respectively (i.e. the argument of each goes from 0 to 1).

Setting either to None defaults to constant function 1.

The slip is set by evaluating the tensor product at the centroid of each subfault.

Can specify a desired seismic moment Mo in which case the slips will be rescaled at the end.

class clawpack.geoclaw.dtopotools.UCSBFault(path=None, **kwargs)

Fault subclass for reading in subfault format models from UCSB

Read in subfault format models produced by Chen Ji’s group at UCSB, downloadable from:

read(path, rupture_type='static')

Read in subfault specification at path.

Creates a list of subfaults from the subfault specification file at path.

Subfault format contains info for dynamic rupture, so can specify rupture_type = ‘static’ or ‘kinematic’ (or ‘dynamic’ for backward compatibility).

clawpack.geoclaw.dtopotools.plot_dZ_colors(x, y, dZ, axes=None, cmax_dZ=None, dZ_interval=None, add_colorbar=True, verbose=False, colorbar_labelsize=10, colorbar_ticksize=10, fig_kwargs={})

Plot sea floor deformation dZ as colormap with contours

clawpack.geoclaw.dtopotools.plot_dZ_contours(x, y, dZ, axes=None, dZ_interval=0.5, verbose=False, fig_kwargs={})

For plotting seafloor deformation dZ

clawpack.geoclaw.dtopotools.rise_fraction(t, rupture_time, rise_time, rise_time_starting=None, rise_shape='quadratic')

A continuous piecewise quadratic or linear function of t that is

  • 0 for t <= rupture_time,

  • 1 for t >= rupture_time + rise_time

Inputs:

  • t (scalar, list, or np.array): times at which to evaluate the rise function.

  • rupture_time (float): time (seconds) when rupture starts.

  • rise_time (float): duration of rupture (seconds).

  • rise_time_starting (float or None): If None, it is internally set to rise_time/2.

  • rise_shape (str): If rise_shape == “quadratic”, the rise time function is piecewise quadratic, continuously differentiable, with maximum slope at t = rupture_time + rise_time_starting and zero slope at t = rupture_time and t = rupture_time + rise_time. If rise_shape == “linear”, the rise time function is piecewise linear, with value 0.5 at t = rupture_time + rise_time_starting.

Outputs:

rf (float or np.array): The rise time function evaluated at t. If the input is a list or tuple of times, returns a numpy array.

clawpack.geoclaw.dtopotools.strike_direction(x1, y1, x2, y2)

Calculate strike direction between two points. Actually calculates “initial bearing” from (x1,y1) in direction towards (x2,y2), following http://www.movable-type.co.uk/scripts/latlong.html