cpt_layering

cpt_layering(qc1Ncs, Ic, depth, **kwargs)

A function that uses the agglomerative clustering algorithm by Hudson et al. (2023) to automatically identify spatially contiguous soil layers with similar q_c1Ncs and I_c values. The algorithm automatically selects the number of layers, and returns the depths at the top and bottom of each layer, along with the corresponding q_c1Ncs and I_c values. Details are contained in the paper cited below.

Arguments

argument	format	description
qc1Ncs	Numpy array, dtype=float	overburden- and fines-corrected cone tip resistance
Ic	Numpy array, dtype=float	soil behavior type index
depth	Numpy array, dtype=float	depth to cpt measurement point

Keyword Arguments (kwargs)

keyword argument	format	description
Nmin=Nmin	integer	minimum number of layers. optional. default = 1.
Nmax=Nmax	integer	maximum number of layers. optional. default = None
tref=tref	float	constant used to define layer thickness parameter in cost function. optional. default = 0.5
wD=wD	float	weight factor for distortion score portion of cost function
wT=wT	float	weight factor for layer thickness portion of cost function

Note

1. Arguments are required and must be specified in order.
2. Keyword arguments are optional, may be specified in any order, and should be specified like "kwarg=value" as shown in the example.
3. If Nmax = Nmin, the algorithm will use a single number of layers. If Nmax > Nmin, the optimal number of layers will be obtained by incrementing the number of layers, starting at Nmin, until the minimum value of the cost function is reached, or the number of layers is equal to Nmax, whichever occurs first. If Nmax = None, it is set equal to the number of data points.

Returns

variable	format	description
ztop	Numpy array, dtype=float	depth to top of each layer
zbot	Numpy array, dtype=float	depth to bottom of each layer
qc1Ncs_lay	Numpy array, dtype=float	qc1Ncs values for each layer
Ic_lay	float	Ic values for each layer

Note

The ztop and zbot values returned by the algorithm result in a spatially contiguous profile, without gaps, such that the sum of layer thicknesses, defined as zbot - ztop, is the same as the profile thickness. Layer boundaries are set at the midpoint between adjacent clusters, as illustrated by the example below.

Examples

Simplified 3-layer system to illustrate basic features of behavior

Jupyter notebook cpt_layering_simple.ipynb

import numpy as np
import ngl_tools.smt as smt
import matplotlib.pyplot as plt

qc1Ncs = np.asarray([10, 10, 10, 20, 20, 20, 30, 30, 30], dtype=float)
Ic = np.asarray([3, 3, 3, 2, 2, 2, 1, 1, 1], dtype=float)
depth = np.asarray([0, 1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
ztop, zbot, qc1Ncs_lay, Ic_lay = smt.cpt_layering(qc1Ncs, Ic, depth, dGWT=1.0, Nmin=1, Nmax=4)

def interleave(v1, v2):
    vout = np.empty(2 * len(v1), dtype=float)
    vout[0::2] = v1
    vout[1::2] = v2
    return vout 

fig, ax = plt.subplots(ncols=2, sharey='row')
ax[0].scatter(qc1Ncs, depth, facecolor='white', edgecolor='C0', s=10, zorder=2, label='data')
ax[1].scatter(Ic, depth, facecolor='white', edgecolor='C0', s=10, zorder=2)
ax[0].plot(interleave(qc1Ncs_lay, qc1Ncs_lay), interleave(ztop, zbot), c='C1', zorder=1, label='clustered')
ax[1].plot(interleave(Ic_lay, Ic_lay), interleave(ztop, zbot), c='C1', zorder=1)

for a in ax:
    a.grid(True, alpha=0.5)

ax[0].set_ylabel('depth')
ax[0].set_xlabel(r'$q_{c1Ncs}$')
ax[1].set_xlabel(r'$I_c$')
ax[0].invert_yaxis()
ax[0].legend()

Applied to real CPT data

To download the files below, right click the link and use "Save Link As".

Data file example_cpt.csv

Jupyter notebook cpt_layering_profile.ipynb

import numpy as np
import json
import pandas as pd
import requests
import io
import matplotlib.pyplot as plt
import smt

def interleave(v1, v2):
    vout = np.empty(2 * len(v1), dtype=float)
    vout[0::2] = v1
    vout[1::2] = v2
    return vout 

# Read example CPT file
df = pd.read_csv('example_cpt.csv')
depth = df['depth'].values
qc = df['qc'].values
fs = df['fs'].values

# Assign groundwater table depth, average unit weight, and compute vertical total and effective stresses
dGWT = 3.0
gamma = 19.0
sigmav = gamma * depth
u = (depth - dGWT) * 9.81
u[u<0] = 0
sigmavp = sigmav - u

# Use inverse filtering algorithm and layering algorithm to compute layered profile
qt_inv, fs_inv, Ic_inv = smt.cpt_inverse_filter(qc, depth, fs=fs, sigmav=sigmav, sigmavp=sigmavp)
Ic, Qtn, Fr = smt.get_Ic_Qtn_Fr(qc, fs, sigmav, sigmavp)
Ic_inv, Qtn_inv, Fr_inv = smt.get_Ic_Qtn_Fr(qt_inv, fs_inv, sigmav, sigmavp)
FC = smt.get_FC_from_Ic(Ic_inv, 0.0)
qc1N, qc1Ncs = smt.get_qc1N_qc1Ncs(qc, fs, sigmav, sigmavp, FC)
qc1N_inv, qc1Ncs_inv = smt.get_qc1N_qc1Ncs(qt_inv, fs_inv, sigmav, sigmavp, FC)
ztop, zbot, qc1Ncs_lay, Ic_lay = smt.cpt_layering(qc1Ncs_inv, Ic_inv, depth, dGWT=dGWT, Nmin=1, Nmax=None)
fig, ax = plt.subplots(ncols=4, sharey='row')

ax[0].plot(qc, depth)
ax[0].plot(qt_inv, depth)
ax[1].plot(fs, depth)
ax[1].plot(fs_inv, depth)
ax[2].plot(Ic, depth)
ax[2].plot(Ic_inv, depth)
ax[3].plot(qc1Ncs, depth)
ax[3].plot(qc1Ncs_inv, depth)
ax[2].axvspan(1,1.31,alpha=0.3,color='darkorange')
ax[2].axvspan(1.31,2.05,alpha=0.3,color='peru')
ax[2].axvspan(2.05,2.6,alpha=0.3,color='palegreen')
ax[2].axvspan(2.6,2.95,alpha=0.3,color='mediumaquamarine')
ax[2].axvspan(2.95,3.6,alpha=0.3,color='navy')
ax[2].axvspan(3.6,4.0,alpha=0.3,color='saddlebrown')
ax[0].set_ylabel('depth [m]')
ax[0].set_xlabel(r'$q_c$ [kPa]')
ax[1].set_xlabel(r'$f_s$ [kPa]')
ax[2].set_xlabel(r'$I_c$ [-]')
ax[3].set_xlabel(r'$q_{c1Ncs}$ [-]')
for a in ax:
    a.grid(True, alpha=0.5)
ax[2].plot(interleave(Ic_lay,Ic_lay), interleave(ztop,zbot), color='black', zorder=1)
ax[3].plot(interleave(qc1Ncs_lay,qc1Ncs_lay), interleave(ztop,zbot), color='black', zorder=1)
ax[0].invert_yaxis()
ax[2].set_xlim(1,4)