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.  

GIS Programming - Module 6 - Geoprocessing with Python

Python Results for All Three Tasks

 We used ArcGIS tools again this week to construct a script, but instead of using the ModelBuilder to put it together, we went directly to PythonWin and entered lines of code with ArcGIS tool functions.

There are very useful Help pages for each tool, that can be found via the ArcMap toolbox, or the Python window in ArcMap.  Correct syntax and parameters, as well as examples, are shown for every tool.








The stand-alone Python script of this week's assignment had three tasks to accomplish on a point shapefile representing locations of hospitals:

1) Add XY coordinates;
2) Create 1000 meter buffers around each hospital point;
3) Dissolve the multiple buffers into a single large buffer.

The figure above shows the results generated in the Interactive Window of PythonWin.  To obtain the results like this, three additional lines of code must be added to the script (one after each function or tool code), using the GetMessages function. This is a good way to keep track of what you've done, and document that the various parts of the script have run successfully.  The parameters that were used are also shown.

Results of Task 1

This figure shows the result of the first task in ArcMap, to add XY coordinates to each record in the hospitals.shp file. A few of the new fields are seen here in the attribute table. Each row is a record representing a hospital; other previously-existing data fields (not shown) are to the left in this table.

Results of Task 2

This shows the resulting output from the Buffer task (2) in the script.  A buffer zone with radius of 1000 meters was added around each small point representing a hospital.       The new polygon shapefile that was generated is called hosp_buffer.shp.  This shapefile has multiple features, one for each hospital point.

Results of Task 3

This figure shows the results in the ArcMap display of the last task executed by the script, to dissolve the multiple buffers into a single feature.  This new output polygon shapefile is called hosp_dis_buff.shp.

This feature represents all areas that lie within 1000 meters of a hospital.









As can be seen from its attribute table below, this shapefile has only a single feature, comprising all 1000 meter buffers for the hospitals.

Single Feature of Dissolve

Below is the detailed Process Summary for writing this script.

1.      I opened a new, blank Python script in PythonWin and saved it, then added the commentary at the top regarding name, date, etc.

2.      Before anything else can work, you must enter this line:  import arcpy, to have access to the functions of ArcGIS site package.  Also enter this line of code: from arcpy import env.  This imports the environment module, so that the environment (including file path) can be set.

3.      Two other lines to enter at the beginning of any script will 1) set the environment (specify the file path) and 2) make it possible for Python to overwrite previously generated results.  These commands in my script were:  

env.workspace = "S:/GISProgramming/Module6/Data"  

env.overwriteOutput = True

4.      The three tasks that this script had to complete were to 1) Add XY coordinates to the hospitals.shp layer 2) Make a 1000 meter buffer around all the hospital point features and 3) Dissolve the buffers into a single feature. 

5.      After the execution of each tool, the script will print the messages from that tool.

6.      To find out how to make the script for each tool, I opened each one in ArcMap, went to Tool Help, and looked at the script syntax and code samples.  Each parameter is also explained there.

7.      There are two parts to the code samples: one is the most basic form, and the other is an entire stand-alone script in PythonWin.  Both are useful.  The most basic form lists the parameters in order, by name.  The more involved script shows how the lines of code should be written, with respect to file paths and layer names.  This part also gives a few ideas of how to improve the script, such as by making extra copies of layers that will be changed. 

8.      Before executing the AddXY function, I did a Copy function on the hospitals.shp layer, to preserve an original version of that layer in the Data folder. This layer is called original_hospitals.shp.      The hospital.shp file itself was then altered by the AddXY tool/function. 

9.      Copy function:  arcpy.Copy_management(in_data, orig_data)

10.   But first (before steps 8-9), I assigned variables to the hospitals.shp (in_data) and the original_hospitals.shp (orig_data).  This was necessarily done before the Copy function was executed. 

11.   I then ran the AddXY function on in_data (the variable assigned to hospitals.shp) which added two new columns to the attribute table, for X and Y coordinates.

12.   The code is:    arcpy.AddXY_management(in_data)

13.   After this ran, I added the line of code:

print (arcpy.GetMessages() + "\n")

14.   This causes the results messages to be printed out for the tool that just ran and 

+ "\n" after the main part of the code causes a blank line to be returned after the messages, to enhance readability. 

15.   I added this line of code: print +"\n" in other parts of the script: this does the same thing. That serves to separate the parts and make the results easier to read. 

16.   After a couple of blank lines were returned, I added the code: print “Task 1 complete.”

17.   Before the next task results, I had the code return two more blank lines (as in step 13).

18.   Task number 2 was to create a 1000 meter buffer around each hospitals.shp point feature.

19.   First, I assigned two more variables:  buffer is for the new output shapefile that will be created in the Results folder.  Its line of code is:

buffer = "S:/GISProgramming/Module6/Results/hosp_buffer"

20.   hosp_buffer is the output file. 

21.   The path could also have read:

buffer = "../Results/hosp_buffer" instead of the entire path.  This is the syntax for a different folder at the same level (ex. Results) within the same higher-up folder as the folder specified in the environments (Data).  In other words, to go from Module6/Data to Module6/Results.

22.   The distance for the buffer is assigned the other variable: distance = "1000 meters"

23.   The next line is the Buffer function code: 

24.   arcpy.Buffer_analysis(in_data, buffer, distance)

25.   Again, the messages for the tool execution are printed using the GetMessages function, blank lines are added, and the line “Task 2 completed.” is returned.

26.   The last task is to dissolve the multiple hospital buffers into a single feature. 

27.   The code is:

28.   arcpy.Dissolve_management(buffer, dissolve)

29.   The default parameter for multi-part is taken, so no further optional parameters are listed.  This causes all the multiple hospital buffer features to be dissolved into one single feature for the hosp_dis_buffer.shp feature. 

30.   Again, the messages are printed for this tool’s execution, using the GetMessages function. 

31.   print “Task 3 complete.” and 2 blank lines are returned.

32.   To signify that the entire script has run, we have one last line:

33.   print “All processes complete.”

34.   A couple more print “\n” lines are added at the bottom, to separate multiple runs in the Interactive Window.  

GIS Programming - Module 5 - Geoprocessing in ArcGIS

This week's assignment really starts to bring things together, between ArcMap and Python.  We used ModelBuilder from within ArcMap to clip a soil type shapefile to the extent of a basin, then selected the soil features with the attribute "Not prime farmland" then erased those selected features from the clipped soils shapefile to create a new shapefile, shown at left.  This shapefile represents the tracts of soils within the basin that are considered some type of suitable farmland (in other words, everything except "Not prime farmland" tracts).

