Python Programming Project - University of South Alabama



Python Programming ProjectIntroductionIn this project we will learn some of the fundamentals of programming in the ArcPython programming environment. The project will consist of 2 different assignments of increasing difficulty:Create a python program that will take the x,y coordinates of 2 points on the command line and calculate the azimuth direction from point 1 to point 2 in degrees, and print the azimuth direction as the resultCreate a python program that can use the coefficients of a 1st, 2nd, and 3rd order trend surface to generate a 3D grid file that can be imported and used with ArcGIS. In addition, an Excel spreadsheet will need to be created that will calculate the polynomial coefficients of the 1st, 2nd, or 3rd order trend surface. The python program should use a command line parameter to control the order of the trend surface output. All of the programs should be adequately documented (%20 of grade) so that the logic of the program is evident. Azimuth CalculatorStep 1: import the “arcpy”, “math”, and “sys” librariesStep 2: Convert the command line strings to numeric values# command line strings are extracted as belowX1param = float(sys.argv[1])# note the “float()” casting converts the string to a floating point number# do the same for y2param, x2param, y2paramNote that if you are copying and pasting from ArcMap that the coordinates will normally have a coma (“,”) separator between hundreds and thousands places, etc. for UTM coordinates. This comma will need to be removed or replaced with a blank before applying the float() function. Also note that unlike other Python counts such as lists the first item on the command line list begins with a 1 index (sys.argv[1]) rather than a 0 index (sys.argv[0]). Step 3: Calculate the change in x and y from (x1,y1) to (x2,y2) and store in variables “dx” and “dy”Step 4: Use “atan2(dy,dx)” to calculate the azimuth angle. You will need an “if” statement to correctly calculate the angle. The “atan2” returns the angle in radians so the result of the function will have to be multiplied by 180/math.pi. The possibility of dx and dy = 0 should be trapped and the azimuth set to 0 in that case. Note that atan2() returns the angle result following these constraints-{in pseudo-code}If dx=0.0 and dy = 0.0 then azimuth = 0 Elseif dx <0.0 and dy >= 0.0 then azimuth = 450 – atan2 (dy,dx) * 180/piElse azimuth = 90 - atan2(dy,dx) * 180/piStep 5: Print the azimuth result so that it appears as the result in the interactive window.As the ultimate test of your azimuth calculator you should be able to use the “identity” tool to capture the (x,y) coordinates of 2 points on a map (use the Talladega Springs project map), feed those coordinates to the command line of the program to calculate the azimuth. Make sure that it “makes sense” relative to your map data. For example, if point 2 is nearly due east of point 1, then feeding point 1 and point 2 to your calculator should yield an azimuth near 90.0 degrees.Of course you don’t have to go to the trouble of loading a big map into ArcMap to test your program. After you have your script typed and saved choose “File > Run” from the main menu and use “0.0 0.0 1.0 1.0” on the command line. The result should be 45.0 degrees. A command line of “0.0 0.0 0.1 1.0” should yield 5.7 degrees. Check for azimuth direction in all quadrants (NE, SE, SW, and NW). Check for special directions such as due west (270) or east (90). Download the “Azimuth.py” program from “Resources” to see a complete program script. Useful Functions and Methods for Azimuth CalculatorThe below functions and methods will prove useful for this problem. To use them make sure to import the arcpy, sys, and math libraries:sys.argv[i]returns the ith parameter of the command line entered when starting the program.float(s)returns the numeric equivalent of the string stored in s. The result of float(“1.57”) would be a floating point number = 1.57.math.atan2(x,y):returns the arctangent angle of the ratio of y/x in radians. The range of atan2(x,y) is pi to –pi (180 to -180 degrees. print obj:prints the contents object “obj” to the Python interactive window.strvar.replace(target,replacement):Every string variable automatically has the method “replace” attached to it so a statement such as strvar.replace(“,”,”“) would replace every occurrence of “,” with a null character (i.e. nothing). strvar.format(number): this string method can precisely format a number within a string. For example if x = 6.456978 the statement:print “:{8.3f}”.format(x) would generate:>>> 6.457Trend Surface ProgrammingThe trend surface project will consist of 2 major components:An Excel spreadsheet that will contain the raw X,Y,Z data and from that calculate the trend surface coefficients for a 1st, 2nd, and 3rd order surface fit. The coefficients will be used in the second part of the project, as will the raw data. The starting spreadsheet containing the well data in the first sheet (“Data” sheet) will be provided. See the “Resources” section of the GIT461 site and download the “TrendSurfTemplate.xlsm” file. In this file the summations and matrix operations needed for the 1st order trend surface are already worked out so you need to add the additional calculations needed for the 2nd and 3rd order trend surfaces.A python program that will use the trend surface coefficients calculated in the Excel spreadsheet to compute and write an ArcGIS-compatible raster grid to a text file. The user should be able to indicate the order (1,2,3) of the surface via a command line parameter. A template program can be downloaded from “Resources” as “TrendSurface_1.py”.The format of a ESRI ASCII raster format is given below:Example ASCII raster:ncols 480nrows 450xllcorner 378923yllcorner 4072345cellsize 30nodata_value -3276843 2 45 7 3 56 2 5 23 65 34 6 32 54 57 34 2 2 54 6 35 45 65 34 2 6 78 4 2 6 89 3 2 7 45 23 5 8 4 1 62 ...Your python program will construct an output file based on a trend surface polynomial that follows the above format.Excel Spreadsheet for Calculating the Trend Surface CoefficientsThe spreadsheet should contain in the first sheet (“Data”) the raw data in these columns:Column A: No.{data number}Column B: X{x coordinate}Column C: Y{y coordinate}Column D: Z{z coordinate}Column E: X^2{x * x)Column F: XY{x * y}Column G: Y^2{y * y)Column H: X^2Y{x * x * y}Column I: XY^2{x * y * y}Column J: XZ{x * z}Column K: YZ{y * z}Column L: X^3{x * x * x}Column M: Y^3{y * y * y}Column N: X^4{x * x * x * x}Column O: Y^4{y * y * y * y}Column P: X^3Y{x * x * x * y}Column Q: X^2Y^2{x * x * y * y}Column R: XY^3{x * y * y * y}Column S: X^Z{x * x * z}Column T: XYZ{x * y * z}Column U: Y^2Z{y * y * z}Column V: X^5{x * x * x * x * x}Column W: X^4Y{x * x * x * x * y}Column X: X^3Y^2{x * x * x * y * y}Column Y: X^2Y^3{x * x * y * y * y}Column Z: XY^4{x * y * y * y * y}Column AA: Y^5{y * y * y * y * y}Column AB: X^6{x * x * x * x * x * x}Column AC: X^5Y{x * x * x * x * x * y}Column AD: X^4Y^2{x * x * x * x * y * y}Column AE: X^3Y^3{x * x * x * y * y * y}Column AF: X^2Y^4{x * x * y * y * y * y}Column AG: XY^5{x * y * y * y * y * y}Column AH: Y^6{y * y * y * y * y * y}Column AI: X^3Z{x * x * x * z}Column AJ: X^2YZ{x * x * y * z}Column AK: XY^2Z{x * y * y * z}Column AL: Y^3Z{y * y * y * z}At the bottom of each column the sum of all of the cells in the column should be calculated with the “sum()” function. Each of the summations should be named “Sum_x” for summation of “x” column, “Sum_X2” for summation of “X^2” column, and so on. Use the name manager to do this (“Formulas” menu > “Name Manager” button). The number of data observations should also be counted and the cell named “N”. This should be the last observation cell in the “No.” column.The second sheet should be named “TrendSurf” and should contain the following matrices:Matrix A1(D3:F5)NSum_XSum_YSum_XSum_X2Sum_XYSum_YSum_XYSum_Y2Matrix B1(H3:H5)Sum_ZSum_XZSum_YZInverse A1(D8:F10)=MINVERSE(D3:F5)1st Order Coeff.(D13:D15)=MMULT(D8:F10,H3:H5)t1_c0t1_c1t1_c2 z(x,y) =c0 + x*c1 + y*c2Note that the trend surface coefficients are calculated by matrix multiplication of the A1 inverse x B1. This will be true also for the 2nd and 3rd order trend surface calculations, the only difference being that the matrices are larger. When using the “minverse()” function highlight the destination cells before selecting it from the “Math & Trig” function set in “Formulas”, indicate the “A” matrix block of cells as input, and then finish the formula dialog window with a “ctrl”+”shift”+”enter” to fill the destination matrix. The same procedure should be used with “mmult()” when multiplying the inverse A and B matrices to calculate the trend surface coefficients.For the 2nd order trend surface:Matrix A2(D18:I23)NSum_XSum_YSum_X2Sum_XYSum_Y2Sum_XSum_X2Sum_XYSum_X3Sum_X2YSum_XY2Sum_YSum_XYSum_Y2Sum_X2YSum_XY2Sum_Y3Sum_X2Sum_X3Sum_X2YSum_X4SumX3YSum_X2Y2Sum_XYSum_X2YSum_XY2Sum_X3YSum_X2Y2Sum_XY3Sum_Y2Sum_XY2Sum_Y3Sum_X2Y2Sum_XY3Sum_Y4Matrix B2(K18:K23)Sum_ZSum_XZSum_YZSum_X^2ZSum_XYZSum_Y^2ZInverse A2(D26:I31)=MINVERSE(D18:I23)2st Order Coeff.(D13:D15)=MMULT(D26:I31,K18:K23)t2_c0t2_c1t2_c2t2_c3t2_c4t2_c5Z(x,y) = c0 + x*c1 + y*c2 + x^2*c3 + x*y*c4 + y^2*c5For the 3rd order trend surface:Matrix A3(D42:M51)NSum_XSum_YSum_X2Sum_XYSum_Y2Sum_X3Sum_X2YSum_XY2Sum_Y3Sum_XSum_X2Sum_XYSum_X3Sum_X2YSum_XY2Sum_X4Sum_X3YSum_X2Y2Sum_XY3Sum_YSum_XYSum_Y2Sum_X2YSum_XY2Sum_Y3Sum_X3YSum_X2Y2Sum_XY3Sum_Y4Sum_X2Sum_X3Sum_X2YSum_X4Sum_X3YSum_X2Y2Sum_X5Sum_X4YSum_X3Y2Sum_X2Y3Sum_XYSum_X2YSum_XY2Sum_X3YSum_X2Y2Sum_XY3Sum_X4YSum_X3Y2Sum_X2Y3Sum_XY4Sum_Y2Sum_XY2Sum_Y3Sum_X2Y2Sum_XY3Sum_Y4Sum_X3Y2Sum_X2Y3Sum_XY4Sum_Y5Sum_X3Sum_X4Sum_X3YSum_X5Sum_X4YSum_X3Y2Sum_X6Sum_X5YSum_X4Y2Sum_X3Y3Sum_X2YSum_X3YSum_X2Y2Sum_X4YSum_X3Y2Sum_X2Y3Sum_X5YSum_X4Y2Sum_X3Y3Sum_X2Y4Sum_XY2Sum_X2Y2Sum_XY3Sum_X3Y2Sum_X2Y3Sum_XY4Sum_X4Y2Sum_X3Y3Sum_X2Y4Sum_XY5Sum_Y3Sum_XY3Sum_Y4Sum_X2Y3Sum_XY4Sum_Y5Sum_X3Y3Sum_X2Y4Sum_XY5Sum_Y6Matrix B3(O42:O51)Sum_ZSum_XZSum_YZSum_X2ZSum_XYZSum_Y2ZSum_X3ZSum_X2YZSum_XY2ZSum_Y3ZInverse A3(D54:M63)=MINVERSE(D42:M51)3rd Order Coeff.(D66:D75)=MMULT(D54:M63,O42:O51)t3_c0t3_c1t3_c2t3_c3t3_c4t3_c5t3_c6t3_c7t3_c8t3_c9 z(x,y)=c0+c1*x+c2*y+c3*x^2+c4*x*y+c5*y^2+c6*x^3+c7*x^2*y+c8*x*y^2+c9*y^3Python Program OutlineStep 1: Initialize constants:import sys# initialize constantstsOrder = 1xll = 15.1yll = 3.8xur = 42.0yur = 29.9dl = " "eoln = "\n"# subdirectory path for outputpath = "C:/temp/"# file name for outputofn = "trend_grd.txt"--------------------------------------------------------------------Step 2: Initialize variables and open output file in write mode--------------------------------------------------------------------# get trend surface ordertsOrder = input(“Enter trend surface order (1,2,or 3):”)# open trend output file in write-modeof = open(path+ofn,"w")# initialize variables# xr = x distance range from lower left to upper right of gridxr = xur - xll# YR = y distance range from lower left to upper right of gridyr = yur - yll# cellsize = distance between adjacent columns and rows in gridcellsize = 0.25ncols = int(round(xr/cellsize))nrows = int(round(yr/cellsize))nodatavalue = 1e+30--------------------------------------------------------------------Step 3: Write header information--------------------------------------------------------------------# write header infoof.write("ncols" + dl + str(ncols+1) + eoln)of.write("nrows" + dl + str(nrows+1) + eoln)of.write("xllcorner" + dl + str(xll) + eoln)of.write("yllcorner" + dl + str(yll) + eoln)of.write("cellsize" + dl + str(cellsize) + eoln)of.write("nodata_value" + dl + str(nodatavalue) + eoln)--------------------------------------------------------------------Step 4: Use a “while” loop to cycle through matrix (x,y) coordinates and calculate the trend surface value (z) at each matrix node. Values are written to file.--------------------------------------------------------------------linect = 6i = 0while i <= nrows: y = yll + (nrows - i) * cellsize j = 0 while j <= ncols: x = xll + (j * cellsize) if tsOrder == 1: z = Trend1(x,y) elif tsOrder == 2: z = Trend2(x,y) elif tsOrder == 3: z = Trend3(x,y) of.write(str(z)+eoln) linect = linect + 1 j = j + 1 i = i + 1print "{0} order trend surface matrix written to: ".format(tsOrder)+ path + ofn# close the output fileof.close()--------------------------------------------------------------------Step 5: Add the Trend1(x,y), Trend2(x,y), and Trend3(x,y) functions. Note that the function definitions should be added to the “top” of the file so they are “seen” before they are called by the main body of the program. --------------------------------------------------------------------# Function definitions# Trend1 function: returns the 1st order trend surface z valuedef Trend1(x,y): c0 = {initialize with results from spreadsheet} c1 = c2 = return c0 + x * c1 + y * c2# Trend2 function: returns the 2nd order trend surface z valuedef Trend2(x,y): c0 = {initialize with results from spreadsheet} c1 = c2 = c3 = c4 = c5 = return c0 + c1 * x + c2 * y + c3 * x**2 + c4 * x * y + c5 * y**2# Trend3 function: returns the 3rd order trend surface z valuedef Trend3(x,y): c0 = {initialize with results from spreadsheet} c1 = c2 = c3 = c4 = c5 = c6 = c7 = c8 = c9 = return c0 + x * c1 + y * c2 + x**2 * c3 + x * y * c4 + y**2 * c5 + x**3 * c6 + x**2 * y * c7 + x * y**2 * c8 + y**3 * c9Using the Trend Surface OutputOnce the program correctly computes the 1, 2 and 3 order trend surface, do the following:Step 1: import each surface into a ArcMap project. Use the “Conversion Tools > To Raster > ASCII to Raster” toolbox function. Name each surface “trend1”, “trend2”, and “trend3” Step 2: Using “raster math” subtract the 1st order surface from the 3rd order surface. Contour this residual surface at a 200 foot interval. Label the contours.Step 3: Use the “Spatial Analyst” > “Reclass” > “Reclassify” toolbox procedure to convert the residuals raster to an integer raster index. Accept the default intervals but edit the breaks so that the jump from interval 5 to 6 occurs at a zero value. This means that the interval that begins at 0, and all higher index values would represent positive residuals (i.e. the 3rd order surface was structurally higher than the 1st order surface). Step 4: Convert the reclassified integer raster into a polygon shape file using the toolbox function “Conversion Tools” > “From Raster” > “Raster to Polygon”. Make sure that the resulting polygon shape file is saved back to your working directory.Step 5: Use the toolbox “Analysis” > “Extract” > “Select” tool to select all of the positive residual polygons into a new polygon shape file named “Pos_Resid”. Remember that the index range selected should be >= to index that began at a 0 residual value (with the default number of intervals this was 6 and higher when I ran the “reclassify” tool). Step 6: Use the toolbox “Conversion Tools” > “To Geodatabase” > “Feature Class to Geodatabase” to convert the “Pos_Resid” polygon shape file to a geodatabase feature class. Remember that you have to first create a geodatabase target file with ArcCatalog. Name this file “trend.mdb”. After sending the positive residuals to the geodatabase you will have the shape area as a field. Use ArcMap to add the floating point fields “Area_sq_ft” , “Volume_cf”, and “Value”. Use calculation queries to fill in the fields with:[Area_sq_ft] = [Shape_area] * 4e6 {each inch on the map is 2000 feet}[Volume_cf] = [Area_sq_ft] *( ([gridcode]-6) * 200 + 100) {calculates volume of formation} [Value_] = [Volume_cf] * (Price_per_cf) * 0.12 {Google the natural gas market value; assume 12% porosity}Step 7: Use the table statistics function to report:Total area of positive residuals (square feet):__________________Total volume of positive residuals (cubic feet):___________________Total value of natural gas contained in formation: __________________Step 8: Produce a hard copy map that displays the residual map as a color zone + labeled contour map as depicted in Figure 1. Use a 1:4 scale for output to letter-size landscape layout. Make sure to label the map with your name in the lower right corner.05753735Figure 1: Appearance of final residual map.Figure 1: Appearance of final residual map.centercenter00 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download