Lab 1 - Calculating Metrics for Spatial Data Quality

This first lab was about calculating metrics for spatial data. One of the first things to get a handle on was the difference between accuracy and precision. Accuracy is basically how close the measured or observed value is to the known or “true” value, where precision is how close or clustered the measured or observed values are to each other. The first part of the lab introduced us to using ArcMap to determine precision and accuracy. Using a set of 50 different points representing a single location mapped 50 times using a GPS unit, we performed tasks to determine the accuracy and precision of the GPS unit. First, we determined the average location by using the Statistics tool to obtain an average X and Y location, and created a feature class, using Editor to place this point on our map. The shapefiles needed to be projected into a coordinate system that uses meters instead of decimal degrees. Next, we needed to find a distance from the average location that corresponded to a percentage of the observations, specifically 50%, 68%, and 95%. To do this, I first performed a spatial join on the waypoints and average location shapefiles, which created a distance field. After sorting the distance field in ascending order, I set my buffers at the distance corresponding to the 25^th value (50%), the 34^th value (68%), and I used the value halfway between the 47^th and 48^th value (95%). This created buffers for the 50^th, 68^th, and 95^th percentile, which shows the precision of the results. The map displaying this is shown below.

Using the 68^th percentile as an accepted value for precision, I compared the device’s accuracy and precision. In the horizontal direction, the distance between the average location and the reference point was determined (using the Measure tool) to be 3.8 meters. The precision value was 4.4 meters. As the distance between the average and the known is less than the precision value, I determined that it is reasonably accurate in the horizontal direction. In the vertical direction, the distance between the average and the reference value was approximately 6 meters, where the precision value is 5.7 meters. This tells me that the device is not as accurate in the vertical direction.

In Part B, we worked with a larger dataset and determined the Root Mean Square Error (RMSE) and created a cumulative distribution function (CDF) of the error. After opening the Excel file and copying the benchmark X and Y values over, I had columns for various points (X and Y values) and the benchmark X and Y values, from which I created new columns and calculated the X and Y errors (difference between the point values and benchmark values), the XY error squared, and the error XY values. From the error XY values, I calculated the mean, median, RMSE, 68^th, 90^th, and 95^th percentiles. I also determined the minimum and maximum.

From the error XY values and the cumulative percent (there are 200 values, so each is 0.5 percent of the total), I created a scatterplot showing the Error_XY versus the cumulative percent, which helps to demonstrate that several of the metrics can be determined from the scatterplot chart.

The most difficult part of the assignment was determining the difference between accuracy and precision. The definition is straightforward enough, but I was having issues mentally picturing the difference when performing the calculations. It made sense to me to say that if the distance between the average location and the reference point was less than the 68% precision value, then the data is accurate, but it’s still a little unclear. Hopefully I can work a little more with the statistical aspect of GIS in the future and gain a better understanding of it.

Module 11 - Sharing Tools

This was our last assignment of the course, and in this one we work with Python tools again. Last week we learned to create tools, and this week we learned to share them. This is somewhat of a continuation from last week. In this lab we modified script so that the parameters would use file paths set in the tool itself instead of from the set file paths in the script. We did this by using the sys.argv[] code. The positions numbers for this start at "1", as opposed to the arcpy.GetParameter() code, which starts at "0".

We learned to edit the item description of the script tool in ArcCatalog. We learned to embed the script into the tool by right-clicking the tool in ArcCatalog and selected "Import Script." One of the more interesting features we can do now that we've imported the script is we can now set a password. The password still allows the tool to run, but the password needs to be input to view or edit the script. The tool in the assignment creates a number of random points within a feature, and then creates a buffer around those points. Below are the results of both my tool dialog and the results in ArcMap.

You can see the feature in purple with the random points (in this case, 50 of them) in black, with the light blue buffers around them.

I really enjoyed this course, and I learned a lot from it. I would have liked to have done a little more of script writing, but this course is designed to use Python with ArcGIS, and I think it teaches that very well. I learned how to use Python to perform geoprocessing tasks, and about toolboxes and tools within Python, working with rasters and different geometries, etc. I think this course is a great start as one moves forward to work more with GIS and Python programming.

