# PlaneTrendSurf: program that calculates an ASCII z-ccordinate ESRI raster for a known planar # orientation. # Input: ASCII text file "plane_in.txt". # Line 1. trend and plunge degrees of the true dip vector (trend in azimuth format). (ex. 135.0, 1.43) # Line 2. Elevation (Z value) of the measured attitude. (ex. 191.0) # Line 3. "X,Y" coordinates of the position of the attitude measurement (ex. "282883.069 , 251535.329"). # Line 4. SW corner (lower left) "X,Y" coordinates of the calculated raster image. # Line 5. "Columns,Rows" in calulated raster (ex. "251, 252"). # Line 6. Grid spacing between rows and columns of raster grid (ex. "21.6"). # ---------------------------------------------------------------------------------------------- # Output: ESRI raster file with following structure: # line 1: "ncols=" {number of columns in raster grid} # line 2: "nrows=" {number of rows in raster grid} # line 3: "xllcorner=" {x coordinate of the lower left (SW) corner of the grid} # line 4: "yllcorner=" {y coordinate of the lower left (SW) corner of the grid} # line 5: "cellsize=" {cell size spacing between columns and rows} # line 6: "nodata value=" {value that represents no data at grid node} #----------------------------------------------------------------------------------------------- import sys, numpy, math #import win32gui #enable WIN32 MessageBox dialog #import win32ui #import win32con # Function deg2rad(d): converts "d" from degrees to radians def deg2rad(d): return d * math.pi / 180.0 #Function rad2deg(r): converts "r" from radians to degrees def rad2deg(r): return r * 180.0 / math.pi # Function dprod(x,y,z,a,b,c): returns the dot-product of the directional components (x,y,z) and # (a,b,c). def dprod(x,y,z,a,b,c): return x*a + y*b + z*c # ------------------ Trend Surface Program --------------------------------------------------- # initialize constants tolerance = 1.0E-9 dl = " " eoln = "\n" # subdirectory path for output path = "C:/QGIS_Data/OutcropPrediction/Prob2/" # file name for output ofn = "Prob2_bot_grd.txt" rfn = "Prob2_bot_in.txt" # open trend output file in write-mode of = open(path+ofn,"w") rf = open(path+rfn,"r") # Enter azimuth, plunge of true dip: att_str = rf.readline() att_list = att_str.split(",") dipazstr = att_list[0] dipaz = float(dipazstr) plungestr = att_list[1] plunge = float(plungestr) # Calculate the pole to the plane directional components, and the directional components to the # X-Z plane if dipaz < 180.0 : pole_az = dipaz + 180.0 else: pole_az = dipaz - 180.0 pole_pl = 90.0 - plunge # directional components to pole to plane = (x1,y1,z1) x1 = math.sin(deg2rad(pole_az)) * math.sin(deg2rad(90.0-pole_pl)) y1 = math.cos(deg2rad(pole_az)) * math.sin(deg2rad(90.0-pole_pl)) z1 = math.cos(deg2rad(90.0-pole_pl)) # directional components of pole to X-Z plane = (x2,y2,z2) x2 = 0.0 y2 = 1.0 z2 = 0.0 # dot product yields theta angle between pole to plane and X-Z plane pole theta1 = math.acos(dprod(x1,y1,z1,x2,y2,z2)) if abs(theta1 - deg2rad(90.0)) < tolerance : theta1 = deg2rad(90.0) # cross product of pole to plane and X-Z plane pole yields the line of intersection # between the 2 planes (x3,y3,z3), and this line is the slope of the plane in the X-Z plane c1_x = (y1*z2-z1*y2)/math.sin(theta1) c1_y = -(x1*z2-z1*x2)/math.sin(theta1) c1_z = (x1*y2-y1*x2)/math.sin(theta1) if c1_z < 0.0 : c1_x = c1_x * -1.0 c1_y = c1_y * -1.0 c1_z = c1_z * -1.0 # check for round off errors if c1_x > 1.0 : c1_x = 1.0 if c1_x < -1.0 : c1_x = -1.0 c1 = -math.tan(math.acos(c1_x)) # directional components of pole to Y-Z plane = (x2,y2,z2) x2 = -1.0 y2 = 0.0 z2 = 0.0 # dot product yields theta angle between pole to plane and Y-Z plane pole theta2 = math.acos(dprod(x1,y1,z1,x2,y2,z2)) # cross product of pole to plane and Y-Z plane pole yields the line of intersection # between the 2 planes (x3,y3,z3), and this line is the slope of the plane in the Y-Z plane c2_x = (y1*z2-z1*y2)/math.sin(theta2) c2_y = -(x1*z2-z1*x2)/math.sin(theta2) c2_z = (x1*y2-y1*x2)/math.sin(theta2) if c2_z < 0.0 : c2_x = c2_x * -1.0 c2_y = c2_y * -1.0 c2_z = c2_z * -1.0 # check for round off errors if c2_y > 1.0 : c2_y = 1.0 if c2_y < -1.0 : c2_y = -1.0 c2 = -math.tan(math.acos(c2_y)) # the elevation (z) value of the measured attitude: elev_str = rf.readline() elev = float(elev_str) # the (x,y) coordinates of the measured attitude: xy_str = rf.readline() crd_list = xy_str.split(",") x_crd = float(crd_list[0]) y_crd = float(crd_list[1]) c0=elev-c1*x_crd-c2*y_crd # (x,y) coordinates of grid SW corner: sw_xy_str = rf.readline() crd_list = sw_xy_str.split(",") sw_x_crd = float(crd_list[0]) sw_y_crd = float(crd_list[1]) # number of columns,rows for grid: grid_str = rf.readline() grid_list = grid_str.split(",") ncols_str = grid_list[0] ncols = int(ncols_str) nrows_str = grid_list[1] nrows = int(nrows_str) # cell size of raster grid: cellsz_str = rf.readline() cellsize = float(cellsz_str) nodatavalue = 1e+30 # write header info of.write("ncols" + dl + str(ncols) + eoln) of.write("nrows" + dl + str(nrows) + eoln) of.write("xllcorner" + dl + str(sw_x_crd) + eoln) of.write("yllcorner" + dl + str(sw_y_crd) + eoln) of.write("cellsize" + dl + str(cellsize) + eoln) of.write("nodata_value" + dl + str(nodatavalue) + eoln) linect = 6 i = 0 while i < nrows: y = sw_y_crd + (nrows - i) * cellsize j = 0 while j < ncols: x = sw_x_crd + (j * cellsize) z = c0 + c1*x + c2*y of.write(str(z)+eoln) linect = linect + 1 j = j + 1 i = i + 1 titlestr = "Results Message" + eoln print titlestr msgstr = "Planar surface raster written to: " + path + ofn + eoln print msgstr # win32ui.MessageBox(msgstr,titlestr,win32con.MB_OK) titlestr = "z = c0 + c1*x + c2*y" + eoln print titlestr msgstr = "c0=" + str(c0) + eoln + "c1=" + str(c1) + eoln + "c2=" + str(c2) + eoln print msgstr # win32ui.MessageBox(msgstr,titlestr,win32con.MB_OK) # close the files of.close() rf.close()