Fisher Statistics, Strike & Dip Version

FisherStatsSD.nb is available online at http://CroninProjects.org/Vince/Codes/FisherStatsSD.nb
© 2011-15 by Vincent S. Cronin.  This code may be used for non - profit educational and scientific purposes with appropriate attribution to the original author (Vincent S. Cronin).

Version 1.0.  Begun January 21, 2011;  revised December 1, 2015

Introduction

Mathematica code to use Fisher statistics (Fisher, 1953) to find the mean and 95% confidence interval of a set of data that characterize several spot orientations of a surface.  The input data are in an Excel spreadsheet, and consist of strikes and dip angles.

Data

Description of input data

Fisher statistics are used to find the mean and characterize the dispersion around the mean of a set of vectors.  They can also be used to find the mean orientation of planes by using dip vectors or the vectors normal to the planes as input data.  The input data include the right-hand-rule strike azimuth, in degrees (range of azimuth is 0 to 360°), and the plunge of the dip vector in degrees (range of dip angle is 0 to 90°).

This notebook assumes that the input values will be integers, because the uncertainty inherent in a Brunton compass or other pocket transit is on the order of 1°.  Hence, the output data are in integers.

Making a data file in Excel

The input data is an Excel spreadsheet (.xls file) in which each record/row has 2 values/columns.  The first column contains the right-hand-rule strike azimuth, in degrees, and the second column has the plunge of the dip vector (i.e., the dip angle), in degrees.  The right-hand-rule strike azimuth is 90° anticlockwise from the trend of the dip vector.

Finding the Excel data file on your computer

Each user will need to ensure that the path to the input data file is correctly specified in the first input line of this notebook (i.e., in the blue box below).  For example, a correct specification for the file "rawInputData.xls" located on the desktop of Vince Cronin's office iMac computer would look like this:

mydata = Import["/Users/vince_cronin/Desktop/rawData.xls"];

and the specification for the same file located in the StructGeol directory on the C drive of Cronin's Dell computer would look like this:

mydata=Import["C:\StructGeol\rawData.xls"];

In the code immediately following this text, be certain that the path to the input data file is correctly specified.

mydata=Import["/Users/vince_cronin/Desktop/rawData.xls"];

inData0=Flatten[mydata,1] ;

Computation

Functions defined in this notebook

The function findVector converts orientation data for a vector, expressed as the trend and plunge of the vector, to the corresponding 3 D Cartesian unit vector coordinates.

findVector[vectTrend_, vectPlunge_]:={Sin[vectTrend Degree]Cos[vectPlunge Degree],Cos[vectTrend Degree]Cos[vectPlunge Degree],-Sin[vectPlunge Degree]};

The function cart2TrendPlunge converts from a 3 D Cartesian unit vector to trend and plunge.  A downward plunge has a positive sign.

cart2TrendPlunge[inVector_]:=Module[{a,unitA,trend,plunge},plunge=-ArcSin[inVector[[3]]](180/π);a={inVector[[1]],inVector[[2]]};unitA={a[[1]]/Norm[a],a[[2]]/Norm[a]};trend=If[(unitA[[1]]<0),(360-(ArcCos[unitA[[2]]](180/π))),(ArcCos[unitA[[2]]](180/π))];{trend,plunge}];

The module equalAreaPlot computes the {x,y} coordinates for an equal-area stereoplot of a vector whose orientation is given in terms of trend and plunge expressed in degrees.   This module is based on the assumption that the stereonet has a unit radius (i.e., radius = 1).  The first element of inTrendPlunge is the trend (azimuth measured clockwise from north), and the second is the plunge angle (measured from horizontal, positive down, negative up).

equalAreaPlot[inTrendPlunge_]:=Module[{radius,xCoord,yCoord},radius=1;xCoord=radius(Sqrt[2])Sin[((π/4))-((((inTrendPlunge[[2]] ((π/180)))/2)))]Sin[inTrendPlunge[[1]] ((π/180))];yCoord=radius(Sqrt[2])Sin[((π/4))-((((inTrendPlunge[[2]] ((π/180)))/2)))]Cos[inTrendPlunge[[1]] ((π/180))];{xCoord,yCoord}];