Participation Assignment #2

The article that I found is an interesting application to GIS. This describes how, in 1996, TransCanada used GIS to select the best route for an expansion pipeline. I found it interesting because this seems to be one of the first times this type of analysis took place using GIS. It's a route selection analysis very similar to what we did early in the semester in the other course with the corridor analysis. They set the project goals, which were threefold. First, minimize risk to public safety. Second, minimize environmental impact. Third, minimize construction and operation costs.

This particular project used information from already existing digital and map sources and interpreted imagery. They used GIS techniques to include “exclusions, constraints, and opportunities.” When minimizing environmental impacts, they weighted environmentally sensitive areas with heavier constraints. Collocations with already existing rights-of-way were weighted with lower constraints, so the GIS could label them as opportunities. They also weighted river and road crossings heavier, as that would increase construction costs. They also created buffers around urban areas to minimize the impact on the population.

Using these criteria, TransCanada developed several route alternatives. The GIS was able to provide the lengths of the routes and detailed crossings information. They also analyzed the impacts of the different criteria on the various alternative routes. Before selecting a preferred route, more data were gathered based on aerial imagery and ground-truthing. This method became a go to method for other entities for route selection, and it is very similar to scenarios we have seen in the other course.

GIS-based process helps TransCanada select best route for expansion line

Module 10 - Creating Custom Tools

In this assignment, we learned to create custom tools and toolboxes. First, we were create a script tool from a provided script, which was to run in ArcMap. I opened the script and modified the script name, last modified date, and modified by comments, and saved it under my name. In ArcMap (ArcToolbox) I created a new toolbox, and within that toolbox, a new script. Creating a new toolbox was rather straightforward, was just right-click -> New Toolbox; very similar to making a new folder. Next we were to set parameters. Creating the parameters is pretty simple tool, as in Properties there is a tool dialog box that makes it easy to not only create parameters, but to modify their settings (data type, multivalue property, etc.). Below is a screenshot of my script tool window.

Now that the parameters are set in the script tool, they needed to be set in the script. I changed the original four variables in the script, replacing them with the arcpy.GetParameter() command. I had a little bit of an issue here, as I was using the arcpy.GetParameterAsText() command initially. This didn’t work, as this only reads parameters as text, which caused errors. As the outputFolder can’t be concatenated with strings, I used the str() function on the outputFolder variable, and ran the tool. Next we were to replace the print statements with arcpy.AddMessage() statements. Initially, I did this and no matter what I did, I could not get the code to work when running it. The problem turned out not to be in the code. What was going on is that I was trying to run the script in PythonWin as a stand-alone script. The arcpy.AddMessage() statements are, I think, designed to be run as a script tool. When I opened ArcMap and ran the tool again, it worked fine. Below is a screenshot of my results window.

To share the toolbox, tool, and script, I selected both the toolbox and the script and sent them to a compressed (zipped) folder.

I feel I learned a lot in this lab, and I think I’ll get a lot of use out of using toolboxes and tools in the future. I think one of the most valuable components of this lab is the idea of sharing the tools and toolboxes. Once I figured out the couple of issues I had and why what I was doing was causing problems, the assignment seemed pretty straightforward. Creating toolboxes is rather user-friendly and should be picked up quickly by most people familiar with Windows. The same can be said for setting the parameters in the dialog box, although it’s important to get the settings such as “Multivalue” correct, or the tool won’t run properly; it is easy to go back and check those, though. Once I understood what the arcpy.GetParameter() and arcpy.AddMessage() functions are doing, I realize they are quite useful and I’ll try to find ways to use them in the future. I also like the idea of sharing toolboxes and scripts, so others can use or modify them.

Module 9 - Working with Rasters

This assignment taught us to work with rasters using Python. We worked with 2 rasters, elevation and landcover, to create a Python script is a single raster image with the following criteria: 1) Forest landcover classification; 2) slope between 5 and 20 degrees; and 3) aspect between 150 and 270 degrees.

