Skip to content

Package structure

All the modules and functions are documented in the following sections.

pynewmarkdisp.newmark

Functions to calculate permanent displacements using empirical correlations.

Permanent displacements are calculated based on the Newmark (1965) sliding block method.

direct_newmark 

direct_newmark(time, accel, ky, g)

Calculate the permanent displacements using the direct Newmark method.

Parameters:

Name Type Description Default
time (n,) ndarray

1D array with the time series of the earthquake record, in [s].

 required
accel (n,) ndarray

1D array with the acceleration series of the earthquake record.

 required
ky float

Critical seismic coefficient (yield coefficient).

 required
g float

Acceleration of gravity, in the same units as accel. For example, if accel is in [m/s²], g should be equal to 9.81; if accel is as a fraction of gravity, g should be equal to 1.

 required

Returns:

Name Type Description
newmark_str dict

Dictionary with the structure from the Newmark's method. The structure includes time, acceleration, velocity, and displacements series, as well as the critical seismic coefficient and accelerarion, and the permanent displacement calculated (in meters).
Source code in pynewmarkdisp/newmark.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
def direct_newmark(time, accel, ky, g):
    """
    Calculate the permanent displacements using the direct Newmark method.

    Parameters
    ----------
    time : (n,) ndarray
        1D array with the time series of the earthquake record, in [s].
    accel : (n,) ndarray
        1D array with the acceleration series of the earthquake record.
    ky : float
        Critical seismic coefficient (yield coefficient).
    g : float
        Acceleration of gravity, in the same units as ``accel``. For example,
        if ``accel`` is in [m/s²], ``g`` should be equal to 9.81; if ``accel``
        is as a fraction of gravity, ``g`` should be equal to 1.

    Returns
    -------
    newmark_str : dict
        Dictionary with the structure from the Newmark's method. The structure
        includes time, acceleration, velocity, and displacements series, as
        well as the critical seismic coefficient and accelerarion, and the
        permanent displacement calculated (in meters).
    """
    length = len(time)
    ay = 9.81 * ky  # Yield acceleration to SI units
    accel = 9.81 * accel / g  # Earthquake acceleration to SI units
    if ay == 0 or np.isnan(
        ay
    ):  # The slope is already unstable under static condition
        vel = np.full(length, np.nan)
        disp = np.full(length, np.nan)
    elif accel.max() > ay:  # Evaluate only if the earthquake exceeds `ay`
        # Velocity
        vel = first_newmark_integration(time, accel, ay)
        # Displacement
        disp = cumul_integral(y=vel, x=time, initial=0)
    else:  # The slope is stable even under seudostatic condition
        vel = np.zeros_like(accel)
        disp = np.zeros_like(accel)

    return {
        "time": time,
        "accel": accel,
        "vel": vel,
        "disp": disp,
        "ky": ky,
        "ay": ay,
        "perm_disp": disp[-1],
    }

empirical_correlations 

empirical_correlations(**kwargs)

Run a function of an empirical correlation for uₚ calculation.

It is possible to call a function from the empir_corr.py module by passing the ID of the function to the opt keyword argument, or by passing the function itself to the correlation_funct keyword argument.

The opt keyword argument follows the numbering of the functions in the Table 1 of Meehan & Vahedifard (2013) (see the empir_corr.py module for more details).

After choosing either opt or correlation_funct, the remaining keyword arguments are passed to the chosen function.

Returns:

Name Type Description
perm_disp float

Permanent displacement calculated by the chosen function using an empirical correlation.
Source code in pynewmarkdisp/newmark.py 
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def empirical_correlations(**kwargs):
    """Run a function of an empirical correlation for uₚ calculation.

    It is possible to call a function from the ``empir_corr.py`` module by
    passing the ID of the function to the ``opt`` keyword argument, or by
    passing the function itself to the ``correlation_funct`` keyword argument.

    The ``opt`` keyword argument follows the numbering of the functions in the
    Table 1 of [Meehan & Vahedifard (2013)](https://doi.org/10.1016/j.enggeo.2012.10.016)
    (see the ``empir_corr.py`` module for more details).

    After choosing either ``opt`` or ``correlation_funct``, the remaining
    keyword arguments are passed to the chosen function.

    Returns
    -------
    perm_disp : float
        Permanent displacement calculated by the chosen function using an
        empirical correlation.
    """
    if "opt" in kwargs:
        correlation_funct = correlations[kwargs["opt"]]
    else:
        correlation_funct = kwargs["correlation_funct"]
    return correlation_funct(**kwargs)

first_newmark_integration 

first_newmark_integration(time, accel, ay)

Perform the first integration of the Newmark method.

The integration is performed using the trapezoidal rule. It accounts for energy dissipation at each peak above the critical acceleration by subtracting the area between $y=a_y$ and $y=a(t)$ after the peak until the velocity is zero.

This funtion is optimized for speed using Numba.

Parameters:

Name Type Description Default
time (n, ) ndarray

1D array with the time series of the earthquake record, in units consistent with accel and ay.

 required
accel (n, ) ndarray

1D array with the acceleration series of the earthquake record, in the same units as ay.

 required
ay int or float

Critical acceleration, in the same units as accel.

 required

Returns:

Name Type Description
vel (n, ) ndarray

1D array with the velocity series obtained by the first integration of the Newmark method. Units are consistent with accel and ay.
Source code in pynewmarkdisp/newmark.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
@njit(fastmath=True, cache=True)
def first_newmark_integration(time, accel, ay):
    """Perform the first integration of the Newmark method.

    The integration is performed using the trapezoidal rule. It accounts for
    energy dissipation at each peak above the critical acceleration by
    subtracting the area between $y=a_y$ and $y=a(t)$ after the peak until the
    velocity is zero.

    This funtion is optimized for speed using Numba.

    Parameters
    ----------
    time : (n, ) ndarray
        1D array with the time series of the earthquake record, in units
        consistent with ``accel`` and ``ay``.
    accel : (n, ) ndarray
        1D array with the acceleration series of the earthquake record, in the
        same units as ``ay``.
    ay : int or float
        Critical acceleration, in the same units as ``accel``.

    Returns
    -------
    vel : (n, ) ndarray
        1D array with the velocity series obtained by the first integration of
        the Newmark method. Units are consistent with ``accel`` and ``ay``.
    """
    length = len(time)
    vel = np.zeros(length)
    for i in np.arange(1, length, 1):
        if accel[i] > ay:
            v = vel[i - 1] + trapz(
                y=accel[i - 1 : i + 1] - ay, x=time[i - 1 : i + 1]
            )
        elif accel[i] < ay and vel[i - 1] > 0:
            v = vel[i - 1] - abs(
                trapz(y=accel[i - 1 : i + 1], x=time[i - 1 : i + 1])
            )
        else:
            v = 0
        vel[i] = max(v, 0)
    return vel

plot_newmark_integration 

plot_newmark_integration(newmark_str, compressed=False)

Plot the double-integration procedure of the direct Newmark method.

Parameters:

Name Type Description Default
newmark_str dict

Dictionary with the structure from the Newmark's method. The structure includes time, acceleration, velocity, and displacements series, as well as the critical seismic coefficient and accelerarion, and the permanent displacement calculated (in meters). This dictionary is returned by the direct_newmark function.

 required
compressed bool

If True, the plot will be compressed in a single row. Otherwise, the plot will be in three rows. Default is False.

 False

Returns:

Name Type Description
fig matplotlib.figure.Figure

