GIS Programming - Module 12 - Final Project

For my final project, I wrote a Python script that crops point shapefiles (quadrangles) of stream sediment geochemistry data to a study area, and isolates those points with a gold concentration of more than 0.02 ppm.

Result of Python Script, displayed in ArcMap

Partial results message from Python script run, showing data points with
high gold concentration (more than 0.02 ppm)

GIS Programming - Module 11 - Sharing Tools

 

The final module of GIS Programming involves making custom tools that we've created from scripts, ideally that are easy to use by other ArcGIS users.

The figure below is a view from the ArcMap display.  (Click on the figure to make it larger.)

Tool Dialog and ArcMap Results for the Random Buffer Tool

When the small scroll-shaped icon, in ArcCatalog, for Random Buffer Tool is double-clicked, it brings up the dialog box for the tool and parameters can be filled in.

As shown on the left of the figure, the result of running this tool is that the specified number of points are generated within the confines of the clip feature (in the case, the lavender polygon for the state of Nuevo Leon), and then buffer zones of the specified size are created around the points.

It is not too difficult to create a tool like this from a script and share it with other users.
A few lines with such code as GetParameter or sys.argv[] have been inserted into the script, which bring the specified parameters from the tool dialog in ArcMap into the script.  A toolbox is created and the script is converted to a script tool.

The parameters can be described and other information can be added to the tool, such as Help and sample scripts, by  right-clicking the  item in ArcCatalog then going to Description.

This is the last module in GIS Programming.  What did I like most about this class?

I thought the final project was a great idea.  I really learned a lot from having to be able to put things together on my own.  The weekly exercises were good for showing us how different processes work and giving us examples.  It was the project, though, where I really learned to put it all together, to create something completely new.  I feel like this will enable me to move forward and write more scripts for future work projects of my own. 

Besides the project, my favorite module was probably Exploring and Manipulating Spatial Data.  It was pretty challenging to understand some of the concepts, but it was very useful.  I was happy to have more time to work on it, and take my time absorbing the material.  

I also liked the two modules on Geoprocessing in ArcGIS and with Python.  They were very useful.

This is also the last module of my last GIS course at UWF!  It has been a ton of work, but I've really learned a lot and enjoyed it.  Hooray!  Thanks a million to all the teachers, TA's and my fellow students!

GIS Programming - Module 10 - Creating Custom Tools

This module examines how a stand-alone Python script that's been written for use with ArcGIS can be modified and turned into a script tool that can be accessed very conveniently from inside ArcMap. It's necessary to make just a few minor changes to the original script, and it can then obtain parameters that are specified by the ArcGIS user from the map environment.  The user does not even need to know anything about Python, and the new script tool can be incorporated into ModelBuilder and other scripts.

For our assignment, we created a new Toolbox from ArcCatalog and within that, a Script Tool that will clip shapefiles to the extent of a clip boundary shapefile.  The tool allows the user  to navigate to the data in its file location.

Figure 1 (below) shows the tool dialog that I created in ArcCatalog,
from the stand-alone Python script.

Figure 1. Clip Tool dialog box, ready for input by the user.

This dialog has spaces for the input path (with a default provided), single or multiple input features that will be clipped, the Clip feature, and the default output path.
It's also possible in the tool settings to add a description, shown on the right of this example, and Tool Help, including what the script looks like.

Figure 2. Dialog box, after being filled by the ArcMap user.

 The user fills in these blanks, as shown in Figure 2 (left) then clicks okay, and the script runs in the background, with no direct involvement by the user.

Figure 3.  Successful Results Message for the Clip Tool,
produced in the ArcMap environment

Although the user need not ever see the actual Python script, the script's progress can be monitored in the foreground of ArcMap,via results messages written into the original script, as shown in Figure 3, left.. These messages indicate the tool's progress (including any errors) to the ArcMap user.

Figure 4.  New shapefiles added to ArcCatalog and
the ArcMap display

Finally, after the tool runs, the resulting clipped features are stored in the Results folder, as seen in ArcCatalog, and can be added to the map display in ArcMap (both shown in Figure 4, left).

The light purple area is the clipping boundary shapefile, the State of Durango, Mexico,  Roads, railways, rivers and urban areas are the input files, and they are shown in the map display as colored lines, clipped to the extent of  the Durango shapefile.  