It is becoming almost habit to start off by importing the necessary modules needed for the script and to set the overwrite output to true. For this assignment, I attempted to create the output in a .tif format, but I got several errors, so I went with the optional part and created a file geodatabase. For the rest of the script, the module exercise was extremely helpful and relevant to the assignment. I created a conditional statement so the script would only run if the Spatial Analyst extension was available. Further down toward the end of the script, there is an else statement printing a message if it is not available. Within the if statement is where the bulk of the script is located, as I only want it to run if the spatial analyst module is available. I used the remap value and reclassify tools to reclassify the landcover raster so that only the forested areas are given a value of "1". I assigned the elevation raster to a variable and used the slope and aspect functions to create temporary rasters of those types. I was impressed with how user friendly the map algebra commands are when used in Python, so performing the slope and aspect calculations weren't an issue. I created temporary rasters based on the 3 criteria listed above and using "less than" and greater than" statements in the script. I combined the 5 rasters (four slope and aspect rasters and the reclassified landcover raster), and saved it to make it a permanent raster. Below is a screenshot of the raster, with the "0" (reddish-brown in color) signifying areas of the landcover raster not classified as forest and areas labeled as "1" (the green areas) showing areas that are forest land cover.

This assignment was far less stressful for me than last week's assignment. I feel I have a greater understanding of this concept, and the exercise was extremely helpful, so if I had a syntax error I could look at the exercise (and the textbook) to determine if I was missing a parentheses or whatever the error was. The only issue I had dealt with the file geodatabase I created. Even setting the overwrite output to true does not allow you to overwrite the geodatabase, which I suppose makes sense. The geodatabase has files within it, and you usually can't overwrite a non-empty folder. I was unable to find any workaround, so after running the program once, you still have to go into the Results folder and delete the geodatabase. I didn't do this as it just occurred to me while writing this, but I wonder if it's possible to have the code delete the files within the geodatabase at the end of the script, and then have Python overwrite the empty geodatabase. I liked working with rasters in this assignment and I was impressed that Python almost seems to be made to work with some of the map algebra operators. I like that a relatively short and simple code can do this much geoprocessing. On another note, the ELearning site appears to be down, so I chose to email the instructor the relevant deliverables.

Module 8 - Working with Geometries

This week’s lab dealt with working with geometries in Python, and it was a challenging one for me, partially because working with arrays and nested loops are not my strong suit and partly due to a pretty intense work schedule this week. In this lab, we were working with the “rivers.shp” file, using loops to iterate through the cursor, and we wrote the data to a .txt file.

After importing the modules and setting the workspace, I created a search cursor for the shapefile I was working with. This called the OID field, the SHAPE geometry object, and the NAME field. I created a .txt file using the write mode, following examples in the text and in our exercise. To write the data to a text file, I needed to create a couple of “for” loops. The first iterates through each row in the shapefile’s attribute table. I created a variable called vertexid within the loop set to 0, so that it resets for each new feature. Within that loop, I created a second loop using the .getPart() method, which iterates through each point in the array of the row. The .getPart() method allows me to access the points in each array. If I did not use this method, another “for” loop would be necessary. I added “1” to the vertexid variable in this loop to keep track of the number of points in each feature.

Here’s where it got tricky for me. I needed to use the write() method to write to the .txt file. I wanted to add the Feature/row ID, the vertex ID, the X and Y coordinates, and the name of the river feature, using “\n” to ensure each vertex is on its own line. I ran into a wall here with syntax errors and exceptions at first. I was pretty sure the loops were correct (at least mostly) because I immediately was able to print correct values for the X and Y coordinates. I ran into a syntax error which gave no input into where or what the error was, so I started deleting code (remembering later I could have commented it out) and attempted to run the program. The issue turned out to be a print statement I had created underneath the write() statement, which wasn’t necessary, so I just left it out. The other error I had was not using arrays properly. I had the output in the .txt file correct, but the name and feature were printing “NAME” and “FEATURE”. Once I changed this to row[2] in the write() statement, the program wrote to the .txt file correctly. The output is shown below.  The first column shows the Feature/row ID, the second column shows the vertex ID for that feature (notice how it resets for each feature), the 3^rd and 4^th columns are the X and Y coordinates of the vertex, and the 5^th and 6^th columns are the feature name and type.