Matplotlib object which might be used to save the figure as a file.
Source code in pynewmarkdisp/newmark.py 
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
def plot_newmark_integration(newmark_str, compressed=False):
    """
    Plot the double-integration procedure of the direct Newmark method.

    Parameters
    ----------
    newmark_str : dict
        Dictionary with the structure from the Newmark's method. The structure
        includes time, acceleration, velocity, and displacements series, as
        well as the critical seismic coefficient and accelerarion, and the
        permanent displacement calculated (in meters). This dictionary is
        returned by the ``direct_newmark`` function.
    compressed : bool, optional
        If ``True``, the plot will be compressed in a single row. Otherwise,
        the plot will be in three rows. Default is ``False``.

    Returns
    -------
    fig : matplotlib.figure.Figure
        Matplotlib object which might be used to save the figure as a file.
    """
    if compressed:
        fig, ax0 = plt.subplots(ncols=1, nrows=1, figsize=[6.5, 3])
        ax1 = ax0.twinx()
        ax2 = ax0.twinx()
        colors = ["#997700", "#BB5566", "#004488"]
    else:
        fig, axs = plt.subplots(
            ncols=1, nrows=3, figsize=[6.5, 5], sharex=True
        )
        ax0, ax1, ax2 = axs
        colors = ["#004488", "#004488", "#004488"]
    # Acceleration
    ax0.plot(newmark_str["time"], newmark_str["accel"], lw=0.75, c=colors[0])
    ay = newmark_str["ay"]
    ky = ay / 9.81
    l1 = ax0.axhline(
        ay,
        c="k",
        ls="--",
        label="$a_\\mathrm{y}=$" + f"{ky:.2f}g",  # = " + f"{ay:.2f} m/s$^2$",
    )
    ax0.set(ylabel="Acceleration, $a$ [m s$^{- 2}$]")
    # Velocity
    ax1.plot(newmark_str["time"], newmark_str["vel"], lw=0.75, color=colors[1])
    ax1.set(ylabel="Velocity, $v$ [m s$^{- 1}$]")
    # Displacement
    lw = 1.2 if compressed else 1.0
    ax2.plot(
        newmark_str["time"],
        newmark_str["disp"],
        lw=lw,
        c=colors[2],
    )
    (p1,) = ax2.plot(
        newmark_str["time"][-1],
        newmark_str["disp"][-1],
        ls="",
        marker="*",
        mfc="white",
        c=colors[2],
        label="$u_\\mathrm{p}=$" + f"{newmark_str['disp'][-1]:.3f} m",
    )
    ax2.set(ylabel="Displacement, $u$ [m]", xlabel="Time, $t$ [s]")
    # Setup
    for ax, color in zip([ax0, ax1, ax2], colors):
        color_lbl = color if compressed else "k"
        ax.yaxis.set_major_formatter(FormatStrFormatter("%.2f"))
        ax.spines["bottom"].set_linewidth(1.5)
        ax.spines["left"].set_linewidth(1.5)
        ax.grid(False)
        ax.yaxis.label.set_color(color_lbl)
        ax.tick_params(axis="y", colors=color_lbl)
    if compressed:  # Adjust the position of the axis and its label
        ax1.spines["left"].set_visible(True)  # Make the spine visible
        ax1.spines["left"].set_position(("outward", 50))
        ax1.yaxis.set_label_position("left")
        ax1.yaxis.set_ticks_position("left")  # Position the ticks on the left
        ax0.grid(True, which="major", linestyle="--")
        fig.legend(
            handles=[l1, p1],
            loc="lower right",
            bbox_to_anchor=(0.9, 0.2),
        )
        ax0.set(xlabel="Time, $t$ [s]")
    else:
        [ax.grid(True, which="major", linestyle="--") for ax in axs]
        ax0.legend(loc="best")
        ax2.legend(loc="lower right")
    fig.tight_layout(pad=0.1)
    return fig

trapz 

trapz(y, x)

Calculate the area under de curve y=f(x) by the trapezoidal rule method.

This funtion is optimized for speed using Numba.

Parameters:

Name Type Description Default
y (n, ) ndarray

1D array with the values of the function y=f(x).

 required
x (n, ) ndarray

1D array with the values of the independent variable x.

 required

Returns:

Name Type Description
area float

Area under the curve y=f(x) calculated by the trapezoidal rule method.
Source code in pynewmarkdisp/newmark.py 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
@njit(fastmath=True, cache=True)
def trapz(y, x):
    """Calculate the area under de curve y=f(x) by the trapezoidal rule method.

    This funtion is optimized for speed using Numba.

    Parameters
    ----------
    y : (n, ) ndarray
        1D array with the values of the function y=f(x).
    x : (n, ) ndarray
        1D array with the values of the independent variable x.

    Returns
    -------
    area : float
        Area under the curve y=f(x) calculated by the trapezoidal rule method.
    """
    return 0.5 * ((x[1:] - x[:-1]) * (y[1:] + y[:-1])).sum()

pynewmarkdisp.spatial

Functions to calculate permanent displacements in a spatial domain.

get_dip_azimuth 

get_dip_azimuth(dem, cellsize=1)

Calculate the azimuth of the slope dip direction from the North.

Positive values are clockwise from the North. The azimuth is returned in [°], ranging from 0 to 360.

Parameters:

Name Type Description Default
dem (m, n) array

2D array with the digital elevation model (DEM) of the area.

 required
cellsize int

Cell size of the raster representing the DEM. Its default value is 1.

 1

Returns:

Name Type Description
azimuth (m, n) array

Spatial distribution of the azimuth of the slope dip direction.
Source code in pynewmarkdisp/spatial.py 
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
def get_dip_azimuth(dem, cellsize=1):
    """Calculate the azimuth of the slope dip direction from the North.

    Positive values are clockwise from the North. The azimuth is returned in
    [°], ranging from 0 to 360.

    Parameters
    ----------
    dem : (m, n) array
        2D array with the digital elevation model (DEM) of the area.
    cellsize : int, optional
        Cell size of the raster representing the DEM. Its default value is 1.

    Returns
    -------
    azimuth : (m, n) array
        Spatial distribution of the azimuth of the slope dip direction.
    """
    dzdy, dzdx = np.gradient(dem, cellsize)
    dzdx *= -1  # Adequate sign for the x component
    return np.degrees(0.5 * np.pi - np.arctan2(dzdy, dzdx)) % 360

get_hillshade 

get_hillshade(dem, sun_azimuth=315, sun_altitude=30, cellsize=1)

Generate a hillshade image from the DEM.

Parameters:

Name Type Description Default
dem (m, n) array

2D array with the digital elevation model (DEM) of the area

 required
sun_azimuth int or float

Azimuth of the sun, in [°]. Its default value is 315.

 315
sun_altitude int or float

Angle from the horizontal plane to the sun, in degrees, indicating the altitude of the sun. Its default value is 30.

 30
cellsize int

Cell size of the raster representing the DEM. Units must be consistent with the units of dem. Its default value is 1.

 1

Returns:

Name Type Description
hillshade (m, n) array

2D array with the hillshade raster.
Source code in pynewmarkdisp/spatial.py 
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
def get_hillshade(dem, sun_azimuth=315, sun_altitude=30, cellsize=1):
    """Generate a hillshade image from the DEM.

    Parameters
    ----------
    dem : (m, n) array
        2D array with the digital elevation model (DEM) of the area
    sun_azimuth : int or float, optional
        Azimuth of the sun, in [°]. Its default value is 315.
    sun_altitude : int or float, optional
        Angle from the horizontal plane to the sun, in degrees, indicating the
        altitude of the sun. Its default value is 30.
    cellsize : int, optional
        Cell size of the raster representing the DEM. Units must be consistent
        with the units of ``dem``. Its default value is 1.

    Returns
    -------
    hillshade : (m, n) array
        2D array with the hillshade raster.
    """
    sun_azimuth_r = np.radians(sun_azimuth)
    sun_altitude_r = np.radians(sun_altitude)
    slope_r = np.radians(get_slope(dem, cellsize))
    azimuth_r = np.radians(get_dip_azimuth(dem, cellsize))
    # Lambertian reflectance
    term1 = (
        np.cos(azimuth_r - sun_azimuth_r)
        * np.sin(slope_r)
        * np.sin(0.5 * np.pi - sun_altitude_r)
    )
    term2 = np.cos(slope_r) * np.cos(0.5 * np.pi - sun_altitude_r)
    reflectance = term1 + term2
    ptp = np.nanmax(reflectance) - np.nanmin(reflectance)
    return 255 * (reflectance - np.nanmin(reflectance)) / ptp