The function round2Int rounds an arbitrary number to the nearest integer value.

round2Int[x_]:=Round[IntegerPart[x*10]/10];

Define the sample size

matrixSize1=Dimensions[inData0];

Variable noInputVectors is the number of vectors in the set.  It is equivalent to the variable "N" in Cronin (2008).

noInputVectors=matrixSize1[[1]];

Convert strike & dip data to dip vector data

The following step converts the input strike-and-dip data to trend-and-plunge vector data for the corresponding dip vectors.  The input strike data must be expressed using the right - hand - rule convention.

inData=Table[{If[(inData0[[i,1]]<270),(inData0[[i,1]]+90),(inData0[[i,1]]-270)],inData0[[i,2]]},{i,noInputVectors}];

Proceed with vector analysis

Variable dirCosines is a table of the direction cosines { l, m, n} for each of the input vectors (equation 1 from Cronin, 2008).

dirCosines=Table[{Cos[inData[[i,2]] Degree]Cos[inData[[i,1]] Degree],Cos[inData[[i,2]] Degree]Sin[inData[[i,1]] Degree],Sin[inData[[i,2]] Degree]},{i,noInputVectors}];

Variable meanDirCos includes the set of three direction cosines for the mean dip vector (equation 2 from Cronin, 2008).

meanDirCos=Total[dirCosines,1];

Variable capR is the length of the mean dip vector (equation 3 from Cronin, 2008).  It is equivalent to the variable "R" in Cronin (2008).

capR=(Sqrt[((meanDirCos[[1]])^2)+((meanDirCos[[2]])^2)+((meanDirCos[[3]])^2)]);

Variable meanUnitVector is the set of coordinates for the mean unit dip vector (equation 4 from Cronin, 2008).

meanUnitVector={((meanDirCos[[1]])/capR),((meanDirCos[[2]])/capR),((meanDirCos[[3]])/capR)};

Variable delta is the plunge angle of the mean dip vector in radians (equation 5 from Cronin, 2008).  It is equivalent to the variable "δ" in Cronin (2008).

delta=ArcSin[meanUnitVector[[3]]];

Variable meanDipTrend is the trend of the mean dip vector in degrees (equation 6 from Cronin, 2008).

meanDipTrend=If[(meanUnitVector[[2]]<0),(((2π)-ArcCos[((meanUnitVector[[1]])/(Cos[ArcSin[meanUnitVector[[3]]]]))])(180/π)),((ArcCos[((meanUnitVector[[1]])/(Cos[ArcSin[meanUnitVector[[3]]]]))])(180/π))];

Variable meanStrike is the trend of the mean strike azimuth measured clockwise from north in degrees (equation 6 from Cronin, 2008), determined using the right-hand-rule convention -- the RHR strike is 90° anticlockwise from the trend of the dip vector.

meanStrike=If[(meanDipTrend<90),(meanDipTrend+270),(meanDipTrend-90)];

Variable k is the precision parameter, which ranges from 0 for a vector set that is strongly noncolinear to infinity for vectors that are perfectly colinear, as when N = R (equation 7 from Cronin, 2008).

k=((noInputVectors-1)/(noInputVectors-capR));

Variable alpha95 is the angular radius of the 95% confidence cone around the mean vector, in radians (equation 8 from Cronin, 2008).  It is equivalent to the variable "(Subscript[α, 95])" in Cronin (2008).

alpha95=ArcCos[1-((((noInputVectors-capR)/capR))(((((1/0.05)))^(1/(noInputVectors-1)))-1))];

Variable theta is an intermediate value in the computation of the uncertainty in the azimuth of the vector (equation 10 from Cronin, 2008).  It is equivalent to the variable "θ" in Cronin (2008).

theta=ArcSin[((Sin[alpha95]Sin[delta])/(Cos[alpha95]Cos[delta]))];

Variable beta is the uncertainty in the azimuth of the vector (equation 12 from Cronin, 2008).  It is equivalent to the variable "β" in Cronin (2008).

