My previous post showed the mathematics behind computationally calculating strike and dip from 3 points. Here I will outline how this may be accomplished in a Python script. This is written with Python3 in mind, but should work in Python2 as well.
If you get lost with the math, check back to my post on the mathematics behind this.
First, we’re going to import math since we’ll be doing some math.
The appropriate input for this function is a list of tuples in the format [(x1, y1, z1), (x2, y2, z2), (x3, y3, z3)]. I define everything in functions so that they can be easily slotted into other programs down the road. The first step is to read this list and place the values into appropriate variables. I specifically designate float() in case the list is parsed out of a text file or entered through the input() method.
def calc_strikedip(pts): ptA, ptB, ptC = pts, pts, pts x1, y1, z1 = float(ptA), float(ptA), float(ptA) x2, y2, z2 = float(ptB), float(ptB), float(ptB) x3, y3, z3 = float(ptC), float(ptC), float(ptC)
Now we define the components u1, u2, and u3 that we get by the cross product of the vectors AB, AC.
u1 = float(((y1-y2)*(z3-z2)-(y3-y2)*(z1-z2))) u2 = float((-((x1-x2)*(z3-z2)-(x3-x2)*(z1-z2)))) u3 = float(((x1-x2)*(y3-y2)-(x3-x2)*(y1-y2)))
These represent pseudo eastings and northings, these are actually coordinates of a new point that represents the normal from the plane’s origin defined as (0,0,0).
If the z value (u3) is above the plane we first reverse the easting then we check if the z value (u3) is below the plane, if so we reverse the northing. I break up the calculation of strike into partA_strike first as the denominator to make it easier for my brain to follow it, it is not strictly necessary.
This is to satisfy the right hand rule in geology where dip is always to the right if looking down strike.
if u3 < 0: easting = u2 else: easting = -u2 if u3 > 0: northing = u1 else: northing = -u1 if easting >= 0: partA_strike = math.pow(easting, 2) + math.pow(northing, 2) strike = math.degrees(math.acos(northing / math.sqrt(partA_strike))) else: partA_strike = northing / math.sqrt(math.pow(easting, 2) + math.pow(northing, 2)) strike = math.degrees(2 * math.pi - math.acos(partA_strike))
When calculating strike, we check if the easting is greater than or equal to 0 to determine if strike should be reversed to obey the right hand rule. If it isn’t then we simply reverse the reading by subtracting the the strike angle calculated from 2π in radians before converting to degrees.
Now we determine dip, again these are split up into numerator and denominator for ease of understanding.
part1_dip = math.sqrt(math.pow(u2, 2) + math.pow(u1, 2)) part2_dip = math.sqrt(math.pow(u1,2) + math.pow(u2,2) + math.pow(u3,2)) dip = math.degrees(math.asin(part1_dip / part2_dip)) return strike, dip
Hopefully this helps someone out. You can download the complete script from my git repository: