This is the sixth blog in a series on Lidar Solutions in ArcGIS. The author of this blog is Clayton Crawford, lead Product Engineer on ESRI's 3D Analyst team in the Software Products group in Redlands.

Updating a portion of a terrain dataset with new measurements

The ability to update a surface is important to people responsible for providing accurate, up to date surfaces and people performing analysis on those surfaces. Updates come in different forms:

• Adding ancillary data (e.g., breaklines)
• Removing or replacing bad data
• Using newer or more accurate data
• Increasing extent with additional data
• Inserting design/modeled data to perform ‘what-if’ analysis

These kinds of updates are best performed on the measurements used to construct a surface rather than on derivatives like raster DEMs. Those can be recreated as needed after the measurement edits have taken place. Terrain datasets support this editing model because they maintain a direct link to the source measurement data. When you modify the measurements, you are automatically modifying the terrain in the same process.

How terrains are edited

Editing a terrain dataset is really about editing measurements. Using standard feature edit tools, you can manipulate the measurements that reside in feature classes that participate in a terrain.

Terrains are made from one or more feature classes with some simple rules for each that control how they get used to shape the terrain surface. For example, a multipoint feature class containing lidar points can get added as mass points, a line feature class containing streams and lake shores is used as a source of breaklines, and a polygon feature class can control the data area boundary.

Most feature classes used to define a terrain are what we call referenced. This means that the terrain maintains a pointer, or handle, to them. The terrain prevents its referenced feature classes from being deleted, and pays attention to any edits that occur on them including the addition, deletion, or modification of feature geometry. You can use the feature editor in ArcMap as well as geoprocessing tools to modify these feature classes. A terrain will automatically flag itself as ‘dirty’ in areas where edits were made. Then the terrain can be rebuilt to bring its pyramid in sync with the updated features. It does this based on the dirty areas so it’s a local, or partial, process; the entire terrain does not need to be reconstructed.

Multipoint feature classes have the option of being embedded. When multipoints are embedded the terrain build process copies the points into pyramid tables held private by the terrain and it becomes the container for the points. The terrain does not reference the source feature class. That can be deleted, allowing you to retrieve what is typically a substantial amount of disk space; approximately 1GB per 150 million points. Terrain specific tools, Add Terrain Points (which can both add and replace) and Remove Terrain Points, are used to edit the embedded points based on an area of interest. These tools also offer the benefit of being BLOB attribute aware so if you have any LAS attributes stored with those multipoints (for more on this topic see blog on creating intensity images) the tools know how to keep the BLOB based values in sync with the points relative to the edits. For example, if a few vertices of a multipoint are deleted from an embedded feature class, the terrain will delete the corresponding BLOB based attribute values for those points.

Appending measurements

Measurements can be added to a terrain via the Append and Add Terrain Points geoprocessing tools. The Append tool operates on regular (referenced) feature classes. The Add Terrain Points tool is used to add or replace points in embedded feature classes. You can also add a feature class to an existing terrain via the Add Feature Class To Terrain tool but be aware that this is treated as a schema edit that invalidates the entire terrain, requiring a full rebuild. If data is to be added incrementally, it’s best to append it to a feature class that already participates in the terrain than add a new feature class to the terrain for each new set of data.

Let’s take a scenario where data is provided in phases. In this example the bare earth lidar points are made available first. The breaklines come later in several deliveries. Knowing this schedule, you can create a terrain referencing the lidar multipoint feature class plus an empty line feature class held in anticipation of the breaklines. See Figure 1 with a zoomed-in view of a terrain made solely with bare earth lidar points.

When some of the breaklines are made available they are added to the terrain by adding them to the line feature class referenced by the terrain. This is done using the geoprocessing Append tool (Figure 2).

After running Append, the terrain will become ‘dirty’ in the areas where the lines were added. To see the dirty areas you add a dirty area renderer from the Symbology tab of the terrain layer Properties dialog (Figures 3a and 3b).