The Select tool, when operated from within ModelBuilder, has a problem in that it always changed the SQL term 'Not prime farmland' to 'Not Prime Farmland'.  This caused none of those features to be selected, because the lower-case letters had been changed to upper-case, and so were not recognized.  The Erase too
l then had nothing to erase, because the Select tool had found nothing to select.  Strangely enough, when the single quotes around Not prime farmland were changed to double quotes, and the model run (with errors, because it didn't like the double quotes), and then the whole SQL statement was entered again, the lower case letters were preserved in the model, and it could run correctly.  Some other people also had this problem and brought it up on this week's discussion board, and found that removing the quotes completely from Not prime farmland, and then putting them back, also fixed the problem.

The model was then converted to a Python script, and a few adjustments were made, such as adding the full file path to the input shapefiles (soils and basin) so that the script could find them, and also adding the line of script arcpy.env.overwriteOutput = True.  This makes it possible to run the script multiple times, although the outputs already exist in ArcCatelog from previous runs.  The Python script can then be converted into a Script Tool within ArcMap, where it can be run like any other tool.

Process Summary Notes:

Part 1, Step 2: Explain how you set up your SoilErase model.
1. First, open a new ArcMap .mxd, and add the shapefiles soils and basin.
2. To create a model, open ArcCatalog, right-click the folder in which you are working (ex. Module5), select New, and then Toolbox.
3. Rename the toolbox appropriately, then right-click it, select New, and then Model.
4. This will cause a ModelBuilder window to open.  In the top menu, select Model¸ then Save As, then name the model.
5. For this model, the task is to clip the Soils shapefile to the extent of the Basin shapefile, then remove the features within Soils that have the attribute “Not prime farmland.”
6. From the ToC on the .mxd, click and drag both Soils and Basin shapefiles into the ModelBuilder window.  These will be the inputs, and they will show up as blue oval boxes, with the names of the shapefiles inside.
7. The Clip and Select tools are dragged from the ArcToolbox  > Analysis Tools > Extract toolset, into the ModelBuilder space.  They will show up as yellow, rectangular boxes.
8. The Connect tool is used to draw arrows from the input shapefiles, to the tools.  Although Select or Clip can done in either order, I did Clip first.
9. As I drew the arrow from Soils to Clip tool, a box appeared, in which either input feature or clip feature could be indicated.  Because the Soils is being clipped, it is the input feature.
10. Then, another arrow is drawn from Basin to Clip, and in this case, the clip feature option is selected.
11. After this is done, a new box appears for the output, supplied with a default name.  This box is a green oval, indicating that it is an output.
12. The Select tool box is set to the right of the new output from the clip tool, and an arrow is drawn from the clip output, to the select tool box.  Parameters must be set; this can be done by double-clicking the Select box.  In the SQL dialog, enter, “FARMLNDCL” = ‘Not prime farmland”
13. A new output box (green oval) appears, which contains the selected features from the last output.
14. The Erase tool is dragged into the ModelBuilder, and an arrow is drawn from the last (Select) output box, to the Erase tool box.  The input feature is the clipped soil shapefile, and the Erase feature should be the Selection of the last output (those tracts that are Not prime farmland.)
15. The new output feature, and the final one, is the soil tracts that fall within the extent of the basin, with those of attribute “Not prime farmland” removed.
16. After this is done, right click the final output, after the Erase tool, and select “Add to Display.”
17. On each of the 2 inputs (blue ovals) and 3 outputs (green ovals), right-click, and check “Model Parameter.”  A “P” will appear next to each of those boxes. (This is not done to the tools.)  This will make it possible to change the input features in the model, as well as the file paths, and make the model useful for general purposes.
18. From the top menu: Model > Save.

You are also strongly encouraged to include any other notes you have on this week’s lecture and lab or your process.
(These extra notes are optional, and are not required for any part of your grade. However, they can greatly improve your learning experience by reinforcing what you’ve learned, and giving you a reference to look back on in the future. Plus, they can make your blog posts more interesting!)
Other notes:


1. After model is exported to Python script form, you must
2. add the “import arcpy” line of code at the top, and then

3. arcpy.env.overwriteOutput=True

4. you must also make sure each shapefile’s entire file path is specified for the variables.

This problem exists for the Select tool, when it is used from within the ModelBuilder: the attribute in the soils shapefile is "Not prime farmlands," with the p and f in lower case.  However, when you enter that into the SQL statement, inside ModelBuilder, either by typing it in, or selecting it from the "Get Unique Values," it always reverts to "Not Prime Farmland," with capital P and F, even if you change it back to lower case and save it.  This is not the correct name for the attribute, so nothing is selected.    I've included a screen shot.
  I tried running this process manually, using the Clip, Select and Erase tools separately, directly from the toolbox.  In this case, the problem DID NOT happen.  The p and f in Not prime farmland did NOT revert to caps, the selection was created, and it was properly erased from the soils shapefile.

took out the quotes, put em back in again, erased the whole statement and added it in again, and now the select works and the Not prime farmlands gets properly erased!

GIS Programming - Module 4 - Debugging and Error Handling

This week in GIS Programming, we practiced locating, fixing, and otherwise dealing with various types of errors in Python scripts. These skills help make that formidable Python much less fearsome!   

There are generally three types of errors that you might make when writing script:

1) syntax errors include spelling and punctuation mistakes.  Making typos, leaving out a colon (:), or indenting incorrectly are very common.
There are also:
2) exceptions, for example, in which modules and functions (and other parts of Python) don't have what they need to run properly.  There are many, many types of exceptions, and luckily they are identified in a very general sense in the error statements in the results, and the line numbers in which they occur can be found as well.