get_idx_at_coords 

get_idx_at_coords(x, y, spat_ref)

Convert coordinates to array indexes.

Parameters:

Name Type Description Default
x int or float

Coordinate in the x-direction.

 required
y int or float

Coordinate in the y-direction.

 required
spat_ref dict

Dictionary with the spatial reference of the raster. It contains the same information of the ASCII file header. See the documentation of the read_ascii function for more information.

 required

Returns:

Name Type Description
indexes tuple

Tuple with the indexes of the array corresponding to the given coordinates.
Source code in pynewmarkdisp/spatial.py 
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
def get_idx_at_coords(x, y, spat_ref):
    """Convert coordinates to array indexes.

    Parameters
    ----------
    x : int or float
        Coordinate in the x-direction.
    y : int or float
        Coordinate in the y-direction.
    spat_ref : dict
        Dictionary with the spatial reference of the raster. It contains the
        same information of the ASCII file header. See the documentation of
        the ``read_ascii`` function for more information.

    Returns
    -------
    indexes : tuple
        Tuple with the indexes of the array corresponding to the given
        coordinates.
    """
    i = int(
        spat_ref["nrows"]
        - np.ceil((y - spat_ref["yllcorner"]) / spat_ref["cellsize"])
    )
    j = int(np.floor((x - spat_ref["xllcorner"]) / spat_ref["cellsize"]))
    return i, j

get_slope 

get_slope(dem, cellsize=1)

Calculate the slope (or inclination) of the terrain.

Parameters:

Name Type Description Default
dem (m, n) array

2D array with the digital elevation model (DEM) of the area.

 required
cellsize int

Cell size of the raster representing the DEM. Units must be consistent with the units of dem. Its default value is 1.

 1

Returns:

Name Type Description
slope (m, n) array

Spatial distribution of the slope, in degrees.
Source code in pynewmarkdisp/spatial.py 
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
def get_slope(dem, cellsize=1):
    """Calculate the slope (or inclination) of the terrain.

    Parameters
    ----------
    dem : (m, n) array
        2D array with the digital elevation model (DEM) of the area.
    cellsize : int, optional
        Cell size of the raster representing the DEM. Units must be consistent
        with the units of ``dem``. Its default value is 1.

    Returns
    -------
    slope : (m, n) array
        Spatial distribution of the slope, in degrees.
    """
    dzdy, dzdx = np.gradient(dem, cellsize)
    dzdx *= -1  # Adequate sign for the x component
    return np.degrees(np.arctan(np.sqrt(dzdx**2 + dzdy**2)))

hzt_accel_from_2_dir 

hzt_accel_from_2_dir(accel_ew, accel_ns, azimuth)

Calculate the resultant acceleration from two horizontal components.

This funtion is optimized for speed using Numba.

Parameters:

Name Type Description Default
accel_ew (n,) ndarray

1D array with the acceleration series in the East-West direction.

 required
accel_ns (n,) ndarray

1D array with the acceleration series in the North-South direction.

 required
azimuth float or int

Azimuth of the slope dip direction, in [°] and measured clockwise from North.

 required

Returns:

Name Type Description
accel_hzt (n, ) ndarray

1D array with the resultant horizontal acceleration series.
Source code in pynewmarkdisp/spatial.py 
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
@njit(fastmath=True, cache=True)
def hzt_accel_from_2_dir(accel_ew, accel_ns, azimuth):
    """Calculate the resultant acceleration from two horizontal components.

    This funtion is optimized for speed using Numba.

    Parameters
    ----------
    accel_ew : (n,) ndarray
        1D array with the acceleration series in the East-West direction.
    accel_ns : (n,) ndarray
        1D array with the acceleration series in the North-South direction.
    azimuth : float or int
        Azimuth of the slope dip direction, in [°] and measured clockwise from
        North.

    Returns
    -------
    accel_hzt : (n, ) ndarray
        1D array with the resultant horizontal acceleration series.
    """
    azimuth_r = np.radians(azimuth)
    return accel_ew * np.sin(azimuth_r) + accel_ns * np.cos(azimuth_r)

load_ascii_raster 

load_ascii_raster(path, dtype=None)

Load a raster file in ASCII format (.ASC extension).

Parameters:

Name Type Description Default
path str

File to be loaded. It must be in Esri ASCII raster format. Its header must be 6 lines long, and the spatial location of the raster is specified by the location of the lower left corner of the lower left cell, as shown below:

NCOLS xxx
NROWS xxx
XLLCORNER xxx
YLLCORNER xxx
CELLSIZE xxx
NODATA_VALUE xxx
row 1
row 2
...
row n
required
dtype data-type

Data type of the returned array. If not given (i.e., dtype=None), the data type is determined by the contents of the file. Its default value is None.

 None

Returns:

Name Type Description
raster (m, n) ndarray

Raster data given as a 2D array.
header dict

Header of the raster file.
Source code in pynewmarkdisp/spatial.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def load_ascii_raster(path, dtype=None):
    """Load a raster file in ASCII format (`.ASC` extension).

    Parameters
    ----------
    path : str
        File to be loaded. It must be in [Esri ASCII raster format](
        https://en.wikipedia.org/wiki/Esri_grid).
        Its header must be 6 lines long, and the spatial location of the raster
        is specified by the location of the lower left corner of the lower
        left cell, as shown below:
        ```
            NCOLS xxx
            NROWS xxx
            XLLCORNER xxx
            YLLCORNER xxx
            CELLSIZE xxx
            NODATA_VALUE xxx
            row 1
            row 2
            ...
            row n
        ```
    dtype : data-type, optional
        Data type of the returned array. If not given (i.e., `dtype=None`), the
        data type is determined by the contents of the file. Its default value
        is `None`.

    Returns
    -------
    raster : (m, n) ndarray
        Raster data given as a 2D array.
    header : dict
        Header of the raster file.
    """
    raster = np.loadtxt(path, skiprows=6)
    header = np.loadtxt(path, max_rows=6, dtype=object)
    header = {
        "ncols": int(header[0, 1]),
        "nrows": int(header[1, 1]),
        "xllcorner": float(header[2, 1]),
        "yllcorner": float(header[3, 1]),
        "cellsize": float(header[4, 1]),
        "nodata_value": int(header[5, 1]),
    }
    raster[np.where(raster == header["nodata_value"])] = np.nan
    raster = raster.astype(dtype) if dtype is not None else raster
    return (raster, header)

map_zones 

map_zones(parameters, zones)

Associate geotechnical parameters to zones.

Parameters:

Name Type Description Default
parameters dict

Dictionary with the geotechnical parameters of each zone. The keys are the zone numbers, and the values are tuples with the zone's friction angle in [°], cohesion in [kPa], and unit weight in [kN/m³], respectively.

 required
zones (m, n) ndarray

2D array with the zone number assigned to each cell. There has to be as many zone numbers as keys in the parameters dictionary.

 required

Returns:

Type Description
phi, c, gamma : (m, n) ndarray, (m, n) ndarray, (m, n) ndarray