At this point the terrain needs to be rebuilt. This is done using either the Build Terrain geoprocessing tool or the Build button on the Update tab of the terrain Properties dialog in ArcCatalog. Once the terrain is re-built the improvement made by the breaklines is evident (Figure 4).

Replace

With lines and polygons in referenced feature classes the replacement of measurements is a two step process. First you delete the old then append the new. If you’re only dealing with a handful of features then you can select and delete them using the Editor in ArcMap. For larger collections rely on geoprocessing tools. For example, use Select By Location followed by Delete Features and Append.

It’s easiest to replace lidar points if they’re embedded. There’s a geoprocessing tool called Add Terrain Points that has a Replace option. This will replace all the points within a given area. So, if you discover something was wrong with a few source point files that were used to build a terrain you can replace them without needing to rebuild the entire terrain from scratch. Figure 5 shows an example where one area, in an otherwise bare earth model, that was inadvertently loaded with first return data.

To fix this problem load the replacement data into a new multipoint feature class. Then run the Add Terrain Points tool with the Replace option. By default, the replacement area will come from the extent of the input feature class (Figure 6).

Once the points have been replaced the terrain needs to be rebuilt to update the affected area. Run the Build Terrain geoprocessing tool or use the Build Terrain button on the Update tab of the terrain Properties dialog in ArcCatalog – either one fixes the terrain.

Conclusion

There are times when updates to a surface model are needed, be it for quality improvement or what-if scenario analysis. It’s hard to make these types of updates to derivative products like raster DEMs without ending up with some anomalies around the update area. It’s more appropriate to modify the source measurement data from which the surface model is derived. For larger datasets, like those coming from lidar, it’s also desirable that datasets only be reprocessed where the updates occur rather than rebuilding everything. Terrain datasets support this by maintaining links to their source measurements in the geodatabase and their use of dirty areas.

That concludes part 6 of Lidar Solutions in ArcGIS. Subscribe to this blog or check back in a month or so for a discussion on how to minimize noise from lidar for contouring and slope analysis.

This blog post was written by David Wynne, Product Engineer on the Geoprocessing Team in ESRI’s Software Products group in Redlands.

The Python Community announced the release of Python 2.5.4 late in December of 2008. Python 2.5.4 is a bug fixes release.

While ArcGIS 9.3.1 was shipped with Python 2.5.1, Python 2.5.4 has now been certified and is supported. In limited testing with ArcGIS 9.3, no issues were encountered.  To take advantage of fixes in Python 2.5.4, it can simply be installed on top of Python 2.5.1. That means, it is not necessary to uninstall Python 2.5.1 first.* 

To see and to download the release notes go here.

* If you uninstalled Python 2.5.1 and you work with the Spatial Statistics geoprocessing tools, you must reinstall Numerical Python (NumPy) from the ArcGIS installation DVDs.

For more information, see FAQ: Can Python 2.5.4 be used with ArcGIS?

This blog post was written by Lauren Scott, Geoprocessing/Spatial Statistics Product Engineer in the Software Products Group at ESRI in Redlands.

There are a number of good resources where you can learn about the new regression analysis tools in ArcGIS 9.3.  The newest resource is a free one-hour Web seminar,  Regression Analysis Basics in ArcGIS 9.3, available for download from the ESRI Virtual Campus.

 Other Regression Analysis resources include:

• Regression Analysis Basics Documentation
• Ordinary Least Squares Regression Tool Help and Learn more about how OLS works
• Geographically Weighted Regression Tool Help and Learn more about how GWR works
• Interpreting OLS tool results
• Interpreting GWR tool results
• Several articles can be found in the Spring 2009 issue of ArcUser magazine [Try the Flip Page view of ArcUser!]:
• Answering Why Questions: An introduction to using regression analysis with spatial data;
• Regression Analysis Components;     and
• Applying Geographically Weighted Regression to a Real Estate Problem.
• GEOconnexion International Magazine, December/January 2009, includes an article titled Regression Analysis Tools for GIS Modeling
• Chapter 5, Analyzing Geographic Relationships, in The ESRI Guide to GIS Analysis, Volume 2 discusses regression analysis and some of the challenges when it is applied to spatial data.

