import warnings
import numpy as np
from apsg.helpers._helper import is_jsonable
from apsg.helpers._math import atan2d
from apsg.helpers._notation import (
fol2vec_rhr,
format_linear,
format_planar,
geo2vec_planar,
parse_quadrant_linear,
parse_quadrant_planar,
vec2geo_linear,
vec2geo_linear_signed,
vec2geo_planar,
vec2geo_planar_signed,
)
from apsg.math._vector import Axial2, Axial3, Vector3
[docs]
class Direction(Axial2):
"""
A class to represent axial (non-oriented) 2D linear feature (direction).
There are different way to create ``Direction`` object:
- without arguments create default ``Direction`` D:0
- with single argument `d`, where:
- `d` could be Vector2-like object
- `d` could be string 'x' or 'y' - principal axes of coordinate system
- `d` could be tuple of (x, y) - vector components
- with 1 argument direction
- with 2 numerical arguments defining vector components
Args:
direction (float): plunge direction of linear feature in degrees.
Examples:
>>> dir2()
>>> dir2('y')
>>> dir2(45)
>>> d = dir2(1, -1)
"""
def __repr__(self):
return f"D:{self.direction:.0f}"
[docs]
class Lineation(Axial3):
"""
A class to represent axial (non-oriented) linear feature (lineation).
There are different way to create ``Lineation`` object:
- without arguments create default ``Lineation`` L:0/0
- with single argument `l`, where:
- `l` could be Vector3-like object
- `l` could be string 'x', 'y' or 'z' - principal axes of coordinate system
- `l` could be quadrant notation string, e.g. 'N45E,30'
- `l` could be tuple of (x, y, z) - vector components
- with 2 arguments plunge direction and plunge
- with 3 numerical arguments defining vector components
Args:
azi (float): plunge direction of linear feature in degrees.
inc (float): plunge of linear feature in degrees.
Examples:
>>> lin()
>>> lin('y')
>>> lin(1,2,-1)
>>> l = lin(110, 26)
>>> l = lin('N45E,30')
"""
def __init__(self, *args, **kwargs):
if len(args) == 1 and isinstance(args[0], str):
try:
args = parse_quadrant_linear(args[0])
except ValueError:
pass # not quadrant text; let Vector3.__init__ handle 'x'/'y'/'z' or raise
super().__init__(*args, **kwargs)
def __repr__(self):
return f"L:{format_linear(*self.geo)}"
[docs]
def cross(self, other):
"""Return Foliation defined by two linear features."""
return Foliation(super().cross(other))
__pow__ = cross
@property
def geo(self):
"""Return tuple of plunge direction and plunge."""
return vec2geo_linear(self)
[docs]
def to_json(self):
"""Return as JSON dict."""
azi, inc = vec2geo_linear_signed(self)
return {
"datatype": type(self).__name__,
"args": (azi, inc),
"kwargs": self._attrs,
}
[docs]
class Foliation(Axial3):
"""
A class to represent non-oriented planar feature (foliation).
There are different way to create ``Foliation`` object:
- without arguments create default ``Foliation`` S:180/0
- with single argument `f`, where:
- `f` could be Vector3-like object
- `f` could be string 'x', 'y' or 'z' - principal planes of coordinate system
- `f` could be quadrant notation string, e.g. 'N30E,40NW'
- `f` could be tuple of (x, y, z) - vector components
- with 2 arguments follows active notation. See apsg_conf.notation
- with 3 numerical arguments defining vector components of plane normal
Args:
azi (float): dip direction (or strike) of planar feature in degrees.
inc (float): dip of planar feature in degrees.
Examples:
>>> fol()
>>> fol('y')
>>> fol(1,2,-1)
>>> f = fol(250, 30)
>>> f = fol('N30E,40NW')
"""
def __init__(self, *args, **kwargs):
if len(args) == 0:
coords = (0, 0, 1)
elif len(args) == 1:
if np.asarray(args[0]).shape == Foliation.__shape__:
coords = np.asarray(args[0])
elif isinstance(args[0], str):
if args[0].lower() == "x":
coords = (1, 0, 0)
elif args[0].lower() == "y":
coords = (0, 1, 0)
elif args[0].lower() == "z":
coords = (0, 0, 1)
else:
try:
strike, dip = parse_quadrant_planar(args[0])
except ValueError:
raise TypeError(
f"Not valid arguments for {type(self).__name__}"
) from None
coords = fol2vec_rhr(strike, dip)
else:
raise TypeError(f"Not valid arguments for {type(self).__name__}")
elif len(args) == 2:
coords = geo2vec_planar(*args)
elif len(args) == 3:
coords = [float(v) for v in args]
else:
raise TypeError("Not valid arguments for Foliation")
self._coords = tuple(coords)
if is_jsonable(kwargs):
self._attrs = kwargs
else:
raise TypeError("Provided attributes are not serializable.")
def __repr__(self):
return f"S:{format_planar(*self.geo)}"
[docs]
def cross(self, other):
"""Return Lineation defined by intersection of planar features."""
return Lineation(super().cross(other))
__pow__ = cross
@property
def geo(self):
"""Return tuple of dip direction and dip."""
return vec2geo_planar(self)
[docs]
def to_json(self):
"""Return as JSON dict."""
azi, inc = vec2geo_planar_signed(self)
return {
"datatype": type(self).__name__,
"args": (azi, inc),
"kwargs": self._attrs,
}
[docs]
def dipvec(self):
"""Return dip vector."""
return Vector3(*vec2geo_planar(self))
[docs]
def strike(self):
"""Return strike as Direction."""
n = self.uv()
return Direction((atan2d(n.y, n.x) + 90) % 360)
[docs]
def pole(self):
"""Return plane normal as vector."""
return Vector3(self)
[docs]
def rake(self, rake):
"""Return rake vector."""
return Vector3(self.dipvec().rotate(self, rake - 90))
[docs]
class Pair:
"""
The class to store pair of planar and linear feature.
When ``Pair`` object is created, both planar and linear feature are
adjusted, so linear feature perfectly fit onto planar one. Warning
is issued, when misfit angle is bigger than 20 degrees.
There are different way to create ``Pair`` object:
- without arguments create default Pair with fol(0,0) and lin(0,0)
- with single argument `p`, where:
- `p` could be Pair
- `p` could be tuple of (fazi, finc, lazi, linc)
- `p` could be tuple of (fx, fy ,fz, lx, ly, lz)
- with 2 arguments f and l could be Vector3 like objects,
e.g. Foliation and Lineation
- with four numerical arguments defining `fol(fazi, finc)` and `lin(lazi, linc)`
Args:
fazi (float): dip azimuth of planar feature in degrees.
finc (float): dip of planar feature in degrees.
lazi (float): plunge direction of linear feature in degrees.
linc (float): plunge of linear feature in degrees.
Attributes:
fvec (Vector3): corrected vector normal to plane.
lvec (Vector3): corrected vector of linear feature.
Examples:
>>> pair()
>>> pair(p)
>>> pair(f, l)
>>> pair(fazi, finc, lazi, linc)
>>> p = pair(140, 30, 110, 26)
"""
__slots__ = ("fvec", "lvec", "misfit", "_attrs")
__shape__ = (6,)
def __init__(self, *args, **kwargs):
if len(args) == 0:
fvec, lvec = Vector3(0, 0, 1), Vector3(1, 0, 0)
elif len(args) == 1 and isinstance(args[0], Pair):
fvec, lvec = args[0].fvec, args[0].lvec
elif len(args) == 1 and np.asarray(args[0]).shape == (4,):
fazi, finc, lazi, linc = (float(v) for v in args[0])
fvec, lvec = Foliation(fazi, finc), Lineation(lazi, linc)
elif len(args) == 1 and np.asarray(args[0]).shape == Pair.__shape__:
fvec, lvec = Vector3(args[0][:3]), Vector3(args[0][-3:])
elif len(args) == 2:
if isinstance(args[0], Vector3) and isinstance(args[1], Vector3):
fvec, lvec = args
else:
raise TypeError("Not valid arguments for Pair")
elif len(args) == 4:
fvec = Foliation(args[0], args[1])
lvec = Lineation(args[2], args[3])
else:
raise TypeError("Not valid arguments for Pair")
fvec = Vector3(fvec)
lvec = Vector3(lvec)
misfit = abs(90 - fvec.angle(lvec))
if misfit > 20:
warnings.warn(f"Warning: Misfit angle is {misfit:.1f} degrees.")
ax = fvec.cross(lvec)
ang = (lvec.angle(fvec) - 90) / 2
self.fvec = Vector3(fvec.rotate(ax, ang))
self.lvec = Vector3(lvec.rotate(ax, -ang))
self.misfit = misfit
if is_jsonable(kwargs):
self._attrs = kwargs
else:
raise TypeError("Provided attributes are not serializable.")
def __repr__(self):
return f"P:{format_planar(*self.fol.geo)}-{format_linear(*self.lin.geo)}"
def __eq__(self, other):
"""Return `True` if pairs are equal, otherwise `False`."""
cls = type(self)
if not isinstance(other, cls):
if np.asarray(other).shape == cls.__shape__:
other = cls(other)
else:
return NotImplemented
return (self.fvec == other.fvec) and (self.lvec == other.lvec)
def __ne__(self, other):
"""Return `True` if pairs are not equal, otherwise `False`."""
return not self == other
def __array__(self, dtype=None, copy=None):
return np.hstack((self.fvec, self.lvec)).astype(dtype)
[docs]
def label(self):
"""Return label."""
return str(self)
[docs]
def to_json(self):
"""Return as JSON dict."""
fazi, finc = vec2geo_planar_signed(self.fvec)
lazi, linc = vec2geo_linear_signed(self.lvec)
return {
"datatype": type(self).__name__,
"args": (fazi, finc, lazi, linc),
"kwargs": self._attrs,
}
[docs]
@classmethod
def random(cls):
"""Random Pair."""
lin, p = Vector3.random(), Vector3.random()
fol = lin.cross(p)
return cls(fol, lin)
[docs]
def rotate(self, axis, phi):
"""Rotates ``Pair`` by angle `phi` about `axis`.
Args:
axis (``Vector3``): axis of rotation
phi (float): angle of rotation in degrees
Examples:
>>> p = pair(fol(140, 30), lin(110, 26))
>>> p.rotate(lin(40, 50), 120)
P:210/83-287/60
Returns:
Pair: The rotated pair.
"""
try:
axis = Vector3(axis)
except TypeError:
raise TypeError("Unsupported argument for rotate. Expecting Vector3")
return type(self)(self.fvec.rotate(axis, phi), self.lvec.rotate(axis, phi))
@property
def rax(self):
return self.fvec.cross(self.lvec)
@property
def rake(self):
"""Return a rake of linear feature on planar feature of ``Pair`` in degrees."""
return -np.degrees(
np.atan2(-self.lvec.dot(self.fol.rake(90)), self.lvec.dot(self.fol.rake(0)))
)
@property
def fol(self):
"""Return a planar feature of ``Pair`` as ``Foliation``."""
return Foliation(self.fvec)
@property
def lin(self):
"""Return a linear feature of ``Pair`` as ``Lineation``."""
return Lineation(self.lvec)
[docs]
class Fault(Pair):
"""
The class to store ``Pair`` with associated sense of movement.
When ``Fault`` object is created, both planar and linear feature are
adjusted, so linear feature perfectly fit onto planar one. Warning
is issued, when misfit angle is bigger than 20 degrees.
There are different way to create ``Fault`` object:
- without arguments create default ``Fault`` with `fol(0,0)` and `lin(0,0)`
- with single argument `p`:
- `p` could be Fault
- `p` could be tuple of (fazi, finc, lazi, linc, sense)
- `p` could be tuple of (fx, fy ,fz, lx, ly, lz)
- with 2 arguments p (Pair object) and sense
- with 3 arguments f, l (Vector3 like objects), e.g. ``Foliation``
and ``Lineation`` and sense
- with 3 arguments fazi, finc, rake
- with 5 arguments fazi, finc, lazi, linc, sense
Args:
fazi (float): dip azimuth of planar feature in degrees.
finc (float): dip of planar feature in degrees.
lazi (float): plunge direction of linear feature in degrees.
linc (float): plunge of linear feature in degrees.
sense (float or str): sense of movement +/-1 hanging-wall down/up. When str,
must be one of 's', 'd', 'n', 'r'.
Attributes:
fvec (Vector3): corrected vector normal to plane.
lvec (Vector3): corrected vector of linear feature.
sense (int): sense of movement (+/-1).
Examples:
>>> f = fault(140, 30, 110, 26, -1)
>>> f = fault(140, 30, 110, 26, 'r')
>>> p = pair(140, 30, 110, 26)
>>> f = fault(p, 'n')
>>> f = fault(fol(120, 80), lin(32, 10), 's')
"""
__shape__ = (7,)
def __init__(self, *args, **kwargs):
if len(args) == 0:
fvec, lvec = Vector3(0, 0, 1), Vector3(1, 0, 0)
elif len(args) == 1 and np.asarray(args[0]).shape == (5,):
fazi, finc, lazi, linc, sense = (float(v) for v in args[0])
fvec, lvec = Foliation(fazi, finc), Lineation(lazi, linc)
sense = self.calc_sense(fvec, lvec, sense)
if sense < 0:
lvec = -lvec
elif len(args) == 1 and isinstance(args[0], Pair):
fvec, lvec = args[0].fvec, args[0].lvec
elif len(args) == 2 and isinstance(args[0], Pair):
fvec, lvec = args[0].fvec, args[0].lvec
if Fault(fvec, lvec).sense != self.calc_sense(fvec, lvec, args[1]):
lvec = -lvec
elif len(args) == 2:
if isinstance(args[0], Vector3) and isinstance(args[1], Vector3):
fvec, lvec = args[0], args[1]
elif len(args) == 3:
if isinstance(args[0], Vector3) and isinstance(args[1], Vector3):
fvec, lvec = args[0], args[1]
if Fault(fvec, lvec).sense != self.calc_sense(fvec, lvec, args[2]):
lvec = -lvec
else:
fvec = Foliation(args[0], args[1])
lvec = fvec.rake(args[2])
elif len(args) == 5:
fvec = Foliation(args[0], args[1])
lvec = Lineation(args[2], args[3])
if Fault(fvec, lvec).sense != self.calc_sense(fvec, lvec, args[4]):
lvec = -lvec
else:
raise TypeError("Not valid arguments for Fault")
super().__init__(fvec, lvec)
if is_jsonable(kwargs):
self._attrs = kwargs
else:
raise TypeError("Provided attributes are not serializable.")
@classmethod
def calc_sense(cls, fvec, lvec, sense) -> int:
if isinstance(sense, str):
p = Pair(fvec, lvec)
if sense.lower() == "s":
if abs(p.rake) < 90:
res = 1
else:
res = -1
elif sense.lower() == "d":
if abs(p.rake) < 90:
res = -1
else:
res = 1
elif sense.lower() == "n":
res = 1
elif sense.lower() == "r":
res = -1
else:
raise ValueError(
f"Invallid sense '{sense}'. Must be one of 's', 'd', 'n', 'r'."
)
return res
else:
return int(sense)
def __repr__(self):
schar = self.sense_str
return (
f"F:{format_planar(*self.fol.geo)}-{format_linear(*self.lin.geo)} {schar}"
)
def __eq__(self, other):
"""Return `True` if pairs are equal, otherwise `False`."""
cls = type(self)
if not isinstance(other, cls):
if np.asarray(other).shape == cls.__shape__:
other = cls(other)
else:
return NotImplemented
return (
(self.fvec == other.fvec)
and (self.lvec == other.lvec)
and (self.sense == other.sense)
)
def __ne__(self, other):
"""Return `True` if pairs are not equal, otherwise `False`."""
return not self == other
def __array__(self, dtype=None, copy=None):
return np.hstack((self.fvec, self.lvec, self.sense)).astype(dtype)
[docs]
def to_json(self):
"""Return as JSON dict."""
fazi, finc = vec2geo_planar_signed(self.fvec)
lazi, linc = vec2geo_linear_signed(self.lvec)
return {
"datatype": type(self).__name__,
"args": (fazi, finc, lazi, linc, self.sense),
"kwargs": self._attrs,
}
[docs]
@classmethod
def random(cls):
"""Random Fault."""
import random
lvec, p = Vector3.random(), Vector3.random()
fvec = lvec.cross(p)
return cls(fvec, lvec, random.choice([-1, 1]))
@property
def sense(self):
return int(np.sign(self.rake))
@property
def sense_str(self):
if abs(self.rax.geo[1]) > 60:
if abs(self.rake) < 90:
schar = "S"
else:
schar = "D"
else:
if self.rake > 0:
schar = "N"
else:
schar = "R"
return schar
[docs]
def p_vector(self, ptangle=90):
"""Return P axis as ``Vector3``."""
return self.fvec.rotate(self.lvec.cross(self.fvec), -ptangle / 2)
[docs]
def t_vector(self, ptangle=90):
"""Return T-axis as ``Vector3``."""
return self.fvec.rotate(self.lvec.cross(self.fvec), ptangle / 2)
@property
def p(self):
"""Return P-axis as ``Lineation``."""
return Lineation(self.p_vector())
@property
def t(self):
"""Return T-axis as ``Lineation``."""
return Lineation(self.t_vector())
@property
def m(self):
"""Return kinematic M-plane as ``Foliation``."""
return Foliation(self.lvec.cross(self.fvec))
@property
def d(self):
"""Return dihedra plane as ``Fol``."""
return Foliation(self.lvec.cross(self.fvec).cross(self.fvec))
[docs]
def angular_misfit(self, sigma):
"""Angular misfit (°) between observed slip and predicted shear-traction direction.
Args:
sigma (Stress3): Stress tensor
Returns:
float: The angular misfit in degrees.
"""
return self.lvec.angle(sigma.fault(self.fvec).lvec)
[docs]
class Cone:
"""
The class to store cone with given axis, secant line and revolution angle
in degrees.
There are different way to create ``Cone`` object according to number
of arguments:
- without args, you can create default``Cone`` with axis ``lin(0, 90)``,
secant ``lin(0, 0)`` angle 360°
- with single argument `c`, where `c` could be ``Cone``, 5-tuple of
`(aazi, ainc, sazi, sinc, revangle)` or 7-tuple of
`(ax, ay ,az, sx, sy, sz, revangle)`
- with 3 arguments, where axis and secant line could be Vector3 like objects,
e.g. Lineation and third argument is revolution angle
- with 5 arguments defining axis `lin(aazi, ainc)`, secant line
`lin(sazi, sinc)` and angle of revolution
Args:
*args: Variable length argument list. See descriptions above.
**kwargs: Additional keyword arguments.
Attributes:
axis (Vector3): axis of the cone
secant (Vector3): secant line
revangle (float): revolution angle
Examples:
>>> cone()
>>> cone(c)
>>> cone(a, s, revangle)
>>> cone(aazi, ainc, sazi, sinc, revangle)
>>> c = cone(140, 30, 110, 26, 360)
"""
__slots__ = ("axis", "secant", "revangle", "_attrs")
__shape__ = (7,)
def __init__(self, *args, **kwargs):
if len(args) == 0:
axis, secant, revangle = Vector3(0, 0, 1), Vector3(1, 0, 0), 360
elif len(args) == 1 and isinstance(args[0], Cone):
axis, secant, revangle = args[0].axis, args[0].secant, args[0].revangle
elif len(args) == 1 and np.asarray(args[0]).shape == (5,):
aazi, ainc, sazi, sinc, revangle = (float(v) for v in args[0])
axis, secant = Lineation(aazi, ainc), Lineation(sazi, sinc)
elif len(args) == 1 and np.asarray(args[0]).shape == Cone.__shape__:
axis, secant, revangle = (
Vector3(args[0][:3]),
Vector3(args[0][3:6]),
args[0][-1],
)
elif len(args) == 2:
if isinstance(args[0], Vector3) and isinstance(args[1], Vector3):
axis, secant = args
revangle = 360
elif isinstance(args[0], Vector3) and np.isscalar(args[1]):
axis = args[0]
azi, inc = axis.geo
secant = Vector3(azi, inc + args[1])
revangle = 360
else:
raise TypeError("Not valid arguments for Cone")
elif len(args) == 3:
if isinstance(args[0], Vector3) and isinstance(args[1], Vector3):
axis, secant, revangle = args
else:
raise TypeError("Not valid arguments for Cone")
elif len(args) == 4:
axis = Lineation(args[0], args[1])
secant = Lineation(args[2], args[3])
revangle = 360
elif len(args) == 5:
axis = Lineation(args[0], args[1])
secant = Lineation(args[2], args[3])
revangle = args[4]
else:
raise TypeError("Not valid arguments for Cone")
self.axis = Vector3(axis)
self.secant = Vector3(secant)
self.revangle = float(revangle)
if self.axis.angle(self.secant) > 90:
self.secant = -self.secant
if is_jsonable(kwargs):
self._attrs = kwargs
else:
raise TypeError("Provided attributes are not serializable.")
def __repr__(self):
azi, inc = vec2geo_linear(self.axis)
return f"C:{format_linear(azi, inc)} [{self.apical_angle():g}]"
def __eq__(self, other):
cls = type(self)
if not isinstance(other, cls):
if np.asarray(other).shape == cls.__shape__:
other = cls(other)
else:
return NotImplemented
return (
(self.axis == other.axis)
and (self.secant == other.secant)
and (self.revangle == other.revangle)
)
def __ne__(self, other):
return not self == other
def __array__(self, dtype=None, copy=None):
return np.hstack((self.axis, self.secant, self.revangle)).astype(dtype)
[docs]
def label(self):
"""Return label."""
return str(self)
[docs]
def to_json(self):
"""Return as JSON dict."""
aazi, ainc = vec2geo_linear_signed(self.axis)
sazi, sinc = vec2geo_linear_signed(self.secant)
return {
"datatype": type(self).__name__,
"args": (aazi, ainc, sazi, sinc, self.revangle),
"kwargs": self._attrs,
}
[docs]
@classmethod
def random(cls):
"""Random Cone."""
axis, secant = Vector3.random(), Vector3.random()
return cls(axis, secant, 360)
[docs]
def rotate(self, axis, phi):
"""Rotates ``Cone`` by angle `phi` about `axis`.
Args:
axis (``Vector3``): axis of rotation
phi (float): angle of rotation in degrees
Examples:
>>> c = cone(lin(140, 30), lin(110, 26), 360)
>>> c.rotate(lin(40, 50), 120)
C:210/83-287/60
Returns:
Cone: The rotated cone.
"""
try:
axis = Vector3(axis)
except TypeError:
raise TypeError("Unsupported argument for rotate. Expecting Vector3")
return type(self)(
self.axis.rotate(axis, phi), self.secant.rotate(axis, phi), self.revangle
)
[docs]
def apical_angle(self):
"""Return apical angle."""
return self.axis.angle(self.secant)
@property
def rotated_secant(self):
"""Return revangle rotated secant vector."""
return self.secant.rotate(self.axis, self.revangle)