beta=ArcTan[((Sin[alpha95]Cos[theta])/((Cos[alpha95]Cos[delta])-(Sin[alpha95]Sin[theta]Sin[delta])))];

Round numerical output to integers

meanDipAzimuth=round2Int[meanDipTrend];

meanStrike=round2Int[meanStrike];

strikeUncertainty=round2Int[beta (180/π)];

meanDipAngle=round2Int[delta(180/π)];

alphaNinetyFive=round2Int[alpha95 (180/π)];

precisionParameterK=round2Int[k];

Graphics

Computation

meanLine={meanDipAzimuth,meanDipAngle};

xMatrix1=(({{1, 0, 0}, {0, Cos[(90-meanDipAngle) Degree], -Sin[(90-meanDipAngle) Degree]}, {0, Sin[(90-meanDipAngle) Degree], Cos[(90-meanDipAngle) Degree]}}));

zMatrix1=(({{Cos[-meanDipAzimuth Degree], -Sin[-meanDipAzimuth Degree], 0}, {Sin[-meanDipAzimuth Degree], Cos[-meanDipAzimuth Degree], 0}, {0, 0, 1}}));

Define a set of points in a circle oriented with a radius of alpha 95 degrees around a vector pointing straight down.

vectorPlunge=90-(round2Int[alpha95 (180/π)]);

initSmCircle=Table[{i,vectorPlunge},{i,0,360,10}];

matrixSize2=Dimensions[initSmCircle];

Variable ptsAroundCircle is the number of points used to define the small circle.

ptsAroundCircle=matrixSize2[[1]];

initSmCircCart=Table[findVector[initSmCircle[[i,1]],initSmCircle[[i,2]]],{i,ptsAroundCircle}];

rotSmCircCart0 =Table[{zMatrix1.xMatrix1.initSmCircCart[[ i]]}, {i, ptsAroundCircle}];

rotSmCircCart=Flatten[rotSmCircCart0,1];

rotSmCircTP=Table[cart2TrendPlunge[rotSmCircCart[[i]]],{i,ptsAroundCircle}];

plotSmCircDataFile=Table[equalAreaPlot[rotSmCircTP[[i]]],{i,ptsAroundCircle}];

part1=Graphics[Circle[{0,0},1]];

part2=Graphics[Line[{{{0,0.01},{0,-0.01}},{{0.01,0},{-0.01,0}},{{0,1},{0,0.95}},{{0,-1},{0,-0.98}},{{1,0},{0.98,0}},{{-1,0},{-0.98,0}}}]];

See also PlotMarkers under ListPlot

part3=Graphics[Line[plotSmCircDataFile]];

part4=Graphics[{PointSize[Small],Red,Point[equalAreaPlot[meanLine]]}];

rawPoints=Table[equalAreaPlot[inData[[i]]],{i,noInputVectors}];

part5=Graphics[{PointSize[Small],Point[rawPoints]}];

Plot the plane and the pole to the plane

Define a set of points that will trace a great - circle arc on the final stereo plot

rawGreatCircle={{270,0},{270,10},{270,20},{270,30},{270,40},{270,50},{270,60},{270,70},{270,80},{270,90},{90,80},{90,70},{90,60},{90,50},{90,40},{90,30},{90,20},{90,10},{90,0}};

ptsAlongGrCirc=19;

initGrCircCart=Table[findVector[rawGreatCircle[[i,1]],rawGreatCircle[[i,2]]],{i,ptsAlongGrCirc}];

rotGrCircCart0=Table[{zMatrix1.xMatrix1.initGrCircCart[[ i]]},{i,ptsAlongGrCirc}];

rotGrCircCart=Flatten[rotGrCircCart0,1];

rotGrCircTP=Table[cart2TrendPlunge[rotGrCircCart[[i]]],{i,ptsAlongGrCirc}];

plotGrCircDataFile=Table[equalAreaPlot[rotGrCircTP[[i]]],{i,ptsAlongGrCirc}];

part6=Graphics[Line[plotGrCircDataFile]];

xMatrix2=(({{1, 0, 0}, {0, Cos[-(meanDipAngle) Degree], -Sin[-(meanDipAngle) Degree]}, {0, Sin[-(meanDipAngle) Degree], Cos[-(meanDipAngle) Degree]}}));

