Python Programming Project

Python Programming Project

Introduction

In 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:

1. 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 result
2. Create 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 Calculator

Step 1: import the “arcpy”, “math”, and “sys” libraries

Step 2: Convert the command line strings to numeric values

# command line strings are extracted as below

X1param = float(sys.argv[1])

# note the “float()” casting converts the string to a floating point number

# do the same for y2param, x2param, y2param

Note 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/pi

Else

azimuth = 90 - atan2(dy,dx) * 180/pi

Step 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 Calculator

The 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 asstrvar.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.457

Trend Surface Programming

The trend surface project will consist of 2 major components:

1. 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.
2. 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 480

nrows 450

xllcorner 378923

yllcorner 4072345

cellsize 30

nodata_value -32768

43 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 Coefficients

The 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)
N / Sum_X / Sum_Y
Sum_X / Sum_X2 / Sum_XY
Sum_Y / Sum_XY / Sum_Y2
Matrix B1 / (H3:H5)
Sum_Z
Sum_XZ
Sum_YZ
Inverse A1 / (D8:F10) / =MINVERSE(D3:F5)
1st Order Coeff. / (D13:D15) / =MMULT(D8:F10,H3:H5)
t1_c0
t1_c1
t1_c2
z(x,y) =c0 + x*c1 + y*c2

Note 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)
N / Sum_X / Sum_Y / Sum_X2 / Sum_XY / Sum_Y2
Sum_X / Sum_X2 / Sum_XY / Sum_X3 / Sum_X2Y / Sum_XY2
Sum_Y / Sum_XY / Sum_Y2 / Sum_X2Y / Sum_XY2 / Sum_Y3
Sum_X2 / Sum_X3 / Sum_X2Y / Sum_X4 / SumX3Y / Sum_X2Y2
Sum_XY / Sum_X2Y / Sum_XY2 / Sum_X3Y / Sum_X2Y2 / Sum_XY3
Sum_Y2 / Sum_XY2 / Sum_Y3 / Sum_X2Y2 / Sum_XY3 / Sum_Y4
Matrix B2 / (K18:K23)
Sum_Z
Sum_XZ
Sum_YZ
Sum_X^2Z
Sum_XYZ
Sum_Y^2Z
Inverse A2 / (D26:I31) / =MINVERSE(D18:I23)
2st Order Coeff. / (D13:D15) / =MMULT(D26:I31,K18:K23)
t2_c0
t2_c1
t2_c2
t2_c3
t2_c4
t2_c5
Z(x,y) = c0 + x*c1 + y*c2 + x^2*c3 + x*y*c4 + y^2*c5

For the 3rd order trend surface:

Matrix A3 / (D42:M51)
N / Sum_X / Sum_Y / Sum_X2 / Sum_XY / Sum_Y2 / Sum_X3 / Sum_X2Y / Sum_XY2 / Sum_Y3
Sum_X / Sum_X2 / Sum_XY / Sum_X3 / Sum_X2Y / Sum_XY2 / Sum_X4 / Sum_X3Y / Sum_X2Y2 / Sum_XY3
Sum_Y / Sum_XY / Sum_Y2 / Sum_X2Y / Sum_XY2 / Sum_Y3 / Sum_X3Y / Sum_X2Y2 / Sum_XY3 / Sum_Y4
Sum_X2 / Sum_X3 / Sum_X2Y / Sum_X4 / Sum_X3Y / Sum_X2Y2 / Sum_X5 / Sum_X4Y / Sum_X3Y2 / Sum_X2Y3
Sum_XY / Sum_X2Y / Sum_XY2 / Sum_X3Y / Sum_X2Y2 / Sum_XY3 / Sum_X4Y / Sum_X3Y2 / Sum_X2Y3 / Sum_XY4
Sum_Y2 / Sum_XY2 / Sum_Y3 / Sum_X2Y2 / Sum_XY3 / Sum_Y4 / Sum_X3Y2 / Sum_X2Y3 / Sum_XY4 / Sum_Y5
Sum_X3 / Sum_X4 / Sum_X3Y / Sum_X5 / Sum_X4Y / Sum_X3Y2 / Sum_X6 / Sum_X5Y / Sum_X4Y2 / Sum_X3Y3
Sum_X2Y / Sum_X3Y / Sum_X2Y2 / Sum_X4Y / Sum_X3Y2 / Sum_X2Y3 / Sum_X5Y / Sum_X4Y2 / Sum_X3Y3 / Sum_X2Y4
Sum_XY2 / Sum_X2Y2 / Sum_XY3 / Sum_X3Y2 / Sum_X2Y3 / Sum_XY4 / Sum_X4Y2 / Sum_X3Y3 / Sum_X2Y4 / Sum_XY5
Sum_Y3 / Sum_XY3 / Sum_Y4 / Sum_X2Y3 / Sum_XY4 / Sum_Y5 / Sum_X3Y3 / Sum_X2Y4 / Sum_XY5 / Sum_Y6
Matrix B3 / (O42:O51)
Sum_Z
Sum_XZ
Sum_YZ
Sum_X2Z
Sum_XYZ
Sum_Y2Z
Sum_X3Z
Sum_X2YZ
Sum_XY2Z
Sum_Y3Z
Inverse A3 / (D54:M63) / =MINVERSE(D42:M51)
3rd Order Coeff. / (D66:D75) / =MMULT(D54:M63,O42:O51)
t3_c0
t3_c1
t3_c2
t3_c3
t3_c4
t3_c5
t3_c6
t3_c7
t3_c8
t3_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^3

Python Program Outline

Step 1: Initialize constants:

import sys

# initialize constants

tsOrder = 1

xll = 15.1

yll = 3.8

xur = 42.0

yur = 29.9

dl = " "

eoln = "\n"

# subdirectory path for output

path = "C:/temp/"

# file name for output

ofn = "trend_grd.txt"

------

Step 2: Initialize variables and open output file in write mode

------

# get trend surface order

tsOrder = input(“Enter trend surface order (1,2,or 3):”)

# open trend output file in write-mode

of = open(path+ofn,"w")

# initialize variables

# xr = x distance range from lower left to upper right of grid

xr = xur - xll

# YR = y distance range from lower left to upper right of grid

yr = yur - yll

# cellsize = distance between adjacent columns and rows in grid

cellsize = 0.25

ncols = int(round(xr/cellsize))

nrows = int(round(yr/cellsize))

nodatavalue = 1e+30

------

Step 3: Write header information

------

# write header info

of.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 = 6

i = 0

while 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 + 1

print "{0} order trend surface matrix written to: ".format(tsOrder)+path + ofn

# close the output file

of.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 value

def 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 value

def 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 value

def 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 * c9

Using the Trend Surface Output

Once 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:

1. [Area_sq_ft] = [Shape_area] * 4e6 {each inch on the map is 2000 feet}
2. [Volume_cf] = [Area_sq_ft] *(([gridcode]-6) * 200 + 100) {calculates volume of formation}
3. [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:

1. Total area of positive residuals (square feet):______
2. Total volume of positive residuals (cubic feet):______
3. 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.

1