GIS Programming - Participation Assignment #2

This assignment was to find and summarize a technical/professional article relating to Python or GIS. The article I am reporting on is:

Visualizing structural geology: From Excel to Google Earth

T.G. Blenkinsop, 2012: Computers and Geosciences v. 45, pp. 52 – 56.

Link to article introduction at DeepDyve.com

This article is of particular interest to me for my work in structural geology and minerals exploration and especially for small companies like the one I work for, which may not have access to expensive GIS software like ArcGIS.  This system uses a virtual globe, specifically Google Earth Pro, and Excel enabled for macros.  (Google Earth Pro used to cost, but now it’s free.)

   The ability to visualize geologic elements in 3D is essential for structural analysis of the Earth’s crust.  An understanding of the deformational and tectonic history of regions and even continents is mad possible when the various planar and linear structures are brought together and their spatial relationships examined.   Measurements done in the field are normally entered into an Excel spreadsheet, with spatial coordinates and/or latitude and longitude, as well as measurements of the orientations of elements such as the limbs and axes of folds, lineation, cleavage and stratigraphic bedding.  The orientations are expressed in terms of the azimuth of a plane’s slope or a line’s direction of plunge, and the degree a linear feature’s plunge or the dip (slope) of a planar feature.  Relationships between various elements’ orientations can help determine how many generations of deformation an area has experienced, and in what order. 

The author has combined a library of 3D symbols such as rods and flat prisms of various shapes and colors, with an Excel workbook called 52K .  The Excel macro uses Visual Basic for Applications to convert field data measurements of the structural elements into a KML (Keyhole Markup Language) for Google Earth Pro.   Any field measurements recorded by hand in the past can be fairly easily entered into Excel spreadsheet form, incorporated into other structural data sets, and presented in a nicely intuitive form on top of a Google Earth landscape, viewed in the 3D mode.  This method also allows structural principles to be presented to students and novices of structural geology, to allow easier understanding of the 3D concepts involved. 

The authors present an example (pictured below) of how, by  measuring very small localized features in numerous locations across an area, general large trends can be represented alongside the Google Earth landscape and compared to other planar and linear features.  In the case of the author’s example, the consistent orientations of various types of structures indicate that the area was deformed in a single event.  This view also helps us to intuitively visualize structures which are not obvious from the field.

In this figure (from Blenkinsop’s 2012 article),  white planes are stratigraphic bedding and red planes are schistosity, or planar deformational parting of the rock fabric. The blue rods represent the hinge or bend of a fold.  The red rods represent the intersection between the two types of planes. 

GIS Programming - Module 9 - Working with Rasters

This week's exercise on working with rasters was very useful but also fairly straightforward and easy to understand.  We applied slope, aspect and raster algebra operations that we've learned in earlier GIS courses to Python scripting and learned how to explore, describe and manipulate rasters using Python.

Two rasters were provided: elevation and land cover.  Slope and Aspect rasters were developed from the elevation, and forested land cover was isolated from the land cover classification.  We determined ideal areas based on slope, aspect and forest cover and the result was a final suitability raster.

The final raster is shown below:
The green areas (value =  1) are considered suitable for the purpose of the analysis:
Slope = 5 - 20 °   Aspect = 150 - 270°  Land cover is forested
Boolean AND logic was used in the script:  all parameters must be true to give output value of 1 in the raster.

Figure 1. Final raster showing suitable areas of slope, aspect and landcover (green)

Below are the results for the script run, consisting of commentary at all stages of the process.

Figure 2.  Interactive Window Results for Script

Process Summary Notes

Which step did you have the most difficulty with? Describe 1) the problem you were having, and 2) the solution or correct steps to fix it. 

1.      This lab was not too difficult for me because I already had quite a lot of experience with rasters from previous UWF courses and already understood the concept of map algebra. 

2.      Creating the raster objects for max and min slope, max and min object and forested land cover required the most thought, because I’d never done scripting for this task before.  The exercise and assignment instructions helped with this, though. 

3.      Also, I had not worked with the Boolean AND operator much before, but instead had added and multiplied the raster values for suitability and least cost analysis.  It makes sense, though, and wasn’t hard to figure out. 