Enjoy!

You may not have noticed, but you can send us feedback about this site. There's a feedback link at the bottom right of the home page, as illustrated at the bottom of this post.

We're particularly interested in hearing from you about blog post topics you'd like to see and suggestions on how we could make this site more useful.

With the Calculate Field tool, you can easily create expressions that use some property of a feature's shape, such as length or area. You can also convert length and area to different units. The illustration below shows calculating the MILES field to the length (in miles) of each line feature.

Details about the Expression syntax

The expression syntax -- !shape.length@miles! -- starts and ends with the special 'decode field' character. This decode field character is different depending on the value of the Expression Type parameter. For Expression Type of PYTHON or PYTHON_9.3, the character is an exclamation point ("!"), as shown in the illustration above. For VB, the decode field characters are "[" and "]".

Inside the decode field characters, you can put the name of any field found on the input table. You can also use the reserved field name "shape" followed by a period and an attribute of the shape. In this case, we're using the length attribute: !shape.length!

Finally, we convert shape.length to miles with the use of the special unit conversion syntax '@<unit>'. Unit can be any of the following: CENTIMETERS | DECIMALDEGREES | DECIMETERS | FEET | INCHES | KILOMETERS | METERS | MILES | MILLIMETERS | NAUTICALMILES | POINTS | YARDS.

If your shape type is polygon, you can use these areal unit keywords: ACRES | ARES | HECTARES | SQUARECENTIMETERS | SQUAREDECIMETERS | SQUAREINCHES | SQUAREFEET | SQUAREKILOMETERS | SQUAREMETERS | SQUAREMILES | SQUAREMILLIMETERS | SQUAREYARDS

More things you can do with the shape object and its attributes

The shape field is actually a geometry object with the following properties.

If, for example, you wanted the X-coordinate of the first point of a line, you would use !shape.firstpoint.x! To find the minimum X-coordinate of the line, use !shape.extent.xmin!

Simple equations

The expression can contain simple equations. For example, to calculate minutes it would take to traverse a street segment going 30 miles per hour, use:

!shape.length@miles! / (30.0 / 60.0)

 Note that:

!shape.length@miles! / (30 / 60)

won't work - when doing real arithmetic you need the decimal point after 30 and 60.

Code blocks

See http://www.esri.com/news/arcuser/0507/files/pythonscript.pdf for information on using code blocks in Calculate Field

This blogs is written by Clayton Crawford, a Product Engineer in the Software Products Group's 3D Team in Redlands. 

This is the third blog in a series on Lidar Solutions in ArcGIS. Links to the first and second posts are below. 

Data Area Delineation from Lidar Points

It's common for lidar or photogrammetric data for a survey to be delivered without a detailed data area boundary. Often, the xy extents of the survey are defined by a tile system that covers an area of interest and the data fills these tiles (Figure 1) or the data are simply guaranteed to cover some minimal extent and there is no explicit or absolute boundary other than what can be inferred (Figure 2). Either way, the area of coverage is usually not a cleanly filled rectangle.

The problem

If a surface is made without declaring the data area up front (i.e., by including a clip polygon when defining a terrain dataset or TIN) then some of what are actually voids around the perimeter get treated as data areas. Analytic results in these areas are unreliable because height estimates are based on samples that can be far away.

Take, for example, the data shown in Figure 3. The graphic on the left depicts a dense collection of lidar points shown in green. The gaps in the interior are water bodies (where lidar is typically omitted). The irregularly shaped data boundary is easy to see but unless an explicit extent is provided, in the form of a clip polygon, TIN and terrain related tools will fill in voids, greatly oversimplifying the actual data extent.

We know areas outside the data collection extent should be excluded from the surface. The problem is coming up with the polygon that provides an accurate representation of this extent. Hand digitizing is laborious. There’s an easier way.