3) The other type of error can be very insidious: the logic error.  In this one, your code might run fine, but the answer you get is wrong, because you've told Python to do something different than what you intended. Python will not know the difference; it's only following directions, and can't read our minds!  This is an example of GIGO: garbage in: garbage out.

For our assignment, we looked at three scripts with syntax and exception errors.  In Scripts 1 and 2, we repaired the errors, but in Script 3, we added a "try-except" loop, which gives Python an alternative path, in order to step around the error.  We didn't actually fix the error, and the results tell us that it's there.

You don't necessarily have to repair an error right away. It is possible to use a "try-except" statement around a block with an error, which isolates and bypasses the error and goes on through the rest of the script.  It "tries" the problematic line, then has an "except" route, which offers an alternative for Python, such as printing out a statement identifying the nature of the error.  This way Python doesn't have to get hung up in the middle of the script.

I found that there are three excellent ways to find errors. The first two are rather rudimentary.  1) You can "comment out" most of the lines in the script.  This means you type # in front of every line except the first one, then run it.  If it runs okay, take the comment # off the next line or loop (for, while, etc) and run it again.  You keep working your way downward through the script, de-commenting the lines one by one, until it runs properly.   This way, you can eliminate the lines that are okay, and find and fix the ones that aren't.
2) Another way to check and see if a block of script is working right, if there are no results following it, is to insert a "print" statement in a new line after it.  I used "print okay."  This way, if "okay" appeared in the interactive window, I knew that everything above that print statement was okay.

3)PythonWin also has debugging tools that work in a similar way, but without having to comment out or add print statements, which can be cumbersome if you have a long script.  You can turn on the debugger, then run the script in a "step-through" fashion one line or block at a time.  At each step, you'll see that it was okay, or that it wasn't.   This way, you can find the errors or confirm that a line is okay.  You can also insert break points, which isolate the parts you want to run.

Below are screenshots of my results from the three assigned scripts. They are shown in the PythonWin Interactive Window, after I had repaired the errors or added a try-except statement: 

Script 1.  These are the results, after the script was corrected.

Script 1 Results

Script 2: This script had 8 errors.  These are the results, after they were identified and fixed.

Script 2 Results

Script 3:
In this script, the error was not corrected.  Rather, a try-except statement was inserted after the error was identified, and the results printed out the type of error occurring in Part A of the script, then went on to Part B and printed the results.  Part B did not contain any errors, and was not directly affected by the error in Part A.  However, until the error in Part A was dealt with, or "caught," the script would not continue to execute, and would not have been able to proceed to Part B and produce results.  

Script 3 Results

GIS Programming - Module 3 - Python Fundamentals Part 2

Results from the Interactive Window for the Module 3 Scripts

In the second week of learning the fundamentals of Python, we examined a script with an if conditional statement, and identified and corrected a couple of errors that were preventing the script from running.


We also wrote two scripts of our own.  The first script created a list of 20 randomly-generated numbers,with values of 0 to 10, then printed the whole list.








In the second script, we specified one "unlucky" number, then had the script remove all instances of that number, if it was there, from the list we'd made in the previous step.  We also had the script print messages that stated if the number was in the list, or not, and if so, how many times.  At the end, the new list without the unlucky number was printed.

Here are some notes from my Process Summary for this week.  This describes generally how the unlucky number can be removed from the script in the last block.

Removing numbers from the list. 

1.      For this part of the code, you already have the list with the random numbers appended into it.

2.      A variable has already been assigned to the unlucky number.

3.      All instances of the unlucky number are removed from the list with a while loop. 

4.      As long as the list still contains instances of the unlucky number, the while loop will continue to repeat and remove the numbers; this is the condition. 

5.      The indented statement within the while loop that actually removes the number is the remove method: the object is the list, and the argument is the variable assigned to the unlucky number. 

6.      In my script, the method is: SeaseList.remove (unlucky).

This script has a system argument: the user must input the unlucky number variable in the Run Script box.