Three 2D arrays with the spatial distribution of the friction angle, cohesion, and unit weight, respectively. Each array has the same shape as zones.
Source code in pynewmarkdisp/spatial.py 
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def map_zones(parameters, zones):
    """Associate geotechnical parameters to zones.

    Parameters
    ----------
    parameters : dict
        Dictionary with the geotechnical parameters of each zone. The keys are
        the zone numbers, and the values are tuples with the zone's friction
        angle in [°], cohesion in [kPa], and unit weight in [kN/m³],
        respectively.
    zones : (m, n) ndarray
        2D array with the zone number assigned to each cell. There has to be as
        many zone numbers as keys in the parameters dictionary.

    Returns
    -------
    phi, c, gamma : (m, n) ndarray, (m, n) ndarray, (m, n) ndarray
        Three 2D arrays with the spatial distribution of the friction angle,
        cohesion, and unit weight, respectively. Each array has the same shape
        as ``zones``.
    """
    phi = np.empty_like(zones)
    c = np.empty_like(zones)
    gamma = np.empty_like(zones)
    for z in parameters.keys():  # Zone: frict_angle, cohesion, gamma
        phi[np.where(zones == z)] = parameters[z][0]
        c[np.where(zones == z)] = parameters[z][1]
        gamma[np.where(zones == z)] = parameters[z][2]
    return phi, c, gamma

plot_spatial_field 

plot_spatial_field(field, dem=None, discrete=False, **kwargs)

Plot a spatial field.

Parameters:

Name Type Description Default
field (m, n) ndarray

2D array with the spatial field to be plotted.

 required
dem (m, n) ndarray or None

2D array with the digital elevation model (DEM) of the area. It is used to plot the hillshade below the spatial field to enhance the visualization. By default its value is None.

 None
discrete bool

If True, the spatial field is plotted as a discrete field. By default its value is False.

 False
**kwargs dict

Additional keyword arguments. spat_ref is a dictionary with the spatial reference of the raster. It contains the same information of the ASCII file header. See the documentation of the read_ascii function for more information. levels is a list or array with the contour levels. cmap is a colormap palette accepted by Matplotlib, or a str with the name of a Matplotlib colormap. label is a list of str with the labels of the unique elements in a discrete field. vmin and vmax define the data range that the colormap covers. labelrot is the rotation angle of the colorbar ticks labels. title is the title of the spatial field.

 {}

Returns:

Name Type Description
fig matplotlib.figure.Figure

Matplotlib object which might be used to save the figure as a file.
Source code in pynewmarkdisp/spatial.py 
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
def plot_spatial_field(field, dem=None, discrete=False, **kwargs):
    """
    Plot a spatial field.

    Parameters
    ----------
    field : (m, n) ndarray
        2D array with the spatial field to be plotted.
    dem : (m, n) ndarray or None, optional
        2D array with the digital elevation model (DEM) of the area. It is used
        to plot the hillshade below the spatial field to enhance the
        visualization. By default its value is ``None``.
    discrete : bool, optional
        If ``True``, the spatial field is plotted as a discrete field. By
        default its value is ``False``.
    **kwargs : dict, optional
        Additional keyword arguments. ``spat_ref`` is a dictionary with the
        spatial reference of the raster. It contains the same information of
        the ASCII file header. See the documentation of the ``read_ascii``
        function for more information. ``levels`` is a list or array with the
        contour levels. ``cmap`` is a colormap palette accepted by
        Matplotlib, or a ``str`` with the name of a Matplotlib colormap.
        ``label`` is a list of ``str`` with the labels of the unique elements
        in a discrete field. ``vmin`` and ``vmax`` define the data range that
        the colormap covers. ``labelrot`` is the rotation angle of the colorbar
        ticks labels. ``title`` is the title of the spatial field.

    Returns
    -------
    fig : matplotlib.figure.Figure
        Matplotlib object which might be used to save the figure as a file.
    """
    m, n = field.shape
    alpha = 1
    spat_ref = kwargs.get(
        "spat_ref", {"xllcorner": 0, "yllcorner": 0, "cellsize": 1}
    )
    extent = (
        spat_ref["xllcorner"],  # x min
        spat_ref["xllcorner"] + spat_ref["cellsize"] * n,  # x max
        spat_ref["yllcorner"],  # y min
        spat_ref["yllcorner"] + spat_ref["cellsize"] * m,  # y max
    )

    fig, ax = plt.subplots(ncols=1, nrows=1, figsize=[5.5, 4.5], sharex=True)
    if dem is not None:  # Generate hillshade in background and contours
        alpha = 0.8  # Transparency for ploting the field over the hillshade
        hsd = get_hillshade(dem, cellsize=spat_ref["cellsize"])
        ax.matshow(hsd, cmap="gray", extent=extent)
        # ax.matshow(dem, cmap="gray", extent=extent, alpha=alpha)
        levels = kwargs.get("levels")
        cs = ax.contour(
            dem,
            colors="k",
            origin="upper",
            linewidths=0.5,
            extent=extent,
            levels=levels,
            alpha=0.8,
        )
        ax.clabel(cs, inline=True, fontsize="x-small")

    if discrete:  # Plot field as a discrete variable with unique values
        ticks = np.unique(field)[~np.isnan(np.unique(field))]  # For colorbar
        n = len(ticks)
        field_r = np.full_like(field, np.nan)  # Empty reclasified field
        ticks_r = np.arange(n)  # Ticks to use in reclasified field
        for i in ticks_r:  #  Fill reclasified field
            field_r[field == ticks[i]] = i
        cmap = plt.colormaps.get_cmap(kwargs.get("cmap")).resampled(n)  # type: ignore
        im = ax.matshow(
            field_r,  # plots the reclasified field instead of the original
            cmap=cmap,
            extent=extent,
            alpha=alpha,
            vmin=-0.5,
            vmax=n - 0.5,
        )  # map reclasified, but use the original labels
        cbar = fig.colorbar(im, ax=ax, shrink=0.75, pad=0.01, ticks=ticks_r)
        ticks = np.round(ticks, 1) if ticks.dtype == np.float64 else ticks
        label = kwargs.get("label", ticks)
        _ = cbar.ax.set_yticklabels(label)
    else:  # Plot field as a continue variable
        cmap = plt.colormaps.get_cmap(kwargs.get("cmap"))  # type: ignore
        vmin = kwargs.get("vmin", np.nanmin(field))
        vmax = kwargs.get("vmax", np.nanmax(field))

        field_isfin = field[np.isfinite(field)]
        if np.nanmin(field_isfin) < vmin and np.nanmax(field_isfin) > vmax:
            extend = "both"
        elif np.nanmin(field_isfin) < vmin:
            extend = "min"
        elif np.nanmax(field_isfin) > vmax:
            extend = "max"
        else:
            extend = "neither"
        im = ax.matshow(
            field,
            cmap=cmap,
            extent=extent,
            alpha=alpha,
            vmin=kwargs.get("vmin"),
            vmax=kwargs.get("vmax"),
        )
        cbar = fig.colorbar(im, ax=ax, shrink=0.75, pad=0.01, extend=extend)

    # Colorbar settings
    lr = kwargs.get("labelrot", 0)
    cbar.ax.tick_params(axis="y", labelrotation=lr, labelsize="small", pad=1.5)
    cbar.set_label(kwargs.get("title", "Scale"), rotation=90, size="large")
    [label.set_va("center") for label in cbar.ax.get_yticklabels()]
    # Axis settings
    ax.xaxis.set_tick_params(labelrotation=0, labelsize="small")
    ax.yaxis.set_tick_params(labelrotation=90, labelsize="small")
    ax.yaxis.set_major_formatter(FormatStrFormatter("%.0f"))
    ax.xaxis.set_major_formatter(FormatStrFormatter("%.0f"))
    [label.set_va("center") for label in ax.get_yticklabels()]
    for spine in ["bottom", "top", "left", "right"]:
        ax.spines[spine].set_linewidth(1.5)
    ax.grid(True, which="major", linestyle="--")
    fig.tight_layout()
    return fig