The solution

What we need to do is synthesize a data boundary from the points that can be used to enforce a proper interpolation zone in the surface (Figure 4).

The key to the equation is the average point spacing. This is the primary variable to use when going after the data area. Surveys usually have explicit minimums on point spacing in order to provide control for interpolators. Areas that don't meet density requirements are exceptions. They usually fall in one of the following categories: water bodies, obscured areas, and ‘holidays’ (send the latter back to the data provider for repair). The vast majority of the data will meet sample density specifications. Point spacing is usually reported in metadata. If you don’t know it try the 3D Analyst Point File Information tool. Alternately, eyeball it using the Measure tool in ArcMap. To learn more about point spacing use the link at the top.

Once you know the point spacing you follow these steps:

1. Rasterize the points with Point To Raster.
2. Assign one value to all data cells with Con.
3. Fill small NoData areas with Expand.
4. Reduce the overall extent of data cells with Shrink.
5. Vectorize the raster with Raster To Polygon.
6. Use a VBA script to remove small holes (interior rings) from polygon(s).

Note that a Spatial Analyst license is required to run a number of these tools.

The remainder of this document goes into more detail on each of these steps.

Point To Raster

Rasterization of the points will help aggregate the area covered by the points and provides a good data structure to work with for subsequent steps. You just need to tell the tool what cell assignment type to use and the output cellsize. Use the COUNT as the assignment type. As far as cellsize goes, specify a value that's several times larger than average point spacing. Otherwise, you’ll get a lot of noise because the points aren’t evenly spaced. From the standpoint of processing efficiency and noise reduction, the larger cellsize you use the better, but there will be a tradeoff with the tightness of fit in the end result. A good place to start is 4 times the average point spacing (Figure 5).

Con

The use of the Con tool in this workflow is to simply turn any and all data cells into cells with one value. This value defines a raster ‘zone’ that will be expanded in the next step. All that’s needed is to take the output from Point To Raster and provide a constant value for a positive expression. All non-zero value cells will be considered positive and be assigned the constant. Since COUNT was used as the cell assignment method during rasterization, any cell with a point in it must have a value greater than zero (Figure 6).

Expand

Unless you used a very coarse cellsize relative to your average point spacing, there's a likelihood many NoData cells remain. Most of these can be eliminated using Expand (Figure 7). You want to remove most of these so the polygon produced during vectorization in a later step isn’t full of a million holes. That would be unnecessarily expensive.

 Expand will grow the zone of interest outward. In our case the zone is all the data cells, coded with a value of 1. This effectively eliminates small gaps found in the interior.

It’s okay if some isolated NoData areas remain in the output. You just don’t want thousands upon thousands. The remainder will be handled in the last step.

Shrink

While Expand eliminates isolated NoData cells it also grows the data area outward. It needs to actually be brought in a little. Clip polygons need to be smaller than the actual point extent, so when terrain or TIN tries to estimate z values along the polygon boundary points can be found on both sides. This is needed to get good z estimates. To reduce the raster’s data boundary use Shrink (Figure 9).

At this point you should have a relatively clean raster with the extent of its data cells slightly within the lidar point extent.

Raster To Polygon

The Raster to Vector tool will convert the raster to a polygon feature class. Make sure the Simplify polygons option is checked. If it’s not, the output will be stair-stepped, rather than smooth, and contain more vertices than necessary (Figure 10).

At this point, the process is almost complete. You need to review the output for correctness. Chances are there’s one more thing that needs to be done: the removal of remaining holes inside your clip polygon.

Remove Interior Rings from Polygons VBA Script

