Stratigraphy operations

The package makes available functions relating to computing things based on stratigraphy.

The functions are defined in deltametrics.strat.

Stratigraphy statistics functions

Metrics and statistics of stratigraphy.

compute_compensation(line1, line2)

Compute compensation statistic betwen two lines.

compute_net_to_gross(sediment_volume[, ...])

Compute net-to-gross for stratigraphy.

compute_thickness_surfaces(top_surface, ...)

Compute deposit thickness.

compute_sedimentograph(sediment_volume[, ...])

Compute the sedimentograph.

Compute stratigraphy routines

Routines to compute stratigraphic volumes.

compute_boxy_stratigraphy_volume(elev, prop)

Process t-x-y data volume to boxy stratigraphy volume.


Process t-x-y data volume to boxy stratigraphy coordinates.

Quick-stratigraphy attributes classes

The “quick” stratigraphy attributes provide a common API that is accessed by DataCube, DataSectionVariable and DataPlanformVariable methods. There are two methods of computing quick stratigraphy, with MeshStratigraphyAttributes as the default.

MeshStratigraphyAttributes(elev, **kwargs)

Attribute set for mesh stratigraphy information, emebdded into a DataCube.


Attribute set for boxy stratigraphy information, emebdded into a DataCube.


Low-level stratigraphy utility functions

The functions outlined in this section are the main functions that do the actual work of computing stratigraphy and preservation, throughout DeltaMetrics. These functions may be useful if you are trying to use parts of DeltaMetrics, but need to customize something for your own use-case.


Compute the preserved elevations of stratigraphy.

Given elevation data alone, we can compute the preserved stratal surfaces. These surfaces depend on the timeseries of bed elevation at any spatial location. We determine preservation by marching backward in time, determining when was the most recent time that the bed elevation was equal to a given elevation.

This function is declared as private and not part of the public API, however some users may find it helpful. The function is heavily utlized internally. Function inputs and outputs are standard numpy ndarray, so that these functions can accept data from an arbitrary source.


elev (ndarray or xr.core.dataarray.DataArray) – The t-x-y volume of elevation data to determine stratigraphy.


  • strata (ndarray) – A t-x-y ndarry of stratal surface elevations.

  • psvd (ndarray) – A t-x-y boolean ndarry of whether a (t,x,y) point of instantaneous time is preserved in any of the final stratal surfaces. To determine whether time from a given timestep is preserved, use psvd.nonzero()[0] - 1.


Compute the preserved timesteps.

The output from _compute_elevation_to_preservation records whether an instance of time, defined exactly at the data interval, is recorded in the stratigraphy (here, “recorded” does not include stasis). This differs from determining which time-intervals are preserved in the stratigraphy, because the deposits reflect the conditions between the save intervals.

While this computation is simply an offset-by-one indexing (psvd[1:, ...]), the function is implemented explicitly and utilized internally for consistency.


True in the preserved time-interval array does not necessarily indicate that an entire timestep was preserved, but rather that some portion of this time-interval (up to the entire interval) is recorded.


psvd (ndarray) – Boolean ndarray indicating the preservation of instances of time. Time is expected to be the 0th axis.


psvd_intervals – Boolean ndarray indicating the preservation of time-intervals, including partial intervals.

Return type:


deltametrics.strat._compute_preservation_to_cube(strata, z)

Compute the cube-data coordinates to strata coordinates.

Given elevation preservation data (e.g., data from _compute_elevation_to_preservation), compute the coordinate mapping from t-x-y data to z-x-y preserved stratigraphy.

While stratigraphy is time-dependent, preservation at any spatial x-y location is independent of any other location. Thus, the computation is vectorized to sweep through all “stratigraphic columns” simultaneously. The operation works by beginning at the highest elevation of the stratigraphic volume, and sweeping down though all elevations with an x-y “plate”. The plate indicates whether sediments are preserved below the current iteration elevation, at each x-y location.

Once the iteration elevation is less than the strata surface at any x-y location, there will always be sediments preserved below it, at every elevation. We simply need to determine which time interval these sediments record. Then we store this time indicator into the sparse array.

So, in the end, coordinates in resultant boxy stratigraphy are linked to t-x-y coordinates in the data source, by building a mapping that can be utilized repeatedly from a single stratigraphy computation.

This function is declared as private and not part of the public API, however some users may find it helpful. The function is heavily utlized internally. Function inputs and outputs are standard numpy ndarray, so that these functions can accept data from an arbitrary source.

  • strata (ndarray) – A t-x-y ndarry of stratal surface elevations. Can be computed by _compute_elevation_to_preservation.

  • z – Vertical coordinates of stratigraphy. Note that z does not need to have regular intervals.


  • strat_coords (ndarray) – An N x 3 array of z-x-y coordinates where information is preserved in the boxy stratigraphy. Rows in strat_coords correspond with rows in data_coords.

  • data_coords (ndarray) – An N x 3 array of t-x-y coordinates where information is to be extracted from the data array. Rows in data_coords correspond with rows in strat_coords.