4.      The main obstacle I had was when I produced my final raster.  It had three colors instead of two, as the assignment example shows.  There were different colors representing 1 and 0, but also white area of NODATA. 

5.      Once I used the Identify tool to see what was what, and read through the lab, ArcGIS Help, and text more, I realized that I needed to remove the NODATA clause in the script of the Reclassify line in the land cover raster. 

6.      I was confused for a while with this, because the reclassed land cover raster has a value of 1 for all forested areas (41, 42, 43) but the other classes retain their original values.  It seemed necessary to assign the other classes to NODATA.

7.      After some further reading, however, I realized that this doesn’t matter, because once the Boolean operation of the 5 parameters is carried out, any value NOT equal to 1, no matter what it is, will result in a value of 0 in the final raster. 

8.      The biggest problem with this script had nothing to do with the raster management but rather with the overwriting of the file geodatabase in repeated runs of the script.  Although I wrote overwrite output, plus uncompressed  and deleted the already-existing file geodatabase in the script, there is still an error which prevents the file geodatabase from being recreated in repeated runs of the script, UNLESS Python is completely closed out between runs. 

9.      I would like to learn how to close down Python from within the script, at the end.  I tried to do this but it didn’t work. 

10.   It is not necessary, however, to delete the .gdb each time from the catalog.

GIS Programming - Module 8 - Working with Geometries

This week's Module covered more tasks that can be accomplished with cursors, and also went into how text files can be created and managed with Python scripting.  This is really useful, because much of the data we deal with in GIS comes to us in the form of text files.

We often must break down shapefiles and feature classes to their most basic level, which is that of the vertices.  This can be done with search cursors and nested for loops.  The initial for loop will return the features or rows, the second will return the parts of features, and the last will produces the individual points or vertices that make up a line or polygon.

In the exercises, we practiced converting text files into shapefiles.  In the assignment, we took a shapefile and populated a new text file with its vertices.  Because the feature class consists of several river features, each of which is comprised of a number of points or vertices, this script requires a nested loop structure, in which rows with OID, X and Y coordinates, and feature name are retrieved  by the search cursor, and inside that loop, the points are retrieved for each array of each feature.  The  X and Y coordinates for each vertex are then written to a text file, along with the OID and Name of each feature for the rivers shapefile.

Below, left, is a portion of the resultant text file.  Each line includes the feature ID number (OID), vertex ID number within that array, X coordinate, Y coordinate, and feature name.

Besides writing to the text file, I also had my script print out the same data for each vertex, after each iteration of the second (inside) for loop.  This way, the progress of the script can be followed as it writes data into the text file.  A part of the Results from the Interactive Window is shown on the right of the text file.

Fig.1.  Text File with vertex data for river.shp

Fig. 2.  Results  (partial) showing script progress

This is the pseudo code for the script.  These are the logical steps, more-or-less in English, for the Python script.

Pseudocode.

import arcpy site package

import env module

env  =  set workspace to …\Data

Set to override old output

fc = river.shp

Search Cursor = retrieve the fields OID, SHAPE@ (all geometry), and NAME for fc shapefile

f = open/create the new text file

Loop 1: for each row in cursor

      v = vertex = 0       Create an ID number for each vertex,   start over at zero for each OID.

      Loop 2:  For each row, get arrays         (for part in row .getPart )  

          v +1        add 1 to the vertex number

         write each point to textfile as a separate line:

         f.write () =  strings:  row OID  +  vertex number + X coord  + Y coord  +  row NAME + line break

f.close textfile

GIS Programming - Module 7 - Exploring and Manipulating Spatial Data

This assignment covered two chapters (6 and 7) in the textbook, and involved more work, over two weeks.  It is so far the most challenging of the modules of this course but potentially the most useful as well.  We worked with lists, dictionaries, and search, input and update cursors.

A dictionary is a list of pairs of objects.  The first of each pair is called the key, and it is associated with a value.  In our example, the keys were county seat cities, and the values were their population values.  A key can be associated with more than one value, but not vice-versa.  A dictionary can be used to look up a value for an object, such as the population for a city.

A cursor is a data access object.  A search cursor can be used to retrieve certain fields from a dataset, with the help of SQL statements.  Our script involved using a search cursor to retrieve the county seat cities in New Mexico, with names and populations, add them to a dictionary, then print out the dictionary.  A delete cursor finds and deletes fields, and/or data based on an SQL statement, from a feature class.  An update cursor can add or change data.

