September 3, 2010

Model to Script to Tool: the "Holy Trinity"

I have been doing some work for a junior mining company recently involving creating Volume/Elevation curves (see my last post). This work is accomplished in ArcGIS using Spatial Analyst and 3D Analyst which works very well but requires many steps. Over the course of this project, however, I have built a number of Python scripts that are now embedded as tools in a new toolbox in ArcMap. So I am going to take this opportunity to briefly discuss how easy it really is to create sophisticated Python geo-processing scripts in Arc and give a little example of how I recommend approaching this task.

First of all I would just like to say that ArcGIS is still the most powerful GIS in my opinion (yes I am still a fan of open source but you can't argue with the benefits of having a large budget to develop good software). I like Arc not so much for it's ability to perform the geo-processing functions but more for the ancillary tools that help you find the appropriate tools and perform these tasks. Arc has built in search capabilities and very good documentation to allow users with a wide variety of skill levels find the appropriate tools. On top of that, there is a very well developed modeling and scripting environment that gives the user the capability to reproduce these tasks in a relatively easy manner.

Ok, now that's been said lets get on with the main purpose of this post, to show how easy it is to create a tool in ArcGIS using the modeling and scripting environments. For this example I thought I'd use the volume elevation curve tool created for the project I mentioned earlier. I will be making reference to a number of steps that may not be clear if you haven't performed this analysis before but you can read my previous post for clarification.

The Model
So over the years I have found that the easiest way to start any Python scripting in Arc is to look at the modeling environment first and see if the task can be done in there. You may ask "why worry about a Python script if you can create a model", but if you need to add some iteration
Figure 1: Surface volume model

to the process (e.g. performing the analysis on successive elevations), for example, then you need to do this with a script. So I started by creating the model shown in Figure 1. The key here is to notice that I have created a number of "model parameters" which will show up as argument variables (sys.argv[]) when the model is exported to Python. These arguments will allow you to attach the fields from a custom tool to them giving the user the flexibility to add their own data, assign a range of elevations to the plane height variable and specify the output location for the delimited text file. You will also notice that there are two height variables (Min Plane Height and Max Plane Height) but only one is actually attached to the model. This is because we are planing to have an iterative model that starts at the bottom of the surface (e.g. 0 cu.m volume) and works it's way up to the top. This requires that we allow the user to assign a lower value and an upper value but because this model is only a single instance of the volume command, it can take only one value. Don't worry though, we will add the iterative properties when we edit the script. Once you have the model the way you like you can export it to a Python script (Model -> Export -> To Script -> Python). Before doing this though, I suggest testing the model to ensure it works to avoid spending an unreasonable amount of time on the scripting.

The Script
Now that we have a script you can open it up in IDLE or whichever Python editing environment you like. The script will look something like Figure 2 after it is first exported.

# ---------------------------------------------------------------------------
# sample_model_script.py
# Created on: Thu Sep 02 2010 04:23:36 PM
# (generated by ArcGIS/ModelBuilder)
# Usage: sample_model_script
# ---------------------------------------------------------------------------

# Import system modules

import sys, string, os, arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create()

# Check out any necessary licenses
gp.CheckOutExtension("3D")

# Load required toolboxes...
gp.AddToolbox("..Toolboxes/3D Analyst Tools.tb
x")

# Script arguments...
Output_Text_File = sys.argv[1]

Dam_Containment_Tin = sys.argv[2]

if Dam_Containment_Tin == '#':

     Dam_Containment_Tin = "defa" # provide a default value if unspecified

Min_Plane = sys.argv[3]


Max_Plane_Height = sys.argv[4]

# Local variables...

# Process: Surface Volume...
gp.SurfaceVolume_3d(Dam_Containment_Tin, Output_Text_File, "BELOW", Min_Height, "1")
Let's run through the components of this script. First of all you have to import the Python libraries necessary to run this script. In this case it is the system (sys), string, operating system (os) and ArcScripting (arcgisscripting) libraries. Right after this you need to create the geoprocessing object that will do the bull work (gp). Next you need to check out the extension(s) required for the type of geoprocessing to be done. In this case we only need 3D Analysist. Then we asign the different arguments to there own reusable variables (e.g. Output_Text_File = sys.argv[1] which assigns the name of the output file, as provided by the user, to it's own variable). In our example we do the same for the input dam containment tin, the lower plane height and the upper plane height. You may find that you need to change the order in which these arguments appear so that the tool variables can be placed in a logical order. So if you want the input TIN to be first followed by the lowest plane height then highest plane height and finally the output text file name then you need to ensure that you set the index number for the sys.argv[] object accordingly. Finally we run the Surface Volume_3D function using the geo-processing object we created earlier.

This script will do most of the work we need, however, we need to create a loop that will iterate through all the elevation changes needed to create a proper VE Curve. In most cases 9 - 12 points are adequate to create a VE Curve in Excel so we will create a loop that creates 9 points.


Height_Interval = (float(Max_Plane) - float(Min_Plane))/9
Height_Plane = float(Min_Plane)
for i in range(1,9):
    gp.SurfaceVolume_3d(
Dam_Containment_Tin, Output_Text_File, "BELOW",
       str(Height_Plane), "1")

    Height_Plane = float(Min_Plane) + (float(Height_Interval) * int(i))


You can see here how we have set up a variable f
or the plane height and calculated it by finding the elevation difference between highest and lowest values and then dividing that by 9 (number of points for the curve). Then we create another variable to represent the current plane height (initially equal to the lowest elevation in the containment TIN). Next we set up a loop that iterates through 9 times for which we calculate the volume below the reference plane and increase it's height each iteration by the Height_Interval. In the end we have a script that, given some user input for each of the variables, will calculate a volume elevation curve and append each value to an output text file.


Creating a Tool
Now that we have the script the way we want it we need to create an ArcToolbox tool. We do this by simply importing it. I recommend adding a new Toolbox first by right clicking on the root of the ArcToolbox window and selecting New Toolbox. You should make sure to give the toolbox an appropriate name so that it is obvious what it contains. Next create a new toolset by right clicking on the new toolbox and selecting New -> Toolset. Finally, add your newly created script by right clicking on the new toolset and selecting Add - > Script, after which you will navigate to the location of your script. Once the script is added you will need to assign parameters to it so that the user can update the variables when the double click on it (see Figure 2). You get to this dialog by right clicking on the script and selecting Properties then click on the Parameters tab. You can add each parameter to the script in the logical order set out in your script.

Figure 2: Add parameters to the script

Remember that the order the parameters show in this dialog is the order they are passed to your code. This means that the first parameter showing here will be accessed in your code through the sys.argv[1] system variable.

You will also need to make sure that you update all the Parameter Properties such as ensuring input and output data sets are assigned accordingly and that Dependencies are also assigned. This last one comes into play in situations such as the requirement to select from a list of fields.

Now when you double click on the script in the toolbox you should see a dialog something like that shown in Figure 3.


Figure 3: Tool dialog for our Script

Now you should have a workable tool that will produce a volume elevation curve from an existing containment TIN.  The intent of this post is not to show you how to do VE Curves but to give you a starting point for creating your own custom script based tools in ArcGIS.  Hope you find it useful.

1 comment:

greenbuild said...

Thanks for posting this methodology. I'm still fairly new to python, and I'm grateful for such clearly explained real world examples.

I wonder if this can all be accomplished as a model in ArcGIS 10 with the new iteration tools available in model builder. In particular the arcpy.IterateCount_mb

I was trying to post the graphic of the model I created, but I guess that is beyond the comment box capabilities. Here's the exported script version:


import arcpy
arcpy.CheckOutExtension("3D")
arcpy.ImportToolbox("Model Functions")

# Script arguments
Dam_Containment_TIN = arcpy.GetParameterAsText(0)
if Dam_Containment_TIN == '#' or not Dam_Containment_TIN:
Dam_Containment_TIN = "C:\\Data\\Hatfield\\surfacetin" # provide a default value if unspecified

Min__Elevation = arcpy.GetParameterAsText(1)
if Min__Elevation == '#' or not Min__Elevation:
Min__Elevation = "1" # provide a default value if unspecified

Max__Elevation = arcpy.GetParameterAsText(2)
if Max__Elevation == '#' or not Max__Elevation:
Max__Elevation = "9" # provide a default value if unspecified

Output_VE_File = arcpy.GetParameterAsText(3)

# Local variables:
Plane_Height = Max__Elevation

# Process: For Min to Max Plane Height
arcpy.IterateCount_mb(Min__Elevation, Max__Elevation, "1")

# Process: Surface Volume
arcpy.SurfaceVolume_3d(Dam_Containment_TIN, Output_VE_File, "BELOW", Plane_Height, "1", "0")