If holes are present in the resulting polygon, grab the RemoveInteriorRings VB script from ArcScripts online (http://arcscripts.esri.com/details.asp?dbid=16019). This will edit out the internal rings leaving just the exterior boundaries. To run the script, follow these steps:

1. Have your polygon feature class output from Raster To Polygon as the first layer in a map document.
2. Go into VBA inside ArcMap via Tools>Macros>Visual Basic Editor.
3. From the Project window right click on the project, select Import File, and bring in the VBA module you downloaded from ArcScripts.
4. Run the RemoveInteriorRings macro.

You should now have a clip polygon that can be used when defining a terrain dataset or TIN. It should conform to the data extent of the points but fall slightly inside them (Figure 11).

That concludes this part of Lidar Solutions in ArcGIS. Subscribe to this blog or check back in a month or so for a discussion on the estimation of forest canopy density and height from lidar.

The first post can be found here: http://blogs.esri.com/Dev/blogs/geoprocessing/archive/2008/10/29/Lidar-Solutions-in-ArcGIS_5F00_part-1_3A00_-Assessing-Lidar-Coverage-and-Sample-Density.aspx 

The second post can be found here: http://blogs.esri.com/Dev/blogs/geoprocessing/archive/2008/12/15/Lidar-Solutions-in-ArcGIS_5F00_part2_3A00_-Creating-raster-DEMs-and-DSMs-from-large-lidar-point-collections.aspx 

I recently had a project that involved creating lots of Python script tools. Since other people were going to use these tools, I wanted my error handling to be robust and informative. This post is about some tips and tricks I discovered during this project.

Error handling basics

The basics of handling Python errors are covered in the 9.3 help topic Error handling with Python. The Python.org document Errors and Exceptions has more detailed information.

Tip #1 - Use the arcgisscripting.ExecuteError exception

In version 9.2, we introduced the arcgisscripting object along with a new exception class, arcgisscripting.ExecuteError. (This arcgisscripting.ExecuteError exception class wasn't documented in version 9.2, so few people knew about it.)  The arcgisscripting.ExecuteError exception is thrown whenever a geoprocessing tool or geoprocessing function encounters an error. What this means is that you can divide errors into two groups, geoprocessing errors (those that throw the arcgisscripting.ExecuteError exception) and everything else.  You can then handle the errors differently, as demonstrated in the code below:

import arcgisscripting
gp = arcgisscripting.create(9.3)
try:
    result = gp.getcount("C:/blah.shp")
    # x = y
        
# Return GEOPROCESSING specific errors
#
except arcgisscripting.ExecuteError:
    gp.AddError(gp.GetMessages(2))
    
# Return any PYTHON or system specific errors
#
except:
    gp.AddError("Python or system error occurred")

The code above is used as the source for a script tool. To keep things simple, the script tool has no parameters. When the script tool is executed, the call to getcount produces an error because the dataset "C:/blah.shp" doesn't exist. As a result, the progress dialog looks as follows:

To prove how geoprocessing errors and Python errors are handled differently, change the two lines of code as follows:

    # result = gp.getcount("C:/blah.shp")
    x = y

Now when the script is executed, an error will occur because the variable 'y' is undefined, and the progress dialog will display as follows:

Tip #2 - Beware of getting error messages from a result object

Before moving on, a quick word about the result object, shown below:

    result = gp.getcount("C:/blah.shp")

If the call to getcount above raises an exception, the result object is null. This means you cannot retrieve error messages from the result object. For example:

import arcgisscripting
gp = arcgisscripting.create(9.3)
try:
    result = gp.getcount("C:/blah.shp")
        
# Return GEOPROCESSING specific errors
# (this method is INCORRECT!)
except arcgisscripting.ExecuteError:
    gp.AddError(result.GetMessages(2))
    
# Return any PYTHON or system specific errors
#
except:
    gp.AddError("Python or system error occurred")

The above code will fail with the message "name 'result' is not defined". This is due to the fact that the result object is null. Since the result object is null, a Python error will be raised since a null object doesn't have a GetMessages() method. Note -- a result object created by calling a geoprocessing service on an ArcGIS Server is never null. Null result objects a created only when a tool is run locally and it raises an error.  For more information about using the result object, see Getting results from a geoprocessing tool.

Tip #3 - Use AddReturnMessage() to preserve links to error codes

import arcgisscripting
gp = arcgisscripting.create(9.3)
try:
    result = gp.getcount("C:\\blah.shp")
        
# Return GEOPROCESSING specific errors
#
except arcgisscripting.ExecuteError:
    for msg in range(0, gp.MessageCount):
        if gp.GetSeverity(msg) == 2:
            gp.AddReturnMessage(msg)
    
# Return any PYTHON or system specific errors
#
except:
    gp.AddError("Python or system error occurred")

Now when the script is executed, the progress dialog will look as follows:

Tip #4 -- Use traceback to return more information about the error

The help system topic Error handling with Python shows you how to use Python's traceback to retrieve the line number causing the error, the name of the .py file, and, for Python errors, a description of the error.

At version 9.3, you can run script tools in-process, as described in Running a script in process. Scripts run in-process execute much faster than scripts run out-of-process, so you always want to run the script in-process. When a script is run in-process, ArcGIS reads the .py file into memory and then executes it. Because ArcGIS puts the .py files into memory, the Python traceback module doesn't know the name of the .py file and will report the file name as "<string>". If you want to report the name of the .py file causing the error (very useful for debugging), you'll have to provide it yourself.

The script below has a local routine, trace(), that uses traceback to return the line number, the file name, and the error description. The routine also contains the the hard-coded file name (shown in bold) which you must change for your script.

import arcgisscripting
gp = arcgisscripting.create(9.3)
try:
    def trace():
        import os, sys, traceback
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]  # script name + line number
        line = tbinfo.split(", ")[1]
        
        # Construct the pathname to this script.  Get the pathname
        #  to the folder containing the script using sys.path[0]. 
        #  Modify line below and replace 'myscript.py' with the name of
        #  this script. 
        #
        filename = sys.path[0] + os.sep + "myscript.py"
        
        # Get Python syntax error
        #
        synerror = traceback.format_exc().splitlines()[-1]
        return line, filename, synerror
    if __name__ == '__main__':  
        # y = x
        result = gp.getcount("C:/blah.shp")
        
except arcgisscripting.ExecuteError:
    # Call our trace routine to retrieve line number and filename.
    #  (returned variable 'err' is ignored since this is a gp error.
    #
    line, filename, err = trace()
    gp.AddError("Geoprocessing error on " + line + " of " + filename + " :")
    for msg in range(0, gp.MessageCount):
        if gp.GetSeverity(msg) == 2:
            gp.AddReturnMessage(msg)
    
except:
    line, filename, err = trace()
    gp.AddError("Python error on " + line + " of " + filename)
    gp.AddError(err)

Below are illustrations of the progress dialog showing the line number and file name for both a geoprocessing error and a Python syntax error (by commenting out the "y = x" statement).

Tip #5 -- Use Finally statement to clean up cursors

Something I discovered by reading the Python doc Errors and Exceptions is the 'finally' statement. No matter what happens in your Python script, code in a 'finally:' block always gets executed. It's a great place to put clean-up actions, such as deleting cursors, as shown in the code below.

import arcgisscripting
gp = arcgisscripting.create(9.3)

def trace():
    import os, sys, traceback
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0] 
    line = tbinfo.split(", ")[1]
    filename = sys.path[0] + os.sep + "testscript.py"
    synerror = traceback.format_exc().splitlines()[-1]
    return line, filename, synerror

