# Tutorial¶

APSG defines several new python classes to easily manage, analyze and visualize orientation structural geology data. Base class `Vec3` is derived from `numpy.array` class and offers several new method which will be explained in following examples.

## Basic usage¶

APSG module could be imported either into own name space or into active one for easier interactive work:

```>>> from apsg import *
```

## Basic operations with vectors¶

Instance of vector object `Vec3` could be created from any iterable object as list, tuple or array:

```>>> u = Vec3([1, -2, 3])
>>> v = Vec3([-2, 1, 1])
```

For common vector operation we can use standard mathematical operators or special methods using dot notation:

```>>> u + v
V(-1.000, -1.000, 4.000)
>>> u - v
V(3.000, -3.000, 2.000)
>>> 3*u - 2*v
V(7.000, -8.000, 7.000)
```

Its magnitude or length is most commonly defined as its Euclidean norm and could be calculated using `abs`:

```>>> abs(v)
2.4494897427831779
>>> abs(u+v)
4.2426406871192848
```

For dot product we can use multiplication operator `*` or `dot` method:

```>>> u*v
-1
>>> u.dot(v)
-1
```

For cross product we can use operator `**` or method `cross`:

```>>> u**v
V(-5.000, -7.000, -3.000)
>>> u.cross(v)
V(-5.000, -7.000, -3.000)
```

To project vector `u` onto vector `v` we can use method `proj`:

```>>> u.proj(v)
V(0.333, -0.167, -0.167)
```

To find angle (in degrees) between to vectors we use method `angle`:

```>>> u.angle(v)
96.263952719927218
```

Method `rotate` provide possibility to rotate vector around another vector. For example, to rotate vector `u` around vector `v` for 45°:

```>>> u.rotate(v, 45)
V(2.248, 0.558, 2.939)
```

## Classes Lin and Fol¶

To work with orientation data in structural geology, APSG provide two classes derived from `Vec3` class. There is `Fol` class to represent planar features by planes and `Lin` class to represent linear feature by lines. Both classes provide all `Vec3` methods, but they differ in way how instance is created and how some operations are calculated, as structural geology data are commonly axial in nature.

To create instance of `Lin` or `Fol` class, we have to provide dip direction and dip, both in degrees:

```>>> Lin(120, 60)
L:120/60
>>> Fol(216, 62)
S:216/62
```

or we can create instance from `Vec3` object using `aslin` and `asfol` properties:

```>>> u.aslin
L:297/53
>>> u.asfol
S:117/37
```

## Vec3 methods for Lin and Fol¶

To find angle between two linear or planar features we can use method `angle`:

```>>> l1 = Lin(110, 40)
>>> l2 = Lin(160, 30)
>>> l1.angle(l2)
41.597412680035468
>>> p1 = Fol(330, 50)
>>> p2 = Fol(250, 40)
>>> p1.angle(p2)
54.696399321975328
```

We can use cross product to construct planar feature defined by two linear features:

```>>> l1**l2
S:113/40
```

or to construct linear feature defined as intersection of two planar features:

```>>> p1**p2
L:278/36
```

Cross product of planar and linear features could be used to construct plane defined by linear feature and normal of planar feature:

```>>> l2**p2
S:96/53
```

or to find perpendicular linear feature on given plane:

```>>> p2**l2
L:276/37
```

To rotate structural features we can use method `rotate`:

```>>> p2.rotate(l2, 45)
S:269/78
```

## Classes Pair and Fault¶

To work with paired orientation data like foliations and lineations or fault data in structural geology, APSG provide two base `Pair` class and derived `Fault` class. Both classes are instantiated providing dip direction and dip of planar and linear measurements, which are automatically orthogonalized. If misfit is too high, warning is raised. The `Fault` class expects one more argument providing sense of movement information, either 1 or -1.

To create instance of `Pair` class, we have to provide dip direction and dip of planar and linear feature, both in degrees:

```>>> p = Pair(120, 40, 162, 28)
>>> p
P:118/39-163/30
>>> p.misfit
3.5623168411508175
>>> type(p)
<class 'apsg.core.Pair'>
```