reclass_raster 

reclass_raster(field, limits)

Reclassify a raster based on a list of limits.

Convert a continuous raster into a categorical one. The limits are used to define the categories.

Parameters:

Name Type Description Default
field (m, n) ndarray

2D array with the spatial field to be reclassified.

 required
limits array_like

Limits to be used in the reclassification. The minimum and maximum values of the field should not be included as the first and last elements of the list, as they are automatically included in the limits.

 required

Returns:

Name Type Description
field_reclass (m, n) ndarray

2D array with the reclassified field.
Source code in pynewmarkdisp/spatial.py 
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
def reclass_raster(field, limits):
    """Reclassify a raster based on a list of limits.

    Convert a continuous raster into a categorical one. The limits are used to
    define the categories.

    Parameters
    ----------
    field : (m, n) ndarray
        2D array with the spatial field to be reclassified.
    limits : array_like
        Limits to be used in the reclassification. The minimum and maximum
        values of the field should *not* be included as the first and last
        elements of the list, as they are automatically included in the limits.

    Returns
    -------
    field_reclass : (m, n) ndarray
        2D array with the reclassified field.
    """
    limits = np.hstack((np.nanmin(field), limits, np.nanmax(field)))
    field_reclass = (len(limits) - 1) * np.ones_like(field, dtype=int)
    for i, lim in enumerate(limits[:-1]):
        field_reclass[
            np.logical_and(field >= limits[i], field <= limits[i + 1])
        ] = i
    return field_reclass

spatial_hzt_accel_from_2_dir 

spatial_hzt_accel_from_2_dir(accel_ew, accel_ns, azimuth)

Calculate spatially the resultant acceleration from two hztl components.

The azimuth is spatially distributed (e.g., obtained using the get_dip_azimuth function), but accel_ew, accel_ns, and time_vect are constant over all the spatial domain.

Parameters:

Name Type Description Default
accel_ew (p,) ndarray

1D array with the acceleration series in the East-West direction.

 required
accel_ns (p,) ndarray

1D array with the acceleration series in the North-South direction.

 required
azimuth (m, n) ndarray

2D array with the spatial distribution of the azimuth of the slope dip direction, in [°] and measured clockwise from North.

 required

Returns:

Name Type Description
accel_hzt (m, n) ndarray

2D array with the resultant horizontal acceleration series calculated at each cell. There is a 1D array with the resultant horizontal acceleration series at each cell (i, j).
Source code in pynewmarkdisp/spatial.py 
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def spatial_hzt_accel_from_2_dir(accel_ew, accel_ns, azimuth):
    """Calculate spatially the resultant acceleration from two hztl components.

    The azimuth is spatially distributed (e.g., obtained using the
    ``get_dip_azimuth`` function), but ``accel_ew``, ``accel_ns``, and
    ``time_vect`` are constant over all the spatial domain.

    Parameters
    ----------
    accel_ew : (p,) ndarray
        1D array with the acceleration series in the East-West direction.
    accel_ns : (p,) ndarray
        1D array with the acceleration series in the North-South direction.
    azimuth : (m, n) ndarray
        2D array with the spatial distribution of the azimuth of the slope dip
        direction, in [°] and measured clockwise from North.

    Returns
    -------
    accel_hzt : (m, n) ndarray
        2D array with the resultant horizontal acceleration series calculated
        at each cell. There is a 1D array with the resultant horizontal
        acceleration series at each cell (i, j).
    """
    row, col = azimuth.shape
    azimuth_r = np.radians(azimuth)
    accel_hzt = np.empty_like(azimuth_r, dtype=object)
    for i, j in product(range(row), range(col)):
        accel_hzt[i, j] = hzt_accel_from_2_dir(
            accel_ew, accel_ns, azimuth[i, j]
        )
    return accel_hzt

spatial_newmark 

spatial_newmark(time, accel, ky, g, azimuth=None)

Calculate the spatial distribution of permanent displacements.

Perform the direct Newmarkl method for each cell of the spatial domain and return the spatial distribution of the permanent displacements.

Parameters:

Name Type Description Default
time ndarray

1D or 2D array with the time data, in [s].

 required
accel tuple, list or ndarray

Object with the acceleration data. If it azimuth is given, accel must have a length of 2 and contain the acceleration in the East-West and North-South directions, respectively; both components of the acceleration might be constant or spatially distributed. If azimuth is None, the accel vector might be spatially distributed (2D array) or constant (1D array).

 required
ky ndarray

2D array with the spatial distribution of the critical seismic coefficient.

 required
g float

Acceleration of gravity, in the same units as accel. For example, if accel is in [m/s²], g should be equal to 9.81; if accel is as a fraction of gravity, g should be equal to 1.

 required
azimuth (m, n) ndarray or None

2D array with the spatial distribution of the azimuth of the slope dip direction, in [°] and measured clockwise from North. By default its value is None.

 None

Returns:

Name Type Description
permanent_disp (m, n) ndarray

2D array with the spatial distribution of permanent displacements.
Source code in pynewmarkdisp/spatial.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
def spatial_newmark(time, accel, ky, g, azimuth=None):
    """Calculate the spatial distribution of permanent displacements.

    Perform the direct Newmarkl method for each cell of the spatial domain and
    return the spatial distribution of the permanent displacements.

    Parameters
    ----------
    time : ndarray
        1D or 2D array with the time data, in [s].
    accel : tuple, list or ndarray
        Object with the acceleration data. If it ``azimuth`` is given,
        ``accel`` must have a length of 2 and contain the acceleration in the
        East-West and North-South directions, respectively; both components of
        the acceleration might be constant or spatially distributed. If
        ``azimuth`` is ``None``, the ``accel`` vector might be spatially
        distributed (2D array) or constant (1D array).
    ky : ndarray
        2D array with the spatial distribution of the critical seismic
        coefficient.
    g : float
        Acceleration of gravity, in the same units as ``accel``. For example,
        if ``accel`` is in [m/s²], ``g`` should be equal to 9.81; if ``accel``
        is as a fraction of gravity, ``g`` should be equal to 1.
    azimuth : (m, n) ndarray or None, optional
        2D array with the spatial distribution of the azimuth of the slope dip
        direction, in [°] and measured clockwise from North. By default its
        value is ``None``.

    Returns
    -------
    permanent_disp : (m, n) ndarray
        2D array with the spatial distribution of permanent displacements.
    """
    ky = num_to_2darray(ky)
    row, col = ky.shape
    permanent_disp = np.full(ky.shape, np.nan)
    for i, j in product(range(row), range(col)):
        if np.isnan(ky[i, j]):
            continue  # Skip step and continue the loop as u_p=np.nan
        time_s, accel_s = _select_case((i, j), time, accel, azimuth)
        newmark_str = direct_newmark(time_s, accel_s, ky[i, j], g)
        permanent_disp[i, j] = newmark_str["perm_disp"]
    return np.around(permanent_disp, 3)

verify_newmark_at_cell 

verify_newmark_at_cell(
    cell, time, accel, g, depth, depth_w, slope, phi, c, gamma, azimuth=None
)

Verify the direct Newmark method at a specific cell.

Parameters:

Name Type Description Default
cell tuple

Tuple with the index of the cell to be verified.

 required
time ndarray

1D or 2D array with the time data, in [s].

 required
accel tuple or list or ndarray

Object with the acceleration data. If it azimuth is given, accel must have a length of 2 and contain the acceleration in the East-West and North-South directions, respectively; both components of the acceleration might be constant or spatially distributed. If azimuth is None, the accel vector might be spatially distributed (2D array) or constant (1D array).

 required
g float