zMatrix2=(({{Cos[-meanDipAzimuth Degree], -Sin[-meanDipAzimuth Degree], 0}, {Sin[-meanDipAzimuth Degree], Cos[-meanDipAzimuth Degree], 0}, {0, 0, 1}}));

rotSmCircCartB =Table[{zMatrix2.xMatrix2.initSmCircCart[[ i]]}, {i, ptsAroundCircle}];

rotSmCircCart2=Flatten[rotSmCircCartB,1];

rotSmCircTP2=Table[cart2TrendPlunge[rotSmCircCart2[[i]]],{i,ptsAroundCircle}];

plotSmCircDataFile2=Table[equalAreaPlot[rotSmCircTP2[[i]]],{i,ptsAroundCircle}];

temp3=Line[plotSmCircDataFile2];

part7=Graphics[{Dashed,Blue,temp3}];

meanPole={If[(meanDipAzimuth<180),meanDipAzimuth+180,meanDipAzimuth-180],(90-meanDipAngle)};

part8=Graphics[{PointSize[Small],Blue,Point[equalAreaPlot[meanPole]]}];

Output

The numerical output below is rounded to the nearest integer value.  All angles and azimuths are expressed in degrees.

meanDipAzimuth

meanStrike

strikeUncertainty

meanDipAngle

alphaNinetyFive

precisionParameterK

The graphics below are lower - hemisphere equal - area projections, with the longer tick mark at the top of the circle directed north and the cross in the middle of the circle (i.e., at the bottom of the projection hemisphere).

The stereo plot below includes a small circle with a radius of alpha95 ((Subscript[α, 95])) degrees showing the projection of the cone - shaped uncertainty region around the mean dip vector (the red dot).

Show[part1,part2,part3,part4]

The stereo plot below includes the dip vectors of the input planar data (black dots) in addition to the features of the previous plot.

Show[part1,part2,part3,part4,part5]

The stereo plot below includes the trace (i.e., a great - circle arc) of the mean/average plane in addition to the features of the previous plot.

Show[part1,part2,part3,part4,part5,part6]

The stereo plot below includes the pole (blue dot) of the mean/average plane and the associated 95% CI uncertainty region (blue dashed circle/ellipse), in addition to the features of the previous plot.

Show[part1,part2,part3,part4,part5,part6,part7,part8]

References Cited and Other Relevant Texts

Batschelet, E., 1981, Circular Statistics in Biology:  London, Academic Press, 371 p.

Borradaile, G., 2003, Statistics of Earth Science Data, Their Distribution in Time, Space, and Orientation:  Berlin, Springer-Verlag, 351 p.

Brunton Company, 2000, Pocket Transit Instruction Manual, The Brunton Company Riverton, WY. 25 p.

Compton, R.R., 1962, Manual of Field Geology:  New York, John Wiley & Sons, 378 p.

Cronin, V.S., 2008, Finding the mean and 95 % confidence interval of a set of strike - and - dip or lineation data: Environmental and Engineering Geoscience, v.14, no.2, p.113 - 119.

Fisher, N.I., Lewis, T., and Embleton, B.J.J., 1987, Statistical Analysis of Spherical Data:  Cambridge, UK, Cambridge University Press, 329 p.

Fisher, R.A., 1953, Dispersion on a sphere: Proceedings Royal Society, London, v. A217, no. 1130. p. 295-305.

Irving, E., 1964, Paleomagnetism and its Application to Geological and Geophysical Problems:  New York, John Wiley & Sons, 399 p.

Mardia, K.V., 1972, Statistics of Directional Data:  London, Academic Press, 357 p.

Opdyke, N.D., and Channell, J.E.T., 1996, Magnetic Stratigraphy:  San Diego, California, Academic Press, 346 p.

Tarling, D.H., 1971, Principles and Applications of Palaeomagnetism:  London, Chapman and Hall, 164 p.

Tauxe, L., 1998, Paleomagnetic Principles and Practices:  Dordrecht, Kluwer Academic Publishers, 299 p.