# -*- coding: utf-8 -*-
'''
Module for defining the classes related to the slip surface, either a circular
or a composite geometry (*e.g.* a tortuous failure surface)
'''
[docs]class CircularSurface:
'''Creates an instance of an object that defines the structure of an
polyline which represents the slip surface of a landslide which geometry
is a circumference-arc. ::
CircularSurface(slopeCoords, dist1, dist2, radius, concave=True)
The arc is defined with two points on the terrain surface and the radius.
That implies there are two possible solutions; to select which one is
wanted, it is necessary to modify the variable ``concave``.
It is possible the arc cuts across the terrain surface in some point
different to its ends, perhaps because of some swedge in the terrrain or
the radius is too long. In that case, the method ``defineStructre`` changes
the attribute ``dist2`` such that it is replaced by the horizontal distance
of the intersection point.
Attributes:
slopeCoords ((2, n `numpy.ndarray`): Coordinates of the vertices
of the polygon within which the slope mass is defined. It
is obtained with the method ``defineboundary`` either from the
classes ``AnthropicSlope`` or ``NaturalSlope`` (module ``slope``).
dist1 (`int` or `float`): First horizontal distance from the leftmost
point of the terrain surface (including the crown) where the arc
is intersected with it.
dist2 (`int` or `float`): Second horizontal distance from the leftmost
point of the terrain surface (including the crown) where the arc
is intersected with it.
radius (`int` or `float`): Length of the circumference-arc radius.
concave (`bool`): Logical variable to define if it is wanted that
the circumference-arc will be concave (upwards), otherwise, it will
be convexe (downwards). Default value is ``True``.
Note:
The class ``CircularSurface`` requires
`numpy <http://www.numpy.org/>`_,
`matplotlib <https://matplotlib.org/>`_ and
`shapely <https://pypi.python.org/pypi/Shapely>`_.
Examples:
>>> from numpy import array
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.slipsurface import CircularSurface
>>> slope = AnthropicSlope(slopeHeight=7.5, slopeDip=[1, 1.5],
>>> crownDist=5, toeDist=5)
>>> surface = CircularSurface(slopeCoords=slope.coords,
>>> dist1=2, dist2=10, radius=9)
>>> surface.__dict__.keys()
dict_keys(['slopeCoords', 'dist1', 'dist2', 'radius', 'concave',
'point1', 'point2', 'center', 'initAngle',
'endAngle', 'coords'])
'''
def __init__(self, slopeCoords, dist1, dist2, radius, concave=True):
'''
CircularSurface(slopeCoords, dist1, dist2, radius, concave=True)
'''
self.slopeCoords = slopeCoords
self.dist1 = min(dist1, dist2)
self.dist2 = max(dist1, dist2)
if self.dist2 > max(slopeCoords[0]):
self.dist2 = max(slopeCoords[0])
self.radius = radius
self.concave = concave
# Definde structure of the arc
self.defineStructre()
[docs] def defineStructre(self):
'''Method to define the structure of the circumference-arc which
represents the slip surface of a landslide.
If the arc cuts across the terrain surface in some point different to
its ends, the attribute ``dist2`` is modified to the horizontal
distance of the intersection point.
The returned angles have values betwen
:math: `\\left[\\pi, -\\pi \\right)`, where the angle equal to zero
coincides with the vector :math: `\\left(1, 0) \\right)`.
Returns:
(`dict`): dictionary with the following outputs.
- **center** (`tuple`): Coordinates of the circumference-arc\
center
- **endAngle** (`float`): Angle in radians of the vector\
that points from the center to the first intersection\
between the terrain surface and the circumference-arc.
- **initAngle** (`float`): Angle in radians of the vector\
that points from the center to the first intersection\
between the terrain surface and the circumference-arc.
- **point1** (`tuple`): Coordinates of the first point that\
intersects the terrain surface
- **point2** (`tuple`): Coordinates of the second point that\
intersects the terrain surface
Examples:
>>> from numpy import array
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.slipsurface import CircularSurface
>>> slope = AnthropicSlope(slopeHeight=7.5, slopeDip=[1, 1.5],
>>> crownDist=5, toeDist=5)
>>> surface = CircularSurface(slopeCoords=slope.coords,
>>> dist1=2, dist2=10, radius=9)
>>> surface.defineStructre()
{'center': (10.881322862689261, 10.831744386868543),
'endAngle': -1.668878272858519,
'initAngle': -2.979016942655663,
'point1': array([2. , 9.375]),
'point2': array([10. , 1.875])}
>>> from numpy import array
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.slipsurface import CircularSurface
>>> slope = AnthropicSlope(slopeHeight=7.5, slopeDip=[1, 1.5],
>>> crownDist=5, toeDist=5)
>>> surface = CircularSurface(slopeCoords=slope.coords,
>>> dist1=2, dist2=10, radius=1)
>>> surface.defineStructre()
ValueError: separation of points > diameter
>>> from numpy import array
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.slipsurface import CircularSurface
>>> slope = AnthropicSlope(slopeHeight=7.5, slopeDip=[1, 1.5],
>>> crownDist=5, toeDist=5)
>>> surface = CircularSurface(slopeCoords=slope.coords,
>>> dist1=2, dist2=10, radius=6)
>>> surface.defineStructre()
ValueError: Radius too short. Increase at least 1.516
'''
import numpy as np
from shapely.geometry import LineString
# Getting the arc coordinates on the terrain surface
if self.slopeCoords[0, 1] > self.slopeCoords[0, -2]:
self.slopeCoords = np.fliplr(self.slopeCoords)
terrainSurfLS = LineString(self.slopeCoords[:, 1:-2].T)
yMin, yMax = self.slopeCoords[1].min()/2, self.slopeCoords[1].max()*2
for name, dist in [('point1', self.dist1),
('point2', self.dist2)]:
vertLine = LineString([(dist, yMin), (dist, yMax)])
intersection = terrainSurfLS.intersection(vertLine)
setattr(self, name, np.array([intersection.x, intersection.y]))
# delta x, delta y between points
dx, dy = self.point2 - self.point1
# dist between points
pointSep = np.linalg.norm(self.point2 - self.point1)
# Minimum radius
minRadius = round(0.5 * pointSep**2 / abs(dx), 3)
# Verification of minimum distance between both points
if pointSep > 2*self.radius:
raise ValueError('separation of points > diameter')
# Verification of minimum radius
elif minRadius > self.radius:
raise ValueError('Radius too short. Increase at least {}'.format(
round(minRadius - self.radius, 3)))
# halfway point
xHalfPoint = (self.point1[0] + self.point2[0])/2
yHalfPoint = (self.point1[1] + self.point2[1])/2
# distance along the mirror line
d = np.sqrt(self.radius**2 - (0.5*pointSep)**2)
# Verification of the minimum radius
if self.concave:
center = (xHalfPoint - d*dy/pointSep, yHalfPoint + d*dx/pointSep)
else:
center = (xHalfPoint + d*dy/pointSep, yHalfPoint - d*dx/pointSep)
setattr(self, 'center', center)
# Getting the angles from the center to the extreme points of the arc
vectInters1 = self.point1 - center
vectInters2 = self.point2 - center
angles = np.arctan2([vectInters1[1], vectInters2[1]],
[vectInters1[0], vectInters2[0]])
setattr(self, 'initAngle', angles[0])
setattr(self, 'endAngle', angles[1])
# Getting the x and y coordinates of 100 points in the arc
anglAux = angles[:]
if self.concave:
anglAux = np.linspace(self.initAngle % (2*np.pi),
self.endAngle % (2*np.pi), 100)
xC = [self.center[0] + self.radius*np.cos(theta) for theta in anglAux]
yC = [self.center[1] + self.radius*np.sin(theta) for theta in anglAux]
setattr(self, 'coords', np.array([xC, yC]))
# Verify if the arc intersects the terrain in some point diferent to
# the initial extreme points. If true, cut it in the intersection-point
lineArc = LineString(self.coords[:, 1:-1].T)
intersection = terrainSurfLS.intersection(lineArc)
if intersection.geom_type == 'Point': # just one intersection
self.dist2 = intersection.x
self.defineStructre()
elif intersection.geom_type == 'MultiPoint': # more than one
self.dist2 = intersection[0].x
self.defineStructre()
return {'center': center, 'point1': self.point1, 'point2': self.point2,
'initAngle': angles[0], 'endAngle': angles[1]}
[docs] def plot(self):
'''Method for generating a graphic of the circumference-arc and the
slope.
Returns:
(`matplotlib.figure.Figure`): object with the matplotlib structure\
of the plot. You might use it to save the figure for example.
Examples:
>>> from numpy import array
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.slipsurface import CircularSurface
>>> slope = AnthropicSlope(slopeHeight=7.5, slopeDip=[1, 1.5],
>>> crownDist=5, toeDist=5)
>>> surface = CircularSurface(slopeCoords=slope.coords,
>>> dist1=2, dist2=10, radius=9)
>>> fig = surface.plot()
.. figure:: https://rawgit.com/eamontoyaa/pybimstab/master/examples/figures/slipsurface_CircularSurface_example1.svg
:alt: slipsurface_CircularSurface_example1
.. only:: html
:download:`example script<../examples/figuresScripts/slipsurface_CircularSurface_example1.py>`.
>>> from numpy import array
>>> from pybimstab.slope import NaturalSlope
>>> from pybimstab.slipsurface import CircularSurface
>>> terrainCoords = array(
>>> [[-2.49, 0.1 , 1.7 , 3.89, 5.9 , 8.12, 9.87, 13.29, 20.29,
>>> 21.43, 22.28, 23.48, 24.65, 25.17],
>>> [18.16, 17.88, 17.28, 15.73, 14.31, 13.58, 13, 3.61, 3.61,
>>> 3.32, 2.71, 2.23, 1.21, 0.25]])
>>> slope = NaturalSlope(terrainCoords)
>>> surface = CircularSurface(slopeCoords=slope.coords,
>>> dist1=7, dist2=20, radius=13)
>>> fig = surface.plot()
.. plot::
.. figure:: https://rawgit.com/eamontoyaa/pybimstab/master/examples/figures/slipsurface_CircularSurface_example2.svg
:alt: slipsurface_CircularSurface_example2
.. only:: html
:download:`example script<../examples/figuresScripts/slipsurface_CircularSurface_example2.py>`.
'''
from matplotlib import pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
# Plot body
ax.plot(self.slopeCoords[0], self.slopeCoords[1], '-k')
ax.plot(self.coords[0], self.coords[1], '-r')
ax.plot(self.center[0], self.center[1], '.r')
# Plot settings
ax.grid(True, ls='--', lw=0.5)
ax.set_aspect(1)
fig.tight_layout()
return fig
# %%
[docs]class TortuousSurface:
'''Creates an instance of an object that defines the structure of an
polyline which represents the slip surface of a landslide which geometry
is a tortuous path surrouding the blocks inside the slope mass. ::
TortuousSurface(bim, dist1, dist2, heuristic='Manhattan',
reverseLeft=False, reverseUp=False, smoothFactor=0,
preferredPath=None, prefPathFact=None)
The surface is defined with two points on the terrain surface and the
heuristic function. It is possible to set a forced path to modify the free
trajectory of the tortuous path with the aim of move it closer to, for
example a circular surface.
Attributes:
bim (`BlocksInMatrix` object): object with the structure of the slope
made of the Blocks-In-Matrix material.
dist1 (`int` or `float`): First horizontal distance from the leftmost
point of the terrain surface (including the crown) where the arc
is intersected with it.
dist2 (`int` or `float`): Second horizontal distance from the leftmost
point of the terrain surface (including the crown) where the arc
is intersected with it.
heuristic (`str`): Name of the geometric model to determine the
heuristic distance. It must be selected either ``Manhattan`` or
``Euclidean``; their description can be found in the ``Astar``
class documentation. `Manhattan` is the default value.
reverseLeft (`bool`): Logical variable to allow or not reverses
movements to the left. Default value is ``False``.
reverseUp (`bool`): Logical variable to allow or not reverses
movements to upward. Default value is ``False``.
smoothFactor (`int`): Value to indicate the B-spline interpolation
order of the smooter function. If is equal to zero, which is the
default value, the surface will not be smoothed.
preferredPath (`numpy.ndarray` or `None`): (2, n) array with the
coordinates of a path where the tortuous surface is going to be
forced; ``None`` is the default value.
prefPathFact (`int` or `float` or `None`): Multiplier of the shortest
distance between the current point and the polyline; ``None`` is
the default value.
Note:
The class ``TortuousSurface`` requires
`numpy <http://www.numpy.org/>`_,
`matplotlib <https://matplotlib.org/>`_ and
`shapely <https://pypi.python.org/pypi/Shapely>`_.
Examples:
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.bim import BlocksInMatrix
>>> from pybimstab.slipsurface import TortuousSurface
>>> slope = AnthropicSlope(slopeHeight=12, slopeDip=[1, 1.5],
>>> crownDist=10, toeDist=10)
>>> bim = BlocksInMatrix(slopeCoords=slope.coords, blockProp=0.25,
>>> tileSize=0.25, seed=123)
>>> surface = TortuousSurface(
>>> bim, dist1=0, dist2=17, heuristic='manhattan',
>>> reverseLeft=False, reverseUp=False, smoothFactor=0,
>>> preferredPath=None, prefPathFact=None)
>>> surface.__dict__.keys()
dict_keys(['bim', 'dist1', 'dist2', 'heuristic', 'reverseLeft',
'reverseUp', 'smoothFactor', 'preferredPath',
'prefPathFact', 'terrainSurfLS', 'point1', 'end1', 'point2',
'end2', 'startIdx', 'goalIdx', 'coords'])
'''
def __init__(self, bim, dist1, dist2, heuristic='manhattan',
reverseLeft=False, reverseUp=False, smoothFactor=0,
preferredPath=None, prefPathFact=None):
'''
TortuousSurface(bim, dist1, dist2, heuristic='Manhattan',
reverseLeft=False, reverseUp=False, smoothFactor=0,
preferredPath=None, prefPathFact=None)
'''
self.bim = bim
self.dist1 = min(dist1, dist2)
self.dist2 = max(dist1, dist2)
if self.dist2 > max(bim.slopeCoords[0]):
self.dist2 = max(bim.slopeCoords[0])
self.heuristic = heuristic
self.reverseLeft = reverseLeft
self.reverseUp = reverseUp
self.smoothFactor = smoothFactor
self.preferredPath = preferredPath
if preferredPath is not None and prefPathFact is None:
prefPathFact = 1
self.prefPathFact = prefPathFact
# Obtain the indexes at the ends of the slip surface
self.getIndexesAtEnds()
# Moving the indexes of the ends when are blocks
if self.bim.grid[self.startIdx] == 1:
self.dist1 += self.bim.tileSize
self.getIndexesAtEnds()
if self.bim.grid[self.goalIdx] == 1:
self.dist2 = self.dist2 - self.bim.tileSize
self.getIndexesAtEnds()
# Obtain the optimum path through the A* algorithm
self.defineStructre()
[docs] def getIndexes(self, coord):
'''Method for obtaining the array indexes of the BIM structure for a
coordinate given in the real scale of the slope stability problem.
The transformation is performed by rounding the division between the
coordinate and the tile size with the ``int_`` function of `numpy`.
That means that always rounds to the left and bottom sides of a tile.
Attributes:
coord (`tuple`): Coordinates of some point in the slope mass or
surface, which is wanted to get them indexes into the BIM
grid-grapth structure.
Returns:
(`tuple`): Indexes of a coordinate from the real-scale problem\
projected to the array that represents the BIM structure of\
the slope; the first value of the tuple is the row and the\
second one is the column of the grid-array respectively.
Examples:
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.bim import BlocksInMatrix
>>> from pybimstab.slipsurface import TortuousSurface
>>> slope = AnthropicSlope(slopeHeight=12, slopeDip=[1, 1.5],
>>> crownDist=10, toeDist=10)
>>> bim = BlocksInMatrix(slopeCoords=slope.coords, blockProp=0.25,
>>> tileSize=0.25, seed=123)
>>> surface = TortuousSurface(bim, dist1=0, dist2=17)
>>> surface.getIndexes(surface.point1)
(66, 00)
>>> surface.getIndexes(surface.point2)
(24, 68)
'''
import numpy as np
xIdx, yIdx = np.int_(np.array(coord) / self.bim.tileSize)
return yIdx, xIdx
[docs] def getCoord(self, indexes):
'''Method for obtaining the real scale problem coordinates of of some
cell in the BIM grid-graph structure of the slope.
The transformation is performed by getting the center of the tile which
contains the coordinates of the point.
Attributes:
indexes (`tuple`): Indexes of some cell in the BIM grid-graph
structure of the slope; the first tuple-value is the ordinate
and the second one is the abscisse.
Returns:
(`tuple`): Indexes of a coordinate from the real-scale problem\
projected to the array that represents the BIM structure of\
the slope; the first tuple value is the abscisse and the\
second one is the ordinate.
Examples:
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.bim import BlocksInMatrix
>>> from pybimstab.slipsurface import TortuousSurface
>>> slope = AnthropicSlope(slopeHeight=12, slopeDip=[1, 1.5],
>>> crownDist=10, toeDist=10)
>>> bim = BlocksInMatrix(slopeCoords=slope.coords, blockProp=0.25,
>>> tileSize=0.25, seed=123)
>>> surface = TortuousSurface(bim, dist1=0, dist2=17)
>>> surface.getCoord((66, 0))
(0.125, 16.446428571428573)
>>> surface.getCoord((24, 68))
(17.125, 5.946428571428573)
'''
x = self.bim.xCells[indexes] + 0.5*self.bim.tileSize
y = self.bim.yCells[indexes] + 0.5*self.bim.tileSize
return x, y
[docs] def getIndexesAtEnds(self):
'''Method for obtaining the array indexes of the BIM grid-grapth
structure for the ends of the slip surface.
Examples:
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.bim import BlocksInMatrix
>>> from pybimstab.slipsurface import TortuousSurface
>>> slope = AnthropicSlope(slopeHeight=12, slopeDip=[1, 1.5],
>>> crownDist=10, toeDist=10)
>>> bim = BlocksInMatrix(slopeCoords=slope.coords, blockProp=0.25,
>>> tileSize=0.25, seed=123)
>>> surface = TortuousSurface(bim, dist1=0, dist2=17)
>>> surface.getIndexesAtEnds()
((66, 0), (23, 68))
'''
import numpy as np
from shapely.geometry import LineString
# Getting the coordinates of the ends
if self.bim.slopeCoords[0, 1] > self.bim.slopeCoords[0, -2]:
self.bim.slopeCoords = np.fliplr(self.bim.slopeCoords)
terrainSurfLS = LineString(self.bim.slopeCoords[:, 1:-2].T)
setattr(self, 'terrainSurfLS', terrainSurfLS)
yMin = self.bim.slopeCoords[1].min() / 2
yMax = self.bim.slopeCoords[1].max() * 2
tileSize = self.bim.tileSize
for name, dist in [('point1', self.dist1),
('end1', self.dist1 - (self.dist1 % tileSize)),
('point2', self.dist2),
('end2', self.dist2 - (self.dist2 % tileSize) +
tileSize)]:
vertLine = LineString([(dist, yMin), (dist, yMax)])
intersection = terrainSurfLS.intersection(vertLine)
setattr(self, name, np.array([intersection.x, intersection.y]))
# Indexes of the start and goal nodes
startIdx = self.getIndexes(self.point1)
try: # Controling when index is out the grid dimension
self.bim.grid[startIdx]
except Exception:
startIdx = (startIdx[0]-1, startIdx[1])
goalIdx = self.getIndexes(self.point2)
# Moving the indexes of the ends when are out the slope
while self.bim.grid[startIdx] == -1 or self.bim.grid[goalIdx] == -1:
# Start index
if self.bim.grid[startIdx] == -1:
startIdx = (startIdx[0]-1, startIdx[1])
if self.bim.grid[goalIdx] == -1:
goalIdx = (goalIdx[0]-1, goalIdx[1])
setattr(self, 'startIdx', startIdx)
setattr(self, 'goalIdx', goalIdx)
return startIdx, goalIdx
[docs] def defineStructre(self):
'''Method to define the structure of the tortuous path which represents
the slip surface of a landslide that occurs in a slope made of BIM.
The surface is generated through the :math:`\\mathrm{A}^\\ast`
algorithm defined in the ``Astar`` module.
Examples:
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.bim import BlocksInMatrix
>>> from pybimstab.slipsurface import TortuousSurface
>>> slope = AnthropicSlope(slopeHeight=12, slopeDip=[1, 1.5],
>>> crownDist=10, toeDist=10)
>>> bim = BlocksInMatrix(slopeCoords=slope.coords, blockProp=0.25,
>>> tileSize=0.25, seed=123)
>>> surface = TortuousSurface(bim, dist1=0, dist2=17)
>>> surface.defineStructre()
array([[ 0. , 0.125 , 0.375 , 0.625 , ... ],
[16.57142857, 16.44642857, 16.19642857, 15.94642857, ... ]])
'''
import numpy as np
from pybimstab.astar import Astar, PreferredPath
from pybimstab.smoothcurve import SmoothCurve
# Obtaining the optimum path of the slip surface through A* algorithm
if self.preferredPath is not None:
preferredPathIdx = np.array(
[self.getIndexes(cr) for cr in self.preferredPath.T]).T
prefPathAstar = PreferredPath(preferredPathIdx, self.prefPathFact)
else:
prefPathAstar = None
astar = Astar(
grid=self.bim.grid, startNode=self.startIdx,
goalNode=self.goalIdx, heuristic=self.heuristic,
reverseLeft=self.reverseLeft, reverseUp=self.reverseUp,
preferredPath=prefPathAstar)
# Transform the optimum path indexes to real-scale coordinates
coords = [self.getCoord(tuple(idx)) for idx in astar.optimumPath.T]
# Appendinding the ends to the path
coords.append(self.end1)
coords.insert(0, self.end2)
# Sorting the path such that the begining is in the left side
coords = np.fliplr(np.array(coords).T)
if self.smoothFactor > 0:
smoothedCoords = SmoothCurve(x=coords[0], y=coords[1],
k=self.smoothFactor, n=500)
coords = smoothedCoords.smoothing
setattr(self, 'coords', coords)
return coords
[docs] def plot(self):
'''Method for generating a graphic of the tortuous slip surface and the
slope.
Returns:
(`matplotlib.figure.Figure`): object with the matplotlib structure\
of the plot. You might use it to save the figure for example.
Examples:
>>> from pybimstab.slope import AnthropicSlope
>>> from pybimstab.bim import BlocksInMatrix
>>> from pybimstab.slipsurface import TortuousSurface
>>> slope = AnthropicSlope(slopeHeight=12, slopeDip=[1, 1.5],
>>> crownDist=10, toeDist=10)
>>> bim = BlocksInMatrix(slopeCoords=slope.coords, blockProp=0.25,
>>> tileSize=0.25, seed=123)
>>> # Not allowing to turn left and up
>>> surface = TortuousSurface(
>>> bim, dist1=0, dist2=17, heuristic='manhattan',
>>> reverseLeft=False, reverseUp=False, smoothFactor=0,
>>> preferredPath=None, prefPathFact=None)
>>> fig = surface.plot()
.. figure:: https://rawgit.com/eamontoyaa/pybimstab/master/examples/figures/slipsurface_TortuousSurface_example1.svg
:alt: slipsurface_TortuousSurface_example1
.. only:: html
:download:`example script<../examples/figuresScripts/slipsurface_TortuousSurface_example1.py>`.
>>> # Allowing to turn left and up (manhattan heusitic function)
>>> surface = TortuousSurface(
>>> bim, dist1=0, dist2=17, heuristic='manhattan',
>>> reverseLeft=True, reverseUp=True, smoothFactor=0,
>>> preferredPath=None, prefPathFact=None)
>>> fig = surface.plot()
.. figure:: https://rawgit.com/eamontoyaa/pybimstab/master/examples/figures/slipsurface_TortuousSurface_example2.svg
:alt: slipsurface_TortuousSurface_example2
.. only:: html
:download:`example script<../examples/figuresScripts/slipsurface_TortuousSurface_example2.py>`.
>>> from numpy import array
>>> from pybimstab.slope import NaturalSlope
>>> from pybimstab.bim import BlocksInMatrix
>>> from pybimstab.slipsurface import CircularSurface, TortuousSurface
>>> terrainCoords = array(
>>> [[-2.49, 0.1, 1.7, 3.89, 5.9, 8.12, 9.87, 13.29, 20.29,
>>> 21.43, 22.28, 23.48, 24.65, 25.17],
>>> [18.16, 17.88, 17.28, 15.73, 14.31, 13.58, 13, 3.61, 3.61,
>>> 3.32, 2.71, 2.23, 1.21, 0.25]])
>>> slope = NaturalSlope(terrainCoords)
>>> bim = BlocksInMatrix(slopeCoords=slope.coords, blockProp=0.3,
>>> tileSize=0.35, seed=123)
>>> preferredPath = CircularSurface(
>>> slopeCoords=slope.coords, dist1=5, dist2=15.78, radius=20)
>>> # With a preferred path and smoothing the surface
>>> surface = TortuousSurface(
>>> bim, dist1=4, dist2=15.78, heuristic='euclidean',
>>> reverseLeft=False, reverseUp=False, smoothFactor=2,
>>> preferredPath=preferredPath.coords, prefPathFact=2)
>>> fig = surface.plot()
.. figure:: https://rawgit.com/eamontoyaa/pybimstab/master/examples/figures/slipsurface_TortuousSurface_example3.svg
:alt: slipsurface_TortuousSurface_example3
.. only:: html
:download:`example script<../examples/figuresScripts/slipsurface_TortuousSurface_example3.py>`.
>>> # Without a preferred path and smoothing the surface
>>> surface = TortuousSurface(
>>> bim, dist1=5, dist2=15.78, heuristic='euclidean',
>>> reverseLeft=False, reverseUp=False, smoothFactor=2,
>>> preferredPath=None)
>>> fig = surface.plot()
.. figure:: https://rawgit.com/eamontoyaa/pybimstab/master/examples/figures/slipsurface_TortuousSurface_example4.svg
:alt: slipsurface_TortuousSurface_example4
.. only:: html
:download:`example script<../examples/figuresScripts/slipsurface_TortuousSurface_example4.py>`.
'''
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import LinearSegmentedColormap as newcmap
# Variables to control the color map and its legend
if np.any(self.bim.grid == -1):
cmap = newcmap.from_list('BIMcmap',
['white', 'lightgray', 'black'], 3)
ticks = [-1+0.333, 0, 1-0.333]
ticksLabels = ['None', 'Matrix', 'Blocks']
else:
cmap = newcmap.from_list('BIMcmap', ['lightgray', 'black'], 2)
ticks = [0.25, 0.75]
ticksLabels = ['Matrix', 'Blocks']
# Plot body
fig = plt.figure()
ax = fig.add_subplot(111)
bar = ax.pcolormesh(self.bim.xCells, self.bim.yCells, self.bim.grid,
cmap=cmap)
if self.preferredPath is not None:
ax.plot(self.preferredPath[0], self.preferredPath[1], ':r',
lw=2.0, label='Preferred\npath')
ax.plot(self.coords[0], self.coords[1], '-r', lw=2.5,
label='Tortuous\nsurface')
ax.plot(self.bim.slopeCoords[0], self.bim.slopeCoords[1], '-k')
# Configuring the colorbar
bar = plt.colorbar(bar, ax=ax, ticks=ticks, pad=0.03,
shrink=0.15, aspect=3)
bar.ax.set_yticklabels(ticksLabels, fontsize='small')
# Plot settings
ax.set_aspect(1)
ax.legend(fontsize='small', bbox_to_anchor=(1.005, 1), loc='best')
ax.grid(True, ls='--', lw=0.5)
ax.set_xlim((-0.02*self.bim.slopeCoords[0].max(),
1.02*self.bim.slopeCoords[0].max()))
ax.set_ylim((-0.02*self.bim.slopeCoords[1].max(),
1.02*self.bim.slopeCoords[1].max()))
fig.tight_layout()
return fig
# %%
'''
BSD 2 license.
Copyright (c) 2018, Universidad Nacional de Colombia, Exneyder Andres Montoya
Araque and Ludger O. Suarez-Burgoa.
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
'''