Acceleration of gravity, in the same units as accel. For example, if accel is in [m/s²], g should be equal to 9.81; if accel is as a fraction of gravity, g should be equal to 1.

 required
depth (m, n) ndarray

2D array with the spatial distribution of the potential sliding surface depth, in [m].

 required
depth_w ndarray

2D array with the spatial distribution of the water table depth, in [m].

 required
slope (m, n) ndarray

2D array with the spatial distribution of the slope inclination, in [°].

 required
phi (m, n) ndarray

2D array with the spatial distribution of the internal friction angle, in [°].

 required
c (m, n) ndarray

2D array with the spatial distribution of the cohesion, in [kPa].

 required
gamma ndarray

2D array with the spatial distribution of the soil unit weight, in [kN/m³].

 required
azimuth (m, n) ndarray or None

2D array with the spatial distribution of the azimuth of the slope dip direction, in [°] and measured clockwise from North. By default its value is None.

 None

Returns:

Name Type Description
newmark_str dict

Dictionary with the results of the direct Newmark method at the cell.
Source code in pynewmarkdisp/spatial.py 
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
def verify_newmark_at_cell(
    cell, time, accel, g, depth, depth_w, slope, phi, c, gamma, azimuth=None
):
    """Verify the direct Newmark method at a specific cell.

    Parameters
    ----------
    cell : tuple
        Tuple with the index of the cell to be verified.
    time : ndarray
        1D or 2D array with the time data, in [s].
    accel : tuple or list or ndarray
        Object with the acceleration data. If it ``azimuth`` is given,
        ``accel`` must have a length of 2 and contain the acceleration in the
        East-West and North-South directions, respectively; both components of
        the acceleration might be constant or spatially distributed. If
        ``azimuth`` is ``None``, the ``accel`` vector might be spatially
        distributed (2D array) or constant (1D array).
    g : float
        Acceleration of gravity, in the same units as ``accel``. For example,
        if ``accel`` is in [m/s²], ``g`` should be equal to 9.81; if ``accel``
        is as a fraction of gravity, ``g`` should be equal to 1.
    depth : (m, n) ndarray
        2D array with the spatial distribution of the potential sliding surface
        depth, in [m].
    depth_w : ndarray
        2D array with the spatial distribution of the water table depth, in
        [m].
    slope : (m, n) ndarray
        2D array with the spatial distribution of the slope inclination, in
        [°].
    phi : (m, n) ndarray
        2D array with the spatial distribution of the internal friction angle,
        in [°].
    c : (m, n) ndarray
        2D array with the spatial distribution of the cohesion, in [kPa].
    gamma : ndarray
        2D array with the spatial distribution of the soil unit weight, in
        [kN/m³].
    azimuth : (m, n) ndarray or None, optional
        2D array with the spatial distribution of the azimuth of the slope dip
        direction, in [°] and measured clockwise from North. By default its
        value is ``None``.

    Returns
    -------
    newmark_str : dict
        Dictionary with the results of the direct Newmark method at the cell.
    """
    i, j = cell
    depth, depth_w, slope = depth[i, j], depth_w[i, j], slope[i, j]
    phi, c, gamma = phi[i, j], c[i, j], gamma[i, j]
    fs_static = factor_of_safety(depth, depth_w, slope, phi, c, gamma, ks=0)
    ky = get_ky(depth, depth_w, slope, phi, c, gamma)
    fs_ky = factor_of_safety(depth, depth_w, slope, phi, c, gamma, ky)
    time_s, accel_s = _select_case((i, j), time, accel, azimuth)
    newmark_str = direct_newmark(time_s, accel_s, ky, g)
    newmark_str["fs_static"] = fs_static
    newmark_str["fs_ky"] = fs_ky
    return newmark_str

pynewmarkdisp.infslope

Functions to calculate the factor of safety of an infinite slope mechanism.

factor_of_safety 

factor_of_safety(depth, depth_w, slope, phi, c, gamma, ks=0)

Calculate the factor of safety for an infinite slope mechanism.

Inputs can be either constant or spatially distributed variables.

Parameters:

Name Type Description Default
depth float or (m, n) ndarray

Depth of the planar suface, in [m].

 required
depth_w float or ndarray

Depth of the watertable, in [m].

 required
slope float or ndarray

Slope inclination, in [°].

 required
phi float or ndarray

Slope inclination, in [°].

 required
c float or ndarray

Cohesion, in [kPa].

 required
gamma float or ndarray

Average unit weigth of the soil, in [kN/m³].

 required
ks float or ndarray

Seismic coefficient for pseudostatic analysis. Its default value is 0, meaning that the analysis is static (no seismic force included).

 0

Returns:

Name Type Description
fs float or ndarray

Factor of safety
Source code in pynewmarkdisp/infslope.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def factor_of_safety(depth, depth_w, slope, phi, c, gamma, ks=0):
    """Calculate the factor of safety for an infinite slope mechanism.

    Inputs can be either constant or spatially distributed variables.

    Parameters
    ----------
    depth : float or (m, n) ndarray
        Depth of the planar suface, in [m].
    depth_w : float or ndarray
        Depth of the watertable, in [m].
    slope : float or ndarray
        Slope inclination, in [°].
    phi : float or ndarray
        Slope inclination, in [°].
    c : float or ndarray
        Cohesion, in [kPa].
    gamma : float or ndarray
        Average unit weigth of the soil, in [kN/m³].
    ks : float or ndarray
        Seismic coefficient for pseudostatic analysis. Its default value is
        ``0``, meaning that the analysis is static (no seismic force included).

    Returns
    -------
    fs : float or ndarray
        Factor of safety
    """
    depth = num_to_2darray(depth)
    depth_w = num_to_2darray(depth_w)
    slope = num_to_2darray(np.radians(slope))  # Slope angle to radians
    phi = num_to_2darray(np.radians(phi))  # Friction angle to radians
    # Calculation
    # depth[np.where(depth == 0)] = np.nan  # Slice height cannot be zero
    head_pressure = (depth - depth_w) * np.cos(slope) ** 2
    head_pressure[np.where(depth_w > depth)] = 0  # Watertable below slip surf.
    pw = head_pressure * 9.81  # kN; Hydrostatic force
    w = depth * gamma * np.cos(slope)  # kPa; Weight of sliding mass
    tau = (w * (np.cos(slope) - ks * np.sin(slope)) - pw) * np.tan(phi) + c
    slip_force = w * (np.sin(slope) + ks * np.cos(slope))
    # slip_force[np.where(beta == 0)] = 0.001  # ~zero for horizontal surfaces
    np.seterr(divide="ignore")
    fs = np.divide(tau, slip_force)
    fs[np.where(fs > 3)] = 3  # Factor of safety restricted to max. 3
    if fs.size == 1:
        fs = fs[0, 0]
    return np.around(fs, 3)

get_ky 

get_ky(depth, depth_w, slope, phi, c, gamma)

Compute the critical seismic coefficient of an infinite slope mechanism.

Inputs can be either constant or spatially distributed variables.

Parameters:

Name Type Description Default
depth float or (m, n) ndarray

Depth of the planar suface, in [m].

 required
depth_w float or ndarray

Depth of the watertable, in [m].

 required
slope float or ndarray

Slope inclination, in [°].

 required
phi float or ndarray

Slope inclination, in [°].

 required
c float or ndarray

Cohesion, in [kPa].

 required
gamma float or ndarray

Average unit weigth of the soil, in [kN/m³].

 required

Returns:

Name Type Description
ky float or ndarray