if __name__ == '__main__':
    # Define variables for a cursor and its row
    #
    cursor, row = None, None
    try:
        cursor = gp.searchcursor("E:/Data/CityOfSanFrancisco.gdb/FireStations")
        row = cursor.Next()
        while row:
          # Code here...
          # for demonstration...this next statement will cause an error
          y = x
          row = cursor.Next()
    except arcgisscripting.ExecuteError:
        line, filename, err = trace()
        gp.AddError("Geoprocessing error on " + line + " of " + filename + " :")
        for msg in range(0, gp.MessageCount):
            if gp.GetSeverity(msg) == 2:
                gp.AddReturnMessage(msg)
    except:
        line, filename, err = trace()
        gp.AddError("Python error on " + line + " of " + filename)
        gp.AddError(err)
    finally:
        # Close cursor and row objects.
        #
        del row, cursor
            
        # Print a message for demo purposes
        #
        gp.AddMessage("\n   ** Cursor and row have been deleted ** \n")

As illustrated below, the cursor gets deleted even if there is an error in the script.

In the main body of the script, the cursor and row object variables are declared and initialized to None (shown in bold in the above code).  This ensures that these two variables exist when the del statement in the finally clause is executed; if you try to delete a non-existent variable, an exception will be thrown.  

About cursors and the del statement