deltametrics.strat._determine_strat_coordinates(elev, z=None, dz=None, nz=None)

Return a valid Z array for stratigraphy based on inputs.

This helper function enables support for user specified dz, nz, or z in many functions. The logic for determining how to handle these mutually exclusive inputs is placed in this function, to ensure consistent behavior across all classes/methods internally.


If none of the optional parameters is supplied, then a default value is used of dz=0.1.


Precedence when multiple arguments are supplied is z, dz, nz.

  • elev (ndarray) – An up-to-three-dimensional array with timeseries of elevation values, where elevation is expected to along the zeroth axis.

  • z (ndarray, optional) – Array of Z values to use for bounding intervals (i.e., points in z), returned unchanged if supplied.

  • dz (float, optional) – Interval in created Z array. Z array is created as np.arange(np.min(elev), np.max(elev)+dz, step=dz).

  • nz (int, optional) – Number of intervals in z, that is, the number of points in z is nz+1. Z array is created as np.linspace(np.min(elev), np.max (elev), num=nz, endpoint=True).

deltametrics.strat._adjust_elevation_by_subsidence(elev, sigma_dist)

Adjust elevation array by subsidence rates. Given the elevation information and information about the distance subsided, the true sedimentation and erosion pattern can be unraveled. The function is written flexibly to handle subsidence distance information provided as either a single constant value, a constant value over a 2-D spatial pattern, a 1-D temporal vector when the distance changes but is applied to the entire domain, and finally a full t-x-y 3-D array with subsidence information for each x-y position in the domain at each instance in time. Critically, when temporal information is provided (either 1-D vector or full 3-D array) values are assumed to be cumulative subsidence distances at each index in time. Alternatively, when constant information is given (1-D float/int, or 2-D x-y array) it is assumed that the constant values reported represent the distance subsided per temporal index. This function is declared as a private function and is not part of the public API. Higher level functions make a call internally to adjust the elevation data based on subsidence rate information when the optional argument sigma_dist is provided.

  • elev (ndarray or xr.core.dataarray.DataArray) – The t-x-y volume of elevation data to determine stratigraphy.

  • sigma_dist (ndarray, xr.core.dataarray.DataArray,) – float, int The subsidence information. If a single float or int is provided, that value is assumed to apply across the entire domain for all times. If a 2-D array is provided, the distances are assumed for each x-y location and are applied across all times. If a 1-D vector is provided, it is assumed that each value applies to all locations in the domain for a given time. If a full 3-D array is provided, it is assumed to follow the same t-x-y convention as the elevation data. Positive distances indicate subsidence, negative distances are uplift.


elev_adjusted – A t-x-y ndarray of the same shape as the input elev array, but the values have all been adjusted by the subsidence information provided by the input sigma_dist.

Return type:


deltametrics.strat._determine_deposit_from_background(sediment_volume, background)

Helper for stratigraphic functions.

Determine the boolean mask of the deposit, i.e., excluding the pre-deposition basin topography.

  • sediment_volume (xarray or ndarray) – The deposit voxel array. First dimension is vertical spatial dimension. In this function, the data here is used only when background is specified as a constant, or for shape if no input is given.

  • background (xarray, ndarray, or float, optional) – Value indicating the background or basin voxels that should be excluded from computation. If an array matching the size of sediment_volume is supplied, this array is assumed to be a binary array indicating the voxels to be excluded (i.e., the deposit is the inverse of background). If a float is passed, this is treated as a “no-data” value, and used to determine the background voxels. If no value is supplied, no voxels are excluded.


deposit – Boolean, with same shape as background, indicating where the volume is deposit (True).

Return type:

ndarray or DataArray


This example shows how choice of background affect the background used for stratigraphic computations.

The first background calculation simply uses the initial basin topography to determine the background voxels, but this ignores erosion and subsequent infilling.

The second background calculation considers this infilling of erosional channels by including the inflill in the deposit, and not in the background.

The third background shows the effect of using a constant value as the background value; hint: this is generally a bad idea unless your data has a background value encoded into it (like -1 or 9999).

golfcube =
golfstrat = dm.cube.StratigraphyCube.from_DataCube(golfcube, dz=0.05)

# background determined from initial basin topography
background0 = dm.strat._determine_deposit_from_background(
    background=golfstrat.Z < golfcube['eta'][0].data)
# background determined from min of bed elevation timeseries
background1 = dm.strat._determine_deposit_from_background(
    background=(golfstrat.Z < np.min(golfcube['eta'].data, axis=0)))
# background determined from a fixed sandfrac value
background2 = dm.strat._determine_deposit_from_background(

fig, ax = plt.subplots(2, 2, figsize=(6, 4))
ax[0, 0].imshow(background0[59], cmap='Greys_r')  # just below initial basin depth
ax[0, 1].imshow(background0[60], cmap='Greys_r')  # just above initial basin depth
ax[1, 0].imshow(background1[59], cmap='Greys_r')  # just below initial basin depth
ax[1, 1].imshow(background2[59], cmap='Greys_r')  # just below initial basin depth

(png, hires.png)