Some geological utilities

copyright 2011 May 21 by Joshua R. Davis

Here are some generally useful utility functions for geology. Trends-and-plunges describe lines, and are given in degrees in geology coordinates (clockwise from north). Strikes-and-dips describe planes, and are given in degrees in geology coordinates (clockwise from north), according to the right-hand rule. Spherical coordinates (0 <= rho, 0 <= phi <= Pi, 0 <= theta <= 2 Pi) are in radians in math coordinates (counterclockwise from east). sphericalFromStrikeDipDeg returns the spherical coordinates of the pole to the plane described by the given strike and dip. sphericalFromCartesian always returns 0 <= rho, 0 <= phi <= pi, 0 <= theta <= 2 pi.

sphericalFromCartesian[{x_, y_, z_}] :=
	Module[{rho, phi, rhoSinPhi, theta},
		rho = Sqrt[x^2 + y^2 + z^2];
		If[rho == 0,
	{0, 0, 0},
	phi = ArcCos[z / rho];
	rhoSinPhi = rho Sin[phi];
	If[rhoSinPhi == 0,
				If[z >= 0,
					{z, 0, 0},
					{-z, Pi, 0}],
				{rho, phi, ArcTan[x, y]}]]];

sphericalFromTrendPlungeDeg[{trendDeg_, plungeDeg_}] :=
{1, (plungeDeg + 90.0) Degree, (90.0 - trendDeg) Degree};

sphericalFromStrikeDipDeg[{strikeDeg_, dipDeg_}] :=
{1, dipDeg Degree, (360 - strikeDeg) Degree};

cartesianFromSpherical[{rho_, phi_, theta_}] := 
{rho * Sin[phi] * Cos[theta], rho * Sin[phi] * Sin[theta], rho * Cos[phi]};

cartesianFromTrendPlungeDeg[x_] :=
cartesianFromSpherical[sphericalFromTrendPlungeDeg[x]];

cartesianFromStrikeDipDeg[x_] :=
cartesianFromSpherical[sphericalFromStrikeDipDeg[x]];

trendPlungeDegFromSpherical[{rho_, phi_, theta_}] :=
{90.0 - theta / Degree, phi / Degree - 90.0};

lowerHemisphere[{x_, y_, z_}] := 
{x, y, z} / If[z <= 0, Norm[{x, y, z}], -Norm[{x, y, z}]];

argMin[array_] :=
	Module[{i, j},
		If[VectorQ[array],
			If[Length[array] == 1,
				{1},
				i = argMin[array[[2;;]]];
				If[array[[1]] <= array[[i[[1]] + 1]],
					{1},
					{i[[1]] + 1}]],
			If[Length[array] == 1,
				Join[{1}, argMin[array[[1]]]],
				i = argMin[array[[1]]];
				j = argMin[array[[2;;]]];
				If[Apply[Part, Join[{array, 1}, i]] <= Apply[Part, Join[{array, j[[1]] + 1}, j[[2;;]]]],
					Join[{1}, i],
					Join[{j[[1]] + 1}, j[[2;;]]]]]]];

Here's an example of how to use argMin.

exampleArray = {{6, 8, -20}, {3, 4, 9}, {2, -27, 12}};
exampleArgMin = argMin[exampleArray]
Apply[Part, Join[{exampleArray}, exampleArgMin]]

The equal angle (i.e. stereographic) projection maps to a disk of radius 1; the equal area projection maps to a disk of radius sqrt(2). Both assume that their arguments lie on the unit sphere.

The hemisphericalPlot function draws an equal-area or equal-area hemispherical plot. The first argument is a list of points to be drawn as open circles. (All points are 3D Cartesian; this function pre-processes them to the unit lower hemisphere for you.) The second argument is a list of points to be drawn as filled circles. You can insert some Mathematica graphics commands into these point lists; for example, inserting Red switches the drawing color to red. The third argument is a list of curves; each curve is a list of points. You may not insert graphics commands into the curve list. There are also two options; the pointSize option is a positive real number, by default 0.02. The projection option is either equalAreaProjection (the default) or equalAngleProjection.

equalAngleProjection[{x_, y_, z_}] := 
	{x / (1 - z), y / (1 - z)};

equalAreaProjection[{x_, y_, z_}] := 
	Module[{scaling}, 
		If[x^2 + y^2 == 0, {0, 0},
			scaling = Sqrt[(2 (1 + z)) / (x^2 + y^2)];
			{scaling * x, scaling * y}]];