The Python del statement deletes a variable. You can use it to delete any Python variable, such as a list or dictionary. However, you rarely need to delete native Python variables.  The del statement is mainly used for deleting non-native variables, such as geoprocessing cursors and rows obtained from the cursor. When you delete a cursor and its row, all pending changes to the row (caused by an update or insert cursor) are flushed and all locks on the dataset are removed.  You should always delete cursors and rows obtained from the cursor.  When deleting a cursor, delete the row object first, then the cursor.

This posting was written by Mark Zollinger, a Product Engineer and long time UNIX user with the Geoprocessing team.

Do you like Solaris or Linux? Do you "do GIS" on one of those platforms? Maybe you've been asked to convert some ArcGIS Workstation AML applications over to ArcGIS Engine. Are you a command-line junkie? (There are more of us than you might think)

If, for one of the above or any other reason, you find yourself wondering if you should try developing GIS on Solaris or Linux, then it's time to stop wondering. Join in with others who include a flavor of *nix in their list of working development platforms.

There are several approaches to GIS development that work very well on Solaris and Linux.

If you are into low-level programming, you should consider Java or C++, both of which can access ArcGIS's ArcObjects API. I've used Java quite a bit over the years, and can tell you that the doors are wide open on what you can do with it. If this sound like the thing for you, then you can move along now, there's nothing to see here. This post isn't about ArcObjects at all.

If you lean towards scripting and higher level abstractions, then geoprocessing's more coarse-grained object model, then you are in the right place. The rest of this post will show you how to set up an IDE for Python development and write a very simple geoprocessing script.

I mentioned a moment ago that I've written my share of Java ArcObjects code. I've got to confess though, that I've always been a fan of scripting languages. Python and GP is where I spend most of my time these days.

Setting up

"Okay", you ask, "How do I start writing GP code in Python on non-Windows platforms?"

You can do Python and GP from the command line, and I often prefer the focus this brings. I've found, however, that most folks prefer to work in a graphic IDE (Integrated Development Environment).

Setting up an IDE that you can use

Probably your best choice is the Eclipse IDE with the Pydev plugin. Eclipse is an Open Source IDE, intended to work with a wide variety of programming languages. Pydev is the plugin that specifically enables Python development in Eclipse.

The latest version of Eclipse, named Ganymede, has not been fully tested by ESRI for use with ArcGIS 9.3, and is not available for some older UNIX releases. We recommend using version 3.3.2, which you can find in the archive at:

To install Eclipse, follow the instructions in the "Installing Eclipse" paragraph at:

(These instructions are for Ganymede, but they basically the same for version 3.3.2)

Once you have Eclipse installed and running, you need to add the Pydev plugin. Just follow the "Getting Started" link on the left side of the Pydev homepage.

Note that you need the Open Source Pydev plugin. The Pydev Extensions plugin is commercial trialware, and is optional.

When you get to the "Configuring the interpreter" page, you are told to choose the Python interpreter that Pydev should use to run Python programs. Use the Python installation you use for Engine development/testing. In most cases, this is the one that got installed with the Engine Runtime:

Note: you might get an error message while trying to choose the interpreter. If you do, click the error dialog's Details button. The message probably includes something like this:

This is an indication that the python interpreter is not correctly set up. The most likely cause is that your ArcGIS Engine environment has not yet been established. Exit Eclipse and source the init_engine shell script that applies to the shell you are using. When you start Eclipse the problem should go away.

Also, item 3 on the "Configuring the interpreter" page asks you to set the SYSTEM PYTHONPATH. When you do this, make sure to include $ARCGISHOME/bin. This matches the PYTHONPATH that gets set when you source init_engine.

Writing a Python script that does GP stuff

Now that you have a development environment, how do you use all this great GP stuff? Let's create a very simple Python script to get you started on the road to GP development. In most Python scripts, the first thing to do is to import the modules that you will need. For GP, we need to import the arcgisscripting module. Let's also import the os and system modules:

To access GP, we need to create a Geoprocessor object. This is the object through which you will call any GP functions or objects

By passing in the 9.3 argument, you are getting the newer version of the Geoprocessor object, which is a little more Python friendly. For this example we won't be using any of the Python friendly features, but it is never too early to start a good habit.

Now that you have a GP object, what do you do with it? Let's set GP's workspace property and print it out, just to prove to ourselves that it really worked. For this example, let's set the workspace to a "data" directory just underneath the current Python directory.

(Of course you can set the workspace to any accessible directory, not just the one returned by sys.path[0]. I'm just using sys.path[0] because I know it's accessible.)

sys.path is a list of strings where Python looks for modules. The first item in the list is the current Python directory. Calling os.path.join is a great way to build up filesystem paths, since it correctly handles all the path separator issues. (More on path separators later.) When you run this script (Run -> Run As -> Python Run), you should see "GP Workspace= /some/path/data" printed in the Console tab near the bottom of the Eclipse window.

You can use Eclipse or the operating system to create the data directory, since it probably doesn't exist yet.

For this example, let's suppose you have point data representing hawk sitings in a national forest. Your study area is a subset of the forest where a particular plant in known to grow. What you need to do is clip your point data using the polygon of your study area.

Once you have copied your data into the aforementioned data directory, your script could execute the clip tool like this:

gp.clip_analysis("nat_for.gdb/hawk_sightings", "study_area.shp", "nat_for.gdb/study_sightings")

Just to be sure it worked, you can print out the messages that the tool generated like this:

Notice that in calling the clip tool, we used forward slashes (/). Back slashes (\) and forward slashes work equally well in ArcGIS on ALL platforms, so you don't need to worry about switching back and forth when you change to Windows. You can even use back slashes on Linux and Solaris. The problem is that in Python, back slashes are a special character, so to use them, you have to treat them differently than most other characters.

- You can escape them: "\\some\\path"
- You can use raw strings: r"\some\path"
- You can use unicode strings: u"\some\path"

If you stick to forward slashes on all platforms, you improve the look of your code and reduce the effort it takes to type it.

So now you have a script that looks like this:

import arcgisscripting
import os, sys

gp = arcgisscripting.create(9.3)

gp.workspace = os.path.join (sys.path[0], "data")
print "GP Workspace = " + gp.workspace

gp.clip_analysis("nat_for.gdb/hawk_sightings", "study_area.shp", "nat_for.gdb/study_sightings")
print (gp.getmessages())

Run it, and it works, creating a new featureclass in the file geodatabase. Now wasn't that easy? To get more complicated, you just string together as many of those tools as you need to get the job done.

Welcome to the world of Geoprocessing with Python!

Links

Troubleshooting

Color issues

Upon using Eclipse on Solaris for the first time, you may see some color problems. For example, the text insertion caret might be white on a white background, making it invisible. Another common symptom is that you cannot see the outline of text input fields. If you run into this or something similar, check the color depth of your display. If it is set to 8 bits, change it to 24.

init_engine

You should have your Engine environment already established in the shell in which you start Eclipse. (source init_engine.sh or source init_engine.csh) If you see inexplicable error messages, this is the most likely cause.