Critical seismic coefficient.
Source code in pynewmarkdisp/infslope.py 
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def get_ky(depth, depth_w, slope, phi, c, gamma):
    """Compute the critical seismic coefficient of an infinite slope mechanism.

    Inputs can be either constant or spatially distributed variables.

    Parameters
    ----------
    depth : float or (m, n) ndarray
        Depth of the planar suface, in [m].
    depth_w : float or ndarray
        Depth of the watertable, in [m].
    slope : float or ndarray
        Slope inclination, in [°].
    phi : float or ndarray
        Slope inclination, in [°].
    c : float or ndarray
        Cohesion, in [kPa].
    gamma : float or ndarray
        Average unit weigth of the soil, in [kN/m³].

    Returns
    -------
    ky : float or ndarray
        Critical seismic coefficient.
    """
    depth = num_to_2darray(depth)
    depth_w = num_to_2darray(depth_w)
    slope = np.radians(slope)  # Slope angle to radians
    phi = np.radians(phi)  # Friction angle to radians
    # Intermediate calculations
    # depth[np.where(depth == 0)] = np.nan  # Slice height cannot be zero
    head_pressure = (depth - depth_w) * np.cos(slope) ** 2
    pw = head_pressure * 9.81  # kN; Hydrostatic force
    w = depth * gamma * np.cos(slope)  # kPa; Weight of sliding mass
    # Ky calculation
    np.seterr(divide="ignore")
    # num = np.divide(np.cos(beta) - pw, w) * np.divide(
    #     np.tan(phi) + c, w
    # ) - np.sin(beta)
    w[np.where(w == 0)] = 0.001
    num = (np.cos(slope) - pw / w) * np.tan(phi) + c / w - np.sin(slope)
    den = np.cos(slope) + np.sin(slope) * np.tan(phi)
    ky = np.divide(num, den)
    # ky = num / den
    ky[np.where(ky > 3)] = 3  # The slope is stable under pseudo-static cond.
    ky[np.where(ky < 0)] = np.nan  # The slope is unstable under static cond.
    # ky[np.where(ky > 1)] = 1  # The slope is stable under sedudostatic cond.
    if ky.size == 1:
        ky = ky[0, 0]
    return np.around(ky, 3)

min_fs_and_depth 

min_fs_and_depth(depth, depth_w, slope, phi, c, gamma, ks=0, nd=10)

Calculate the minimum factor of safety and the depth at which it occurs.

As for the factor_of_safety function, inputs can be either constant or spatially distributed variables.

It evaluates the factor of safety for a range of nd depths between the terrain surface and the slip surface and return the minimum factor of safety and the depth at which it occurs. This is useful to find the critical slip surface, which will be located at the first depth where the factor of safety drops below 1 or where it reaches its minimum value.

Parameters:

Name Type Description Default
depth float or (m, n) ndarray

Depth of the planar suface, in [m].

 required
depth_w float or ndarray

Depth of the watertable, in [m].

 required
slope float or ndarray

Slope inclination, in [°].

 required
phi float or ndarray

Slope inclination, in [°].

 required
c float or ndarray

Cohesion, in [kPa].

 required
gamma float or ndarray

Average unit weigth of the soil, in [kN/m³].

 required
ks float or ndarray

Seismic coefficient for pseudostatic analysis. Its default value is 0, meaning that the analysis is static (no seismic force included).

 0
nd int

Number of depths to evaluate between the terrain surface and the slip surface. The default is 10.

 10

Returns:

Name Type Description
min_fs float or ndarray

Minimum factor of safety
depth_at_min_fs float or ndarray

Depth at which the minimum factor of safety occurs, in [m].
Source code in pynewmarkdisp/infslope.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def min_fs_and_depth(depth, depth_w, slope, phi, c, gamma, ks=0, nd=10):
    """Calculate the minimum factor of safety and the depth at which it occurs.

    As for the ``factor_of_safety`` function, inputs can be either constant or
    spatially distributed variables.

    It evaluates the factor of safety for a range of ``nd`` depths between the
    terrain surface and the slip surface and return the minimum factor of
    safety and the depth at which it occurs. This is useful to find the
    critical slip surface, which will be located at the first depth where the
    factor of safety drops below 1 or where it reaches its minimum value.

    Parameters
    ----------
    depth : float or (m, n) ndarray
        Depth of the planar suface, in [m].
    depth_w : float or ndarray
        Depth of the watertable, in [m].
    slope : float or ndarray
        Slope inclination, in [°].
    phi : float or ndarray
        Slope inclination, in [°].
    c : float or ndarray
        Cohesion, in [kPa].
    gamma : float or ndarray
        Average unit weigth of the soil, in [kN/m³].
    ks : float or ndarray
        Seismic coefficient for pseudostatic analysis. Its default value is
        ``0``, meaning that the analysis is static (no seismic force included).
    nd : int, optional
        Number of depths to evaluate between the terrain surface and the slip
        surface. The default is ``10``.

    Returns
    -------
    min_fs : float or ndarray
        Minimum factor of safety
    depth_at_min_fs : float or ndarray
        Depth at which the minimum factor of safety occurs, in [m].
    """
    depth = num_to_2darray(depth)
    # depth[np.where(depth == 0)] = np.nan  # Slice height cannot be zero
    min_fs = factor_of_safety(depth, depth_w, slope, phi, c, gamma, ks)
    depth_at_min_fs = depth.copy()
    for ith_depth_fract in np.linspace(0, 1, nd + 1)[1:]:
        ith_depth = ith_depth_fract * depth
        ith_fs = factor_of_safety(ith_depth, depth_w, slope, phi, c, gamma, ks)
        mask = np.where(ith_fs < min_fs)
        min_fs[mask] = ith_fs[mask]
        depth_at_min_fs[mask] = ith_depth[mask]
    return min_fs, depth_at_min_fs

num_to_2darray 

num_to_2darray(num)

Turn an int or float input into a 2D array with shape (1, 1).

Parameters:

Name Type Description Default
num int or float

Variable to be converted to a 2D array

 required

Returns:

Name Type Description
num_array ndarray

2D array with shape (1, 1) whose only element is the input variable.
Source code in pynewmarkdisp/infslope.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def num_to_2darray(num):
    """Turn an int or float input into a 2D array with shape (1, 1).

    Parameters
    ----------
    num : int or float
        Variable to be converted to a 2D array

    Returns
    -------
    num_array : ndarray
        2D array with shape (1, 1) whose only element is the input variable.
    """
    return (
        np.array(num).reshape((1, 1)) if isinstance(num, (int, float)) else num
    )

pynewmarkdisp.empir_corr

Functions to calculate permanent displacements from empirical correlations.

Most of the functions in this module are taken from Table 1 in Meehan & Vahedifard (2013), who compiled a several empirical correlations or simplified methods for calculating the permanent displacement of a slope. The keys of the correlations dictionary are the same as the numbering of the cited table. However, additional functions were added to the dictionary after consulting the original references.

ambraseys_and_menu_88 

ambraseys_and_menu_88(ay, a_max, **kwargs)

Ambraseys and Menu (1988) empirical correlation for uₚ.

Eq. 10 in Meehan & Vahedifard (2013).

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def ambraseys_and_menu_88(ay, a_max, **kwargs):
    """Ambraseys and Menu (1988) empirical correlation for uₚ.

    Eq. 10 in Meehan & Vahedifard (2013).

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * 10**0.90 * ((1 - car) ** 2.53 * car ** (-1.09))

bray_and_travasarou_07_eq6 

bray_and_travasarou_07_eq6(ay, a_max, magnitude, **kwargs)

Bray and Travasarou (2007) empirical correlation for uₚ.

Eq. 10 in Meehan & Vahedifard (2013). Eq. 6 in Bray and Travasarou (2007).

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required
magnitude float

Earthquake magnitude.

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def bray_and_travasarou_07_eq6(ay, a_max, magnitude, **kwargs):
    """Bray and Travasarou (2007) empirical correlation for uₚ.

    Eq. 10 in Meehan & Vahedifard (2013). Eq. 6 in Bray and Travasarou (2007).

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.
    magnitude : float
        Earthquake magnitude.

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * np.exp(
        -0.22
        - 2.83 * np.log(ay)
        - 0.333 * (np.log(ay)) ** 2
        + 0.566 * np.log(ay) * np.log(a_max)
        + 3.04 * np.log(a_max)
        - 0.244 * (np.log(a_max)) ** 2
        + 0.278 * (magnitude - 7)
    )

