This is my first time coding in python, so bear with me. I wrote a script to calculate an analytical solution to find the orientation of a plane based on three points. I'm using ArcGIS 10.1.
from __future__ import division
import arcpy
import numpy as np
# Input points
x1 = float(arcpy.GetParameterAsText(1)) # Point 1
y1 = float(arcpy.GetParameterAsText(2))
z1 = float(arcpy.GetParameterAsText(3))
x2 = float(arcpy.GetParameterAsText(4)) # Point 2
y2 = float(arcpy.GetParameterAsText(5))
z2 = float(arcpy.GetParameterAsText(6))
x3 = float(arcpy.GetParameterAsText(7)) # Point 3
y3 = float(arcpy.GetParameterAsText(8))
z3 = float(arcpy.GetParameterAsText(9))
# Plane equation
A = (y1*z2) + (z1*y3) + (y2*z3) - (z2*y3) - (z3*y1) - (z1*y2)
B = (z2*x3) + (z3*x1) + (z1*x2) - (x1*z2) - (z1*x3) - (x2*z3)
C = (x1*y2) + (y1*x3) + (x2*y3) - (y2*x3) - (y3*x1) - (y1*x2)
D = (z1*y2*x3) + (z2*y3*x1) + (z3*y1*x2) - (x1*y2*z3) - (y1*z2*x3) - (z1*x2*y3)
E = np.sqrt(A**2 + B**2 + C**2)
AZ = np.arctan(A/B) * (180/np.pi) # Preliminary azimuth for dip direction
dip = np.arcsin(-np.cos((np.pi/2) + np.arccos(C/E))) * (180/np.pi) # True dip
# Test for quadrant
alpha = A/E
beta = B/E
# Place dip-direction in proper quadrant
if alpha > 0 and beta > 0:
DD = AZ
elif alpha > 0 and beta < 0:
DD = 180 + AZ
elif alpha < 0 and beta < 0:
DD = 180 + AZ
elif alpha < 0 and beta > 0:
DD = 360 + AZ
# Right-hand rule strike
if DD - 90 < 0:
RHR = 360 + (DD - 90)
else:
RHR = DD - 90
# Centroid coordinates
Cx = (x1+x2+x3)/3
Cy = (y1+y2+y3)/3
Cz = (z1+z2+z3)/3
# Insert centroid point
cursor = arcpy.InsertCursor(arcpy.GetParameterAsText(0)) # Feature class must already exist
row = cursor.newRow()
row.Shape = arcpy.Point(Cx,Cy,Cz)
# Populate table
row.X = Cx
row.Y = Cy
row.Z = Cz
row.DD = DD
row.dip = dip
row.RHR = RHR
cursor.insertRow(row)
del cursor, row
# Refresh MXD
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
This works, except that the user has to hover each point and manually enter the coordinates. Ideally, I would like: (1) the user initially define the target feature layer in a pop-up menu, and hit OK; (2) click on the map to extract the first XY map coordinates; (3) pop-up queries the user for the first Z, hit OK; (4) click on the map to extract the second XY map coordinates; (5) pop-up queries the user for the second Z, hit OK; (6) click on the map to extract the last XY map coordinates; (7) pop-up queries the user for the last Z, hit OK; then(8) the script runs to completion.I found the ESRI tool tutorial that allows the map XY to be displayed on a message box after a mouse-click, but cannot figure out how to pass the variable.