Planar and linear features are accessible using `fol` and `lin` properties:

```>>> p.fol
S:118/39
>>> p.lin
L:163/30
>>> type(p.fol)
<class 'apsg.core.Fol'>
>>> type(p.lin)
<class 'apsg.core.Lin'>
```

To rotate `Pair` instance we can use `rotate` method:

```>>> p.rotate(Lin(45, 10), 60)
P:314/83-237/61
```

Instantiation of `Fault` class is similar, we just have to provide argument to define sense of movement:

```>>> f = Fault(120, 60, 110, 58, -1)  # -1 for normal fault
>>> f
F:120/59-110/59 -
```

Note the change in sense of movement after `Fault` rotation:

```>>> f.rotate(Lin(45, 10), 60)
F:312/62-340/59 +
```

For simple fault analyses `Fault` class also provide `p`, `t`, `m` and `d` properties to get PT-axes, kinematic plane and dihedra separation plane:

```>>> f.p
L:315/75
>>> f.t
L:116/14
>>> f.m
S:27/85
>>> f.d
S:290/31
```

## Group class¶

`Group` class serve as a homogeneous container for `Lin`, `Fol` and `Vec3` objects. It allows grouping of features either for visualization or batch analysis:

```>>> g = Group([Lin(120,60), Lin(116,50), Lin(132,45), Lin(90,60), Lin(84,52)],
>>>           name='L1')
>>> g
L1: 5 Lin
```

To simplify interactive group creation, you can use function `G`:

```>>> g = G('120 60 116 50 132 45 90 60 84 52', name='L1')
```

Method `len` returns number of features in group:

```>>> len(g)
5
```

Most of the `Lin`, `Fol` and `Vec3` methods could be used for `Group` as well. For example, to measure angles between all features in group and another feature, we can use method `angle`:

```>>> g.angle(Lin(110,50))
array([ 11.49989817,   3.85569115,  15.61367789,  15.11039885,  16.3947936 ])
```

To rotate all features in group around another feature, we can use method `rotate`:

```>>> gr = g.rotate(Lin(150, 30), 45)
```

To show data in list you can use `data` property:

```>>> gr.data
[L:107/35, L:113/26, L:126/30, L:93/26, L:94/18]
```

Property `R` gives mean or resultant of all features in group. Note that `Lin` and `Fol` are axial in nature, so resultant vector is not reliable. You can use `ortensor` property.:

```>>> g.R
L:110/55
```

`Group` class offers several methods to infer spherical statistics as spherical variance, Fisher’s statistics, confidence cones on data etc.:

```>>> g.var
0.02337168447438509
>>> g.fisher_stats
{'k': 34.229454059110871, 'a95': 13.264029905117329, 'csd': 13.844747281750971}
>>> g.delta
12.411724720740544
```

To calculate orientation tensor of all features in group, we can use `ortensor` property.:

```>>> g.ortensor
Ortensor: L1 Kind: LLS
(E1:0.954,E2:0.04021,E3:0.005749)
[[ 0.07398181 -0.09605477 -0.14324311]
[-0.09605477  0.28446118  0.42092899]
[-0.14324311  0.42092899  0.64155701]]
```

## Ortensor class¶

`Ortensor` class represents orientation tensor of set of planar or linear features. Eigenvalues and eigenvectors could be obtained by methods `eigenvals` and `eigenvects`. Eigenvectors could be also represented by linear or planar features using properties `eigenlins` and `eigenfols`. Several properties to describe orientation distribution is also impleneted, e.g. Woodcock’s `shape` and `strength` or Vollmer’s `P`, `G`, `R` and `C` indexes.:

```>>> ot = Ortensor.from_group(g)
>>> ot.eigenvals
(0.954038468659639, 0.04021274946196462, 0.005748781878396576)
>>> ot.eigenvects.data
[V(0.192, -0.542, -0.818), V(-0.981, -0.082, -0.176), V(-0.028, -0.836, 0.547)]
>>> ot.eigenlins.data
[L:110/55, L:5/10, L:268/33]
>>> ot.eigenfols.data
[S:290/35, S:185/80, S:88/57]
>>> ot.strength, ot.shape
(5.111716009046873, 1.6278666609093613)
>>> ot.P, ot.G, ot.R
(0.91382571919767419, 0.068927935167136606, 0.017246345635188932)
```

## StereoNet class¶

When `StereoNet` class instance is created with arguments, they are immidiatelly plotted and shown. Most of the objects provided by APSG could be plotted. Here we use the `Group` object and principal planes and lines of `Ortensor`:

```>>> StereoNet(g, ot.eigenfols, ot.eigenlins)
``` When `StereoNet` class instance is created without arguments, several methods and properties could be used for additional operations. To finalize plot we can use `show` method (not needed in jupyter notebooks)

Any `Fol`, `Lin`, `Group` object could be visualized as plane, line or pole in stereographic projection using `StereoNet` class:

```>>> s = StereoNet()
>>> s.plane(Fol(150, 40))
>>> s.pole(Fol(150, 40))
>>> s.line(Lin(112, 30))
>>> s.show()
``` A cones (or small circles) could be plotted as well:

```>>> s = StereoNet()
>>> g = Group.randn_lin(mean=Lin(40, 15))
>>> s.line(g, 'k.')
>>> s.cone(g.R, g.fisher_stats['a95'], 'r')  # confidence cone on resultant
>>> s.cone(g.R, g.fisher_stats['csd'], 'g')  # confidence cone on 63% of data
>>> s.show()
``` To make density contours plots, a `contour` and `contourf` methods are available:

```>>> s = StereoNet()
>>> g = Group.randn_lin(mean=Lin(40, 20))
>>> s.contourf(g, 8, legend=True, sigma=2)
>>> s.line(g, 'g.')
>>> s.show()
``` Except `Group`, APSG provides `PairSet` and `FaultSet` classes to store `Pair` or `Fault` datasets. It can be inicialized by passing list of `Pair` or `Fault` objects as argument or use class methods `from_array` or `from_csv`:

```>>> p = PairSet([Pair(120, 30, 165, 20),
>>>              Pair(215, 60, 280,35),
>>>              Pair(324, 70, 35, 40)])
>>> p.misfit
array([ 2.0650076 ,  0.74600727,  0.83154705])
>>> StereoNet(p)
``` `StereoNet` has two special methods to visualize fault data. Method `fault` produce classical Angelier plot:

```>>> f = FaultSet([Fault(170, 60, 182, 59, -1),
>>>               Fault(210, 55, 195, 53, -1),
>>>               Fault(10, 60, 15, 59, -1),
>>>               Fault(355, 48, 22, 45, -1)])
>>> s = StereoNet()
>>> s.fault(f)
>>> s.line(f.p, label='P-axes')
>>> s.line(f.t, label='T-axes')
>>> s.plane(f.m, label='M-planes')
>>> s.show()
``` `hoeppner` method produce Hoeppner diagram and must be invoked from `StereoNet` instance:

```>>> s = StereoNet()
>>> s.hoeppner(f, label=repr(f))
>>> s.show()
``` Note that `fault` method is used, when data are passed directly to `StereoNet` instance:

```>>> f = Fault(120, 60, 110, 58, -1)
>>> StereoNet(f, f.m, f.d, f.p, f.t)
``` ## StereoGrid class¶

`StereoGrid` class allows to visualize any scalar field on StereoNet. Internally it is used for plotting contour diagrams, but it exposes `apply_func` method to calculate scalar field by any user-defined function. Function must accept three element `numpy.array` as first argument passed from grid points of `StereoGrid`.

Following example shows how to plot resolved shear stress on plane from given stress tensor. `StereoGrid.apply_func` method is used to calculate this value over all directions and finally values are plotted by `StereoNet`:

```>>> S = Stress([[-10, 2, -3],[2, -5, 1], [-3, 1, 2]])
>>> d = StereoGrid()
>>> d.apply_func(S.shear_stress)
>>> s = StereoNet()
>>> s.contourf(d, 10, legend=True)
>>> s.show()
``` The `FaultSet` provide also `amgmech` method which provide access to Angelier dihedra method. Results are stored in `StereoGrid`. Default behavior is to calculate counts (positive in extension, negative in compression):

```>>> f = FaultSet.examples('MELE')
>>> StereoNet(f.angmech())
``` Setting method to ‘probability’, maximum likelihood estimate is calculated.:

```>>> f = FaultSet.examples('MELE')
>>> StereoNet(f.angmech(method='probability'))
``` ## Fabric plots¶

APSG provides several fabric plots to visualize Tensor-type objects (`Ortensor`, `Ellipsoid`). `FlinnPlot` class provide classical Flinn’s deformation diagram,

```>>> g1 = Group.examples('B2')
>>> g2 = Group.examples('B4')
>>> g3 = Group.uniform_lin(name='Uniform')
>>> FlinnPlot(g1, g2, g3);
``` `RamsayPlot` class provide Ramsay modification of Flinn’s deformation diagram,

```>>> RamsayPlot(g1, g2, g3);
``` `VollmerPlot` class provide triangular fabric plot (Vollmer, 1989),

```>>> VollmerPlot(g1, g2, g3);
``` and `HsuPlot` class provide Hsu fabric diagram using natural strains.

```>>> HsuPlot(g1, g2, g3);
``` ## Cluster class¶

`Cluster` class provide access to scipy hierarchical clustering. Distance matrix is calculated as mutual angles of features within Group keeping axial and/or vectorial nature in mind. `Cluster.explain` method allows to explore explained variance versus number of clusters relation. Actual cluster is done by `Cluster.cluster` method, using distance or maxclust criterion. Using of `Cluster` is explained in following example. We generate some data and plot dendrogram:

```>>> g1 = Group.randn_lin(mean=Lin(45,30))
>>> g2 = Group.randn_lin(mean=Lin(320,56))
>>> g3 = Group.randn_lin(mean=Lin(150,40))
>>> g = g1 + g2 + g3
>>> cl = Cluster(g)
>>> cl.dendrogram(no_labels=True)
``` Now we can explore evolution of within-groups variance versus number of clusters on Elbow plot (Note change in slope for three clusters):

```>>> cl.elbow()
``` Finally we can do clustering and plot created clusters:

```>>> cl.cluster(maxclust=3)
>>> cl.R.data  # Restored centres of clusters
[L:146/39, L:43/26, L:323/59]
>>> StereoNet(*cl.groups)
``` ## Some tricks¶

Double cross products are allowed but not easy to understand.

For example `p**l**p` is interpreted as `p**(l**p)`: a) `l**p` is plane defined by `l` and `p` normal b) intersection of this plane and `p` is calculated:

```>>> p = Fol(250,40)
>>> l = Lin(160,25)
>>> StereoNet(p, l, l**p, p**l, l**p**l, p**l**p)
``` `Pair` class could be used to correct measurements of planar linear features which should spatialy overlap:

```>>> pl = Pair(250, 40, 160, 25)
>>> pl.misfit
18.889520432245405
>>> s = StereoNet()
>>> s.plane(Fol(250, 40), 'b', label='Original')
>>> s.line(Lin(160, 25), 'bo', label='Original')
>>> s.plane(pl.fol, 'g', label='Corrected')
>>> s.line(pl.lin, 'go', label='Corrected')
>>> s.show()
``` `StereoNet` has method `arrow` to draw arrow. Here is example of Hoeppner plot for variable fault orientation within given stress field:

```>>> S = Stress([[-8, 0, 0],[0, -5, 0],[0, 0, -1]]).rotate(Lin(90,45), 45)
>>> d = StereoGrid(npoints=300)
>>> s = StereoNet()
>>> s.tensor(S)
>>> for dc in d.dcgrid:
>>>     f = S.fault(dc)
>>>     s.arrow(f.fvec, f.lvec, f.sense)
>>> s.show()
``` 