hynesgriffin_and_Franklin_1984_mean 

hynesgriffin_and_Franklin_1984_mean(ay, a_max, **kwargs)

Hynes-Griffin and Franklin (1984) empirical correlation for uₚ (mean).

Eq. 8 in Meehan & Vahedifard (2013). It is a functional form of the mean resulting from regression analysis of chart-based solution from the original paper.

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def hynesgriffin_and_Franklin_1984_mean(ay, a_max, **kwargs):
    """Hynes-Griffin and Franklin (1984) empirical correlation for uₚ (mean).

    Eq. 8 in Meehan & Vahedifard (2013). It is a functional form of the mean
    resulting from regression analysis of chart-based solution from the
    original paper.

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * 10 ** (
        -0.287
        - 2.854 * car
        - 1.733 * car**2
        - 0.702 * car**3
        - 0.116 * car**4
    )

hynesgriffin_and_Franklin_1984_ub 

hynesgriffin_and_Franklin_1984_ub(ay, a_max, **kwargs)

Hynes-Griffin and Franklin (1984) empirical corr. for uₚ (upper bound).

Eq. 7 in Meehan & Vahedifard (2013). It is a functional form of the upper bound resulting from regression analysis of chart-based solution from the original paper.

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def hynesgriffin_and_Franklin_1984_ub(ay, a_max, **kwargs):
    """Hynes-Griffin and Franklin (1984) empirical corr. for uₚ (upper bound).

    Eq. 7 in Meehan & Vahedifard (2013). It is a functional form of the upper
    bound resulting from regression analysis of chart-based solution from the
    original paper.

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * 10 ** (
        0.804
        - 1.847 * car
        - 0.285 * car**2
        + 0.193 * car**3
        + 0.078 * car**4
    )

jibson_07_eq10 

jibson_07_eq10(ay, a_max, arias_int, **kwargs)

Jibson (2007) empirical correlation for uₚ.

Eq. 19 in Meehan & Vahedifard (2013). Eq. 10 in Jibson (2007).

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required
arias_int float

Arias intensity, in [m/s].

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def jibson_07_eq10(ay, a_max, arias_int, **kwargs):
    """Jibson (2007) empirical correlation for uₚ.

    Eq. 19 in Meehan & Vahedifard (2013). Eq. 10 in Jibson (2007).

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.
    arias_int : float
        Arias intensity, in [m/s].

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * 10 ** (
        0.561 * np.log10(arias_int) - 3.833 * np.log10(car) - 1.474
    )

jibson_07_eq6 

jibson_07_eq6(ay, a_max, **kwargs)

Jibson (2007) empirical correlation for uₚ.

Eq. 16 in Meehan & Vahedifard (2013). Eq. 6 in Jibson (2007).

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
def jibson_07_eq6(ay, a_max, **kwargs):
    """Jibson (2007) empirical correlation for uₚ.

    Eq. 16 in Meehan & Vahedifard (2013). Eq. 6 in Jibson (2007).

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * 10 ** (0.215) * ((1 - ay / a_max) ** 2.341 * car ** (-1.438))

jibson_07_eq9 

jibson_07_eq9(ay, arias_int, **kwargs)

Jibson (2007) empirical correlation for uₚ.

Eq. 18 in Meehan & Vahedifard (2013). Eq. 9 in Jibson (2007).

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
arias_int float

Arias intensity, in [m/s].

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
def jibson_07_eq9(ay, arias_int, **kwargs):
    """Jibson (2007) empirical correlation for uₚ.

    Eq. 18 in Meehan & Vahedifard (2013). Eq. 9 in Jibson (2007).

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    arias_int : float
        Arias intensity, in [m/s].

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    return 1e-2 * 10 ** (
        2.401 * np.log10(arias_int) - 3.481 * np.log10(ay) - 3.23
    )

saygili_and_rathje_08_eq5 

saygili_and_rathje_08_eq5(ay, a_max, **kwargs)

Saygili and Rathje (2008) empirical correlation for uₚ.

Eq. 20 in Meehan & Vahedifard (2013). Eq. 5 in Saygili and Rathje (2008).

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
def saygili_and_rathje_08_eq5(ay, a_max, **kwargs):
    """Saygili and Rathje (2008) empirical correlation for uₚ.

    Eq. 20 in Meehan & Vahedifard (2013). Eq. 5 in Saygili and Rathje (2008).

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * np.exp(
        5.52
        - 4.43 * car
        - 20.93 * car**2
        + 42.61 * car**3
        - 28.74 * car**4
        + 0.72 * np.log(a_max)
    )

saygili_and_rathje_08_eq6col2 

saygili_and_rathje_08_eq6col2(ay, a_max, v_max, **kwargs)

Saygili and Rathje (2008) empirical correlation for uₚ.

Eq. 21 in Meehan & Vahedifard (2013). Eq. 6 in Saygili and Rathje (2008) with parameters from Table 1, column 2.

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required
v_max float

Peak horizontal ground velocity, in [cm/s].

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
def saygili_and_rathje_08_eq6col2(ay, a_max, v_max, **kwargs):
    """Saygili and Rathje (2008) empirical correlation for uₚ.

    Eq. 21 in Meehan & Vahedifard (2013). Eq. 6 in Saygili and Rathje (2008)
    with parameters from Table 1, column 2.

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.
    v_max : float
        Peak horizontal ground velocity, in [cm/s].

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * np.exp(
        -1.56
        - 4.58 * car
        - 20.84 * car**2
        + 44.75 * car**3
        - 30.5 * car**4
        - 0.64 * np.log(a_max)
        + 1.55 * np.log(v_max)
    )

saygili_and_rathje_08_eq6col4 

saygili_and_rathje_08_eq6col4(ay, a_max, arias_int, **kwargs)

Saygili and Rathje (2008) empirical correlation for uₚ.

Not in Meehan & Vahedifard (2013), but it is equivalent to Eq. 21. Eq. 6 in Saygili and Rathje (2008) with parameters from Table 1, column 4.

Parameters:

Name Type Description Default
ay float

Critical acceleration, in fraction of g.

 required
a_max float

Peak horizontal ground acceleration, in fraction of g.

 required
arias_int float

Arias intensity, in [m/s].

 required

Returns:

Name Type Description
perm_disp float

Permanent displacement, in [m].
Source code in pynewmarkdisp/empir_corr.py 
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
def saygili_and_rathje_08_eq6col4(ay, a_max, arias_int, **kwargs):
    """Saygili and Rathje (2008) empirical correlation for uₚ.

    Not in Meehan & Vahedifard (2013), but it is equivalent to Eq. 21.
    Eq. 6 in Saygili and Rathje (2008) with parameters from Table 1, column 4.

    Parameters
    ----------
    ay : float
        Critical acceleration, in fraction of g.
    a_max : float
        Peak horizontal ground acceleration, in fraction of g.
    arias_int : float
        Arias intensity, in [m/s].

    Returns
    -------
    perm_disp : float
        Permanent displacement, in [m].
    """
    # ay[np.where(ay > a_max)] = a_max  # There is no excess acceleration
    car = ay / a_max  # Critical acceleration ratio
    return 1e-2 * np.exp(
        2.39
        - 5.24 * car
        - 18.78 * car**2
        + 42.01 * car**3
        - 29.15 * car**4
        - 1.56 * np.log(a_max)
        + 1.38 * np.log(arias_int)
    )