dtopotools module for moving topography

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

The notebook geoclaw/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

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)

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 ” “.

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 = None

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

property corners

Coordinates of the corners of the fault plane.

depth = None

Depth of subfault based on coordinate_specification in meters.

dip = None

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 = None

Latitutde of the subfault based on coordinate_specification.

length = None

Length of subfault in meters.

longitude = None

Longitude of the subfault based on coordinate_specification.

mu = None

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 = None

Rake of subfault movement in degrees.

rise_shape = None

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

rise_time = None

Total rise time of subfault (sec)

rise_time_starting = None

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

rupture_time = None

Time subfault starts to rupture (sec)

rupture_type = None

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 = None

Slip on subfault in strike direction in meters.

strike = None

Strike direction of subfault in degrees.

width = None

Width of subfault in meters.

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