fringePrimitives[plotRadius_, padding_] :=
	{Circle[{0.0, 0.0}, plotRadius],
	Table[Line[{
			{plotRadius * Cos[angle], plotRadius * Sin[angle]},
{(plotRadius + padding) * Cos[angle], (plotRadius + padding) * Sin[angle]}}],
{angle, 0, 2 Pi, Pi / 18}],
	Text["000",{0,plotRadius + 2 * padding}],
	Text["090",{plotRadius + 2 * padding,0}],
	Text["180",{0,-(plotRadius + 2 * padding)}],
	Text["270",{-(plotRadius + 2 * padding),0}]};

pointPrimitives[projection_, pointList_, primitive_] :=
	Map[
		Function[item,
			If[ListQ[item],
			primitive[projection[lowerHemisphere[item]]],
			item]],
		pointList];

Options[hemisphericalPlot] = {pointSize -> 0.02, projection -> equalAreaProjection};

hemisphericalPlot[circleList_, diskList_, curveList_, OptionsPattern[]] :=
	Graphics[Join[
			{Black},
			fringePrimitives[If[OptionValue[projection] == equalAreaProjection, Sqrt[2], 1], 0.1],
			{Black},
			pointPrimitives[OptionValue[projection], circleList, Function[point, Circle[point, OptionValue[pointSize]]]],
			{Black},
			pointPrimitives[OptionValue[projection], diskList, Function[point, Disk[point, OptionValue[pointSize]]]],
			{Black},
			Map[Function[curve, Line[Map[OptionValue[projection], Map[lowerHemisphere, curve]]]], curveList]],
		AspectRatio -> 1];

Here's an example of plotting data. The filled data points are explicitly made up, in Cartesian coordinates. The unfilled data points are randomly generated in spherical coordinates. Two silly lines are also drawn, for no good reason.

black = {{1, 3, 2}, {2, 2, 2}, {5, 1, 3}, {4, 6, 2}, {1, 0, 2}, {0, 3, 1}, {2, 5, 0}, {3, 7, 4}, {-1, 2, 0}, {2, 5, 1}, {-2, 5, 2}, {3, 1, 1}, {0, -1, 3}, {1, -2, 1}, {-2, -1, -1}};
black = Map[lowerHemisphere, black];
white = Table[{1, RandomReal[{0, Pi}], RandomReal[{0, 2 Pi}]}, {i, 1, 10}];
white = Map[cartesianFromSpherical, white];
hemisphericalPlot[white, black, {{{1, 2, 3}, {4, 5, 6}}, {{-1, -2, 3}, {-4, -5, 6}, {-7, 8, 9}}}]

While we're at it, we might as well do Wulff and Schmidt nets. The routines below take two angles in radians. The first is the grid spacing. The second is an angle that determines how many samples to use in constructing the curves that make up the net. Make it equal to the grid spacing to get a net made up of straight lines; if your grid spacing is large, then this may look crude, so change the second angle to a smaller number.

lowerHemisphereGrid[gridSkip_, sampleSkip_] :=
	Module[{curves, curve, phi, theta},
		curves = {};
		(* Construct the parallels. *)
		phi = gridSkip;
		While[phi < Pi,
			curve = {};
			theta = 0;
			While[theta ≤ Pi,
				curve = Append[curve, {Sin[phi] Cos[theta], Cos[phi], -Sin[phi] Sin[theta]}];
				theta += sampleSkip];
			curves = Append[curves, curve];
			phi += gridSkip];
		(* Construct the meridians, including the two boundary meridians: *)
		theta = 0;
		While[theta <= Pi,
			curve = {};
			phi = 0;
			While[phi <= Pi,
				curve = Append[curve, {Sin[phi] Cos[theta], Cos[phi], -Sin[phi] Sin[theta]}];
				phi += sampleSkip];
			curves = Append[curves, curve];
	theta += gridSkip];
		curves];

wulffNet[gridSkip_, sampleSkip_] :=
hemisphericalPlot[{}, {}, lowerHemisphereGrid[gridSkip, sampleSkip], projection -> equalAngleProjection];

schmidtNet[gridSkip_, sampleSkip_] :=
hemisphericalPlot[{}, {}, lowerHemisphereGrid[gridSkip, sampleSkip]];

Here is an example:

wulffNet[10 Degree, 5 Degree]
schmidtNet[10 Degree, 5 Degree]