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.