This was a really good lab, as I learned about working with geometries and got more experience working with arrays. Both the instructor and the other students on the discussion boards were a big help with steering me in the right direction with my code. I wish I had had more time to really work with it, as opposed to an hour or two at a time, but even so, I learned a lot and I think it’s less confusing and less intimidating to work with than with FORTRAN or even IDL code.

Lab 8 - Damage Assessment

In this lab, we carried out a building damage assessment based on aerial imagery taken before and after Superstorm Sandy impacted the northeastern United States. Although not especially powerful as far as some hurricanes can be, the extreme size of the storm, the fact that it collided with a powerful weather system moving into the northeast, and its impact on several of the U.S.’s largest cities (including a historic storm surge), combined to make it the second costliest cyclone to hit the United States since 1900.

First I created a base map using the World Countries, US boundary, and FEMA States shapefiles. The FEMA shapefile was created from a selection of impacted states from the US boundary shapefile. After importing the Sandy storm track table, I created a polyline showing the storm’s path from the original hurricane track points. I learned to create marker symbols, and symbolized the hurricane with a hurricane symbol. I used a color scheme from green to red, with green showing where Sandy was less powerful and red at the storm’s most powerful, in this case a Category 2 hurricane. I labeled the hurricane track points and added graticules and the essential map features. Below is my final map.

Next I needed to prepare the data for the damage assessment. This was achieved by creating a new file geodatabase and new raster mosaics, one for the pre-storm imagery and one for the post-storm imagery. I briefly learned to compare the before and after imagery by learning the Flicker and Swipe tools in the Effects toolbar. One of the most important aspects of the assignment was creating attribute domains, which are used to constrain values allowed in an attribute table or feature class. Shown below is a screenshot of domains I created for the damage assessment. Damage ranges from 0 (no damage) to 4 (total destruction).

I added the county parcels and the study area layers to our map. This part was quite time consuming, as I had to digitize each structure from the pre-storm image within the study area, and then compare the pre-storm and post-storm imagery to determine the level of damage the structure received. This was often difficult to determine, as I was using aerial imagery and it was difficult to see much of the damage unless it was very severe. It was also often difficult to determine the building type (residential, commercial, or industrial).

Next I wanted to examine the relationship between the locations of the damaged structures and their proximity to the coastline. I digitized the coastline near the study area and used Select by Location to determine the number of damaged structures within 100 m, 100-200 m, and 200-300 m from the coastline. As one would imagine, structures closer to the coastline were more likely to be destroyed or suffer catastrophic damage, and those further from the coastline were more likely to suffer little to no damage. Structures blocked from the winds and/or water had a better chance of surviving intact as well. Below is a table of the structural damage.

                                                 Count of structures within distance categories from the coastline

Structural Damage Category
	0-100 m	100-200 m	200-300 m
No Damage	0	0	7
Affected	0	6	7
Minor Damage	0	26	22
Major Damage	2	6	2
Destroyed	10	10	4
Total	12	48	42

I used the Attribute Transfer tool to copy the attributes from the county parcels layer to the structure damage layer. To populate the new fields, I manually matched the parcel to the damaged structure point. Once I had done this for all the damage points, I exported the data to an Excel spreadsheet and labeled the points. The purpose of this was to help determine who owns the damaged structures for emergency management and insurance purposes.

This lab seems to cover most of what damage assessment is used for. Although manually digitizing is rather tedious, it helps determine the extent of the damage. In a real situation however, I would be more comfortable creating the initial assessment using aerial imagery, but actually going through the area to visually confirm the assessment values would be important as well.