In our assignment, we used a search cursor to find the names those cities in New Mexico that are county seats, along with their populations.  Then, using a for loop, we updated a new, empty dictionary, adding pairs of key and value with the name and population of each city.

The following figures show the results from my script.  Because this script is fairly long and involved, it details the steps of the process, and prints statements after each process is accomplished, and finally prints the resultant dictionary of county seats and their populations.

(Because these results are long, they're shown in three figures, although in actuality the whole thing was produced in a single run of the script, in the Interactive Window.)

Figure 1. Part 1: Results for Creation of Geodatabase, List,
Initiation of Copy to Geodatabase

Figure 2. Part 2: Results of Copy  (continued)

Figure 3. Part 3:  Results for Populate Dictionary
with Search Cursor for loop.  Print Dictionary.

Process Summary Details.

Question: Which step did you have the most difficulty with? 
Describe:
1) the problem you were having, and 
2) the solution or correct steps to fix it.

I had some trouble with Step 7, populating the dictionary with the search cursor.  My main trouble was understanding the logic of the Search Cursor, and exactly what it’s supposed to do.  Once I read the text and Help topics and some of the Discussion post answers, I figured out the purpose and then was able to put together the logic in the form of pseudo code.  I then understood the idea of populating the dictionary with the rows, and began to figure out what exactly the row variable represents.    
My first version looked like this:
cursor = arcpy.da.SearchCursor(fc,["NAME", "FEATURE", "POP_2000"],'"FEATURE" = \'County Seat\'')
for row in cursor:
    print row[0], row[1], row[2]
This was before I had gotten to the step where I created the dictionary, and the result was simply a list of the names, feature type and population for those county seat cities.  
Then I created the empty dictionary, but had some trouble figuring out where to put the update statement (from the assignment instructions) in relation to the search cursor statement.  For several tries, the script ran okay, but I didn’t get the dictionary to populate.  It remained empty. Unfortunately, I didn’t save this incorrect version of the script and can’t remember exactly what it looked like. I initially put the dictionary update statement outside the for loop, and put print statements for the rows inside.  This ran okay, but because the rows weren’t being iterated through the update function, the dictionary wasn’t being updated at all.  
After looking at the ArcGIS Help page, and the text, and the discussion posts, it finally began to dawn on me just what the search cursor is supposed to do, and realized that the update statement needed to be on the inside of the for  loop of the cursor statement.  At the same time, I also started paying attention to the fact that the [0] and [1] and [2] etc. represent index positions, reviewed that concept a little,  and saw that I needed to have the correct ones for each key and value that I wanted to add to the dictionary from the search cursor rows.  Another thing I realized, is that row is actually a variable.    After studying the syntax some more, I took an experimental stab at putting the update statement inside the for loop, and luckily it worked.  I found that the trickiest part was configuring the part inside the parentheses properly.  It was necessary to combine the syntax of the row variable with that of the dictionary.
I have to emphasize here that part of my understanding came AFTER this, when I saw that it worked.  
At this point, I still had the print county_seats  statement inside the for loop as well, and this resulted in the dictionary being printed repeatedly, and longer, every time a new key and value pair was added.  I’ve had enough practice at this point, however, that I knew right away that the print statement needed to be un-indented, so I did that.  
My final script that works properly is shown here:

for row in cursor:
    county_seats.update({row[0]: row[2]})  
print county_seats

row 0, the key, is first index position: NAME.  row 2, the value, is third index position: POP_2000.  The dictionary is inside the curly brackets, with the key and value separated by a colon ( : ).
It would be very useful to save scripts that run okay, but don’t produce the correct results.  I plan to do this in the future, and include commentary about what the script fails to do, as well as what it does do.  


  I initially had a section at the beginning of the script that looked for the existence of the geodatabase, and deleted it if it was there.  I did this because I thought that was what caused the error statement that the .gdb already exists.  Then I realized it’s okay if the geodatabase already exists, and that it will be overwritten, but only if ArcMap is closed, and there’s no lock on the geodatabase.  When I figured this out, I got rid of the Exist – Delete part of the script.