Introduction to Geographic Information Systems

ERE 450/550

Lab Exercise 3

Lab due for ALL sections on Friday October 20 at 5pm in Bray 410

Purpose

Elizabeth Nifkin has requested a memo, map, and table communicating the results of the soil erosion question: Where are areas of potential soil erosion problems for the city of Syracuse? For course purposes you will also do a spatial model (flow chart) and answer sheet for this exercise.

You will use the Revised Universal Soil Loss Equation (RUSLE) to compute the annual soil loss: A(ton/acre/year) = RKLSCP This analysis will be done using the NAD27 datum and the UTM Zone 18N projected coordinate system. The variables of the model are as follows: the R factor from rain energy, K factor from soils, LS factor from topography, C factor from land cover, P factor from erosion control or management practices (we’ll assume that there are none), and finally our geographic unit of analysis, the City of Syracuse. The information for doing this lab was obtained from New York State DEC’s stormwater management resources, the work of Helena Mitasova et. al. (1997, 1998), and the USDA’s USLE handbook (Wischmeier and Smith, 1978).

You are doing this work for Lizzy, but while she is GIS savvy, she will pass the memo on to those who don’t use GIS, so be careful of technical language.

You thought you had done enough data prep in Lab 2, but we’re still a ways off from calculating areas of erosion. First, we need to think about how this will work. First, which model to use? Vector or Raster? Because we are going to multiply layers together, the easiest model is the Raster model. Now let’s think about the variables of the RUSLE model. You need an R-value or an erosion factor from rain. For this one we don’t need a layer, because the R-value for Onondaga county is 80. The K-factor is the erodibility of the soil. You will use a table from the STATSGO data set, do some calculations, and then join that table to the clipped soils polygons. From there you will create a raster layer of K-values. Using the DEM you will create slope, another related layer (flow accumulation), and an equation to produce the LS-factor. The C-factor layer will be created by reclassing the land cover data, and finally you will multiply it all times a binary raster layer of the city. We will consider the P factor to be 1.

Build and Manage Database

Good news. You already did a bunch of data prep, and the prepared data is on the Greenfield/PROJDATAEG folder.

You’ll need the following data sets to create appropriate data for the (RUSLE) model.

· Digital Elevation Model of Syracuse - SYRDEM in the DEMs folder.

· Clipped STATSGO soils data – CLIPSTATSGO in the SOILS folder (plus a table you’ll get later from my folder)

· Clipped Land Cover 1992 NLCD – NLCDSYR in the COVER folder.

· City of Syracuse boundary – SYRUTM27 in the CENSUS folder.

Very Important!!! If you stop in the middle of this lab, make sure you save your project, but more importantly make sure you take your whole RUSLE folder from the C drive and save it to your portable media device, that way you’ll be able to pick up where you left off. SAVING A PROJECT DOES NOT SAVE THE DATA WITH IT!!!!! You do not have to save the data from my folders in your folder, but you should save the critical datasets created from the originals.

Projection, Resampling, Reclassing, and Math (CFACTOR)

1. Set up working folder(s) for this project using ArcCatalog. After ArcCatalog Opens make sure it is connected to the C drive of the computer and to my PROJDATA folder on the network. Highlight the C drive and under the file menu select new folder, and name the new folder RUSLE.

Because we want everything using the NAD27 datum and the UTM Zone 18 meters coordinate system, you’ll have to reproject the clipped land cover (it’s currently NAD83 with the USGS Albers projection). This step deals with rectification and registration…or in other words making sure our layers line-up in the same place in space (within the same spatial reference). You will also resample the land cover raster to “fit” the new spatial reference and the other raster layers.

2. From ArcToolbox find and select Project Raster from Data Management Tools\Projections and Transformations\Raster

3. For input use the NLCDSYR (did you know you can click and drag the file from the catalog tree to the input box?)

4. Make sure the output is going back into your RUSLE folder and name it LCUTM27.

5. For output coordinate system find and select the NAD 1927 UTM Zone 18N.prj in Projected Coordinate Systems\UTM\NAD 1927. Click Apply and Ok.

6. Select Nearest for resampling technique.

7. For Optional Cell Size type in 30.

8. Click Ok.

9. From here on out, I will assume that you will close the processing window when the process is done.

10. Start ArcMap, turn on the Spatial Analyst Extension (Tool menu – Extension – Spatial Analyst) and the Spatial Analyst toolbar (View menu – Toolbars-Spatial Analyst), and add all the layers that use the NAD27 datum and UTM Zone 18 meters projected coordinate system.

Question 1: What are the names of the 4 datasets that use the NAD 1927 and UTM Zone 18 coordinate system?

11. “Save As” the project into your RUSLE folder. Name it RUSLE.mxd.

12. In the table of contents move the land cover layer below the DEM.

13. Zoom in on the edge between the land cover data and DEM near their overlap. Make sure you zoom so close you can see the individual cells.

Question 2: Do the cells line up perfectly between LCUTM27 and SYRDEM? Are the cells the same size between the two data sets?

Obviously, you have to do a little more work so that they line-up exactly (a little more of that rectification, registration, and resampling). You know that expression, “trust, but verify?” Well actually, I tend not to trust much of anything in GIS. I learned the following method a long time ago by trial and error and experimentation so that I could make sure the layers lined-up. It’s a little quirky, but it works.

14. Activate ArcToolbox in ArcMap if it isn’t already there.

15. Find and select the Resample tool from Data Management Tools\Raster.

16. For input use the LCUTM27, for output, name it RSLC (is it going back into your RUSLE folder? – This is another thing I’m going to expect that you do…KEEP TRACK OF WHERE THESE THINGS GO AND ARE SAVED BEFORE YOU CLICK AWAY!!!!!!! ALWAYS CHECK THE PATH!!!!!)

17. Select Nearest for the resampling technique. Look at the help window to answer the question below.

Question 3: What are the other two resampling techniques? What does NEAREST refer to? (copy what it says in the help window)

Question 4: Okay think about what you are doing for a minute. You are about to take a raster data set that has an original resolution of 30 meters and turn it into a raster data set that has a resolution of 10 meters. How many 10m by 10m cells would make up a 30m by 30m resolution cell?

Question 5: Ignore the spatial reference change for a minute. Compared to the original 30m land cover dataset, will you have higher, lower, or the same accuracy with the resampled land cover data?

18. Click on Environments. Click on the general settings link and change the output extent to SYRDEM (use the browse button to select the data set and the window will change to Same as SYRDEM.) If this doesn’t work you may use the SYRUTM27 layer….that’s what we used to clip SYRDEM anyway.

19. Wait until the coordinates appear below, then change (drop down) the output extent to “As specified below.” The coordinates should not change; if they do, repeat the step above.

20. Select SYRDEM for the “Snap Raster.” (use the browse button to select SYRDEM and the window will change to Snap to SYRDEM)

21. Now use the Raster Analysis Settings link and change the Cell Size the same as the SYRDEM (the window will change like above and the window below will change to 10 meters).

22. Click ok to close Environment Settings.

23. If the optional cell size still says 30 change it to same as SYRDEM as you did above. (No, I don’t know why some things refresh and others don’t… but again I don’t trust the software).

24. Click ok to run the tool.

25. Take a look at this layer vs. SYRDEM now. I know it’s hard to see, but they should line-up perfectly.

Now, you’re going to change the existing values of the cells in the land cover data into the C-factor values. This will have a few steps so be deliberate. Because we cannot just reclass integers to the real numbers of C-factor, you will have to first reclass to integers, convert the dataset from integer data to float data (real numbers) and finally divide those numbers by 1000.

There are two places to access reclassify. There is a tool in the toolbox, but for this exercise, use the reclassify command from the Spatial Analyst Toolbar.

26. From the Spatial Analyst Toolbar find and click on the button with “Spatial Analyst” on it and it will drop down, select Reclassify.

27. The input raster is RSLC. The reclass field is Value. Name the output C1000 with the RUSLE folder (you checked that path, right?).

28. Click on the “Unique” button to change the old values to the NLCD values (11 through 92).

29. Use the following table to enter the numbers in the table in the reclassify window. You will have to click your mouse within the New Value column to change the number.

Old Value	New Value
11	0
21	50
22	30
23	60
32	600
41	6
42	4
43	5
81	100
82	220
85	100
91	2
92	20

30. Click ok to reclassify.

31. Now find and select the Float tool from ArcToolbox within Spatial Analyst Tools\Math. What you are about to do is to convert this integer data to float data.

32. The input is C1000, and name the output C1000F (CHECK FOLDER!!). Click Ok.

33. Now find and select the Divide tool from ArcToolbox within Spatial Analyst Tools\Math. Now you’re going to create the actual CFACTOR layer.

34. The input 1 is C1000F, and name the output CFACTOR (folder?). For the second input type in 1000, and click ok.

Question 6: What is the resolution of CFACTOR? What was the resolution of the original land cover data you used to create CFACTOR? (This will be important later) (Don’t forget units)

35. At this point you may delete RSLC, C1000, and C1000F from your RUSLE folder. You don’t have to, but if you need to save space on your portable storage you may want to. Deleting data sets can sometimes be a little difficult. First, make sure these layers are “removed” from ArcMap (right click on the layer and select Remove). Second, go to ArcCatalog and right click on the layer name and select delete.

Great! Two parts of the equation done, and a few more to go.

Table Calculations, Joining, and Rasterization (KFACTOR)

Now, you’re going to play with the tables in the vector (feature classes) data for soils so that when we convert it to raster (rasterize), it will have the real number values we want. This is how you’ll create the KFACTOR layer.

1. Turn the Editor toolbar on (View menu to Toolbars to Editor, check marks mean that it’s “on”)

2. Add LAYER.dbf to ArcMap from the SOILS folder in my PROJDATA folder. It is buried deep in the path JCC\Census90andStatgo\STATSGO_Statewide\NY\SPATIAL

3. If you cannot see the table you just added, click on the Source tab on the table of contents (bottom).

4. Open the LAYER table by right-clicking on the layer in the table of contents and select open with the table icon next to it. We’re going to want data from the field KFFACT, but it is based on the soil layer (horizon), so we’ll cheat a little and use an average of KFFACT.

5. Highlight (click on) the field for MUID (It should turn blue).

You will use MUID because this field is common to both the Layer table and to the clipped STATSGO data. We need a field of common attributes to put separate tables in our database together. This concept is a database concept, but it is especially important in GIS where we link tables with no spatial information to tables with spatial reference (the mapped features).

6. Right click while you’re on the field heading and select Summarize. The Summarize window will open.

7. The first window (1.) should show MUID, if it doesn’t select it with the drop down menu.

8. For (2.) scroll down to KFFACT (NOT KFACT), click on the + button to see the statistical choices, select Average, by clicking the box (should see checkmark).

9. Name the table sumsoil.dbf, and make sure your output table is going where you want it to go so we can find it later (RUSLE folder!). Say Yes to add it to the map. Open the table and take a look.

10. Now you’re going to join the summary table to the soils attribute table by right clicking on CLIPSTATSGO (the soils data layer), selecting Joins and Relates, and then Join.

11. In the window number 1, the field you’ll use is MUID.

12. For the number 2 window, select the summary table (sumsoil.dbf).

13. For window number 3 the field you’ll use is again MUID.

14. Click Advanced and select the radio button (option) for Keep Only Matching Records. Click Ok. Click Ok.

15. Open the attribute table for the soils data (CLIPSTATSGO), and you should see your new Average_KFFACT field.

Again, the idea here is the database concept that if we have common fields (in this case MUID) or keys we can combine data from lots of other tables to our spatially referenced polygons. Next, you have to convert the soils polygons into a raster layer.

16. Find and select the Feature to Raster tool from ArcToolbox Conversion Tools\To Raster.

17. Use the drop down button (and not the Browse button) next to input feature to select the layer CLIPSTATSGO; for field, find and select sumsoils.Ave_KFFACT; name the output raster KFACTOR…okay I should probably let go of reminding you about making sure that this data set is going into the correct folder, but to be honest, it’s a mistake I still make. Don’t let the machine default on you; you determine your own path!!!

18. Now click on Environments and do the output extent, snap raster, and cell size like you did for the land cover resampling above (not the reproject). Click ok to exit the settings and make sure the cell size is also set to same as SYRDEM and click Ok.

Question 7: What is the resolution of KFACTOR (Don’t forget units)?

Wow, so let’s see…. you have R, K, and C. Now, do the LS or a.k.a “slope” factor.

Derivative Mapping, Hydro Tools, Calculator (LSFACTOR)

1. From the Spatial Analyst toolbar select Slope from Surface Analysis

2. The input surface is the SYRDEM, the output measurement is Degrees. Leave the Z factor and output cell size alone (the cell size should still be 10 like the original DEM), name the output raster SLOPE.

Question 8: What is the resolution of SLOPE based on the metadata? (Don’t forget units)

Question 9: Okay, now I tell you that SLOPE was calculated using the elevation values from its eight neighboring cells. Or, in other words, the value results from a “smoothed” 9-cell collection. What is the more accurate resolution of slope now that you know how it processed? (UNITS?) (Hint: three 10m resolution cells by three 10m resolution cells)

3. Find and activate the Flow Direction Tool from ArcToolbox (Spatial Analyst Tools\Hydrology)

4. The input surface raster is SYRDEM, the output is FLOWDIR (always make sure that these data sets are going into your working folder). Leave the other options empty. (If the edge option is turned on (w/check), turn it off). Click ok.

Question 10: Same old question here, but this time we’ll combine them. What is the resolution of FLOWDIR according the metadata? (and units) But, it used the neighborhood (the 8 surrounding cells) like slope. So what is the more accurate resolution based on your answer to 9 above? (units?)

5. Find and activate the Flow Accumulation Tool (near Flow Direction).The input surface raster is FLOWDIR, the output is FLOWAC. Don’t do the weight raster option. Click OK.

For the following step, be prepared; the raster calculator is VERY fussy.

6. Open the Raster Calculator from the Spatial Analyst toolbar (If you can’t see all the functions, make sure you can see all the operations like trig by clicking on the button in the lower right-hand corner with the >>) and enter the following expression…CAREFULLY!!! Then press Evaluate.

Pow([FLOWAC] * 10 / 22.1, 0.6) * Pow(Sin([SLOPE] * 0.01745) / 0.09, 1.3) * 1.6

7. The new dataset CALCULATION should appear in ArcMap’s TOC, right click on the layer and choose Make Permanent, name it LSFACTOR and save it to the RUSLE folder. Change the name of the layer in the table of contents (TOC), by clicking twice on the name and type in LSFACTOR over CALCULATION, and press Enter.

8. You may delete FLOWDIR, FLOWAC, and SLOPE from your RUSLE folder.

Rasterization (Last one - SYRBND)

You’re going to rasterize the boundary file of Syracuse to create a binary layer that will limit the soil loss equation results to the specific city boundary.

1. Open the attribute table for SYRUTM27 (Remember? Right click on the layer name in the TOC and select that option).

2. On the table window, click on Options and select Add Field. Name the field BINVAL, choose Short Integer for Type, type in 4 for Precision, and click Ok.

3. After the new field appears in your table, right click on the heading (BINVAL) and select Calculate. Read the warning (no biggie, if you have to correct an error, you can always delete the field and do it again). Click Yes.

4. When the Field Calculator window opens, you’ll note that the field name is already there with the = next to it. In the window below BINVAL = type in 1, and click ok. After the 1 shows up in the new field of your table. Close the table window.

5. Select the Feature to Raster tool. Click and drag or use the drop down to select SYRUTM27 as your input. The field is BINVAL. The output is SYRBND. Set the Environmental Settings as you did with the soils data set. (If you have trouble with the extent, you may use the SYRUTM27 layer because that’s what you used to originally clip the SYRDEM). Make sure that that the output resolution is the same as SYRDEM or 10. Proceed as you have before with this tool.

Awesome! You’ve got them all. Now comes the easy part, the actual calculation.

Analysis

The following calculation will produce a value for each 10m by 10m cell. The units of this value are Tons/Acre/Year because the R value and K value are both in English units, while the LS and C factors are dimensionless. That’s basically what you are about to do, this calculation will be done for each cell, that’s why we have to be sure that our layers are perfectly lined-up in the same place in space.

Raster Calculator

1. Open the raster calculator from the Spatial Analyst toolbar, multiply, and evaluate the following datasets as an expression:

80 * [KFACTOR] * [LSFACTOR] * [CFACTOR] * [SYRBND]

2. Make the CALCULATION permanent, Name it ETPAPY. Your layer should look like the LS Factor only clipped for the city of Syracuse boundary. (To see this, you may have to turn off the other layers in the TOC. Remember if it has a check, it’s on, and if it doesn’t, it’s off.

3. Change the name of the layer in the TOC.

4. You may delete KFACTOR, LSFACTOR, CFACTOR, and SYRBND from your RUSLE folder.

Reclass

Now you’re going to do a reclass to highlight the areas of relative priority.

The city is going to be concerned about any area with a potential rate of 1 ton/acre/year soil loss, but those areas with a rate of 5 ton/acre/year is of even more concern. You’re going to reclass these areas of priority as 1 low, 2 medium, and 3 high.

1. Find and activate Reclassify from the Spatial Analyst toolbar.
2. The input raster is ETPAPY, the reclass field is value.
3. Click on the Classify button. In the classification window change classes to 3 (upper right), and click ok.
4. Change old values and new values according to the table below.

Old Values	New Values
0 – 1	1
1 – 5	2
5 - 20000	3
No Data	No Data

4. Name the output raster PRIORITY

Question 11: Is the ranking that you just did ordinal or interval data?

Congrats. You have done a pretty hefty analysis. Now on to the output requested by Lizzy.

Question 12: What is the resolution of PRIORITY according to the metadata? What is the more accurate resolution of PRIORITY given the combination of 30m and 10m resolution datasets? (Don’t forget units)

Output

Compose a memo to Elizabeth Nifkin that reports your work and findings. Remember to let your table and map do most of the “talking.” No more than one page!

Create a table on a separate piece of paper that provides the area (in Hectares) for each priority (Low, Medium, and High) and the total tons/acre/year for each of those priority zones. Your table should have 3 fields and 4 rows (including the header). Use rounded values (to the 1’s) and no decimals for the erosion and hectare values.

1. From ArcToolbox find Zonal Statistics As Table tool (in the Spatial Analysts Tools\Zonal.

2. The input feature or raster zone data is PRIORITY, the zone field is Value, the input value raster is ETPAPY. Name the output table EbyPRIORITY.dbf.

3. Take a look. The sum is the total tons/acre/year for each priority area, the count (cell count) is our area for each priority zone. Huh? Well, each cell is 10 by 10m so count * 100 is area in square meters. For conversion, 10,000 square meters = 1 hectare.

4. Use EbyPRIORITY.dbf to create the table requested using Word or Excel. Title the table something short but descriptive of the information….don’t forget to reference this table in your memo with its title). Please put your name and the date on it too.

Compose a map for a black and white printer showing the priority areas

1. Remove all layers except PRIORITY by right clicking on them and selecting remove.

2. Center the PRIORITY layer in the ArcMap map view by zooming or by right-clicking on the layer and selecting zoom to layer.

3. Symbolize the PRIORITY layer by right clicking and selecting Properties and then the Symbology tab. Change the colors in the boxes with color by double clicking in the box under Symbol. Use the gray colors, a dark one (but not black) for priority 3, a lighter one for 2, and the lightest one or even white for priority 1. Now under label change by typing over the number so that 1 is Low, 2 Medium, and 3 High.

4. Turn on the Layout Toolbar. (View menu-Toolbars-Layout)

5. Now change the right hand tile of the ArcMap window from the Data View to the Layout View (look for the tiny button in the lower left corner of the map view that looks like a piece of paper. You will see what it will look like on a piece of paper with a neatline already drawn around it.

6. From the Insert menu select Legend and a wizard will open (make sure you pay attention to the different options as you go through the legend). In the first window make sure your map layer (priority) are also in the legend items. If not move them over using the arrows. Click next. Change the legend title to Erosion Priority and click next. Don’t change anything on the next window, and click Next. Don’t change anything on the next window and click next. Don’t change anything on this window and click Finish.

7. Now use the Insert menu to add a title, a scale-bar, north arrow, and a text box.

8. Go to the file menu and make sure the print settings are set for the printer you will use.

9. Move the map elements around to make them look good. Double click in the title box and provide a short but descriptive title (use this as reference in your memo). Double click in the text box and put your name, date, and attribute the data to USGS, NYSDEC, U.S. Census, and USDA (you may use the acronyms).

10. To save your map, use export from the file menu and select one of the options. I like using .emf for insertion into word documents, but sometimes it’s just easier to save it as a jpg, tiff or bmp.

11. Print it from ArcGIS, but if you have to come back you can insert the exported file from step 11 into a word document and then print it.

12. After you get the copy of the map to turn in. Play with the map elements and layout toolbar to check out all the options.

Name _____________________________________________ Lab Section ________

ERE 450/550 Lab Exercise 3 Answer Sheet – Lab due 10/20 at 5pm

Hours spent doing the lab_______________ Score_____/30

Question 1: What are the names of the 4 datasets that use the NAD 1927 and UTM Zone 18 coordinate system? (4)

Question 2: Do the cells line up perfectly between LCUTM27 and SYRDEM? Are the cells the same size between the two data sets? (2)

Question 3: What are the other two resampling techniques? What does NEAREST refer to? (copy what it says in the help window) (3)

Question 4: Okay think about what you are doing for a minute. You are about to take a raster data set that has an original resolution of 30 meters and turn it into a raster data set that has a resolution of 10 meters. How many 10m by 10m cells would make up a 30m by 30m resolution cell? (1)

Question 5: Ignore the spatial reference change for a minute. Compared to the original 30m land cover dataset, will you have higher, lower, or the same accuracy with the resampled land cover data? (1)

Question 6: What is the resolution of CFACTOR? What was the resolution of the original land cover data you used to create CFACTOR? (This will be important later) (Don’t forget units) (4)

Question 7: What is the resolution of KFACTOR (Don’t forget units)? (2)

Question 8: What is the resolution of SLOPE based on the metadata? (Don’t forget units) (2)

Question 9: Okay, now I tell you that SLOPE was calculated using the elevation values from its eight neighboring cells. Or, in other words, the value results from a “smoothed” 9-cell collection. What is the more accurate resolution of slope now that you know how it processed? (UNITS?) (Hint: three 10m resolution cells by three 10m resolution cells) (2)

Question 10: Same old question here, but this time we’ll combine them. What is the resolution of FLOWDIR according the metadata? (and units) But, it used the neighborhood (the 8 surrounding cells) like slope. So what is the more accurate resolution based on your answer to 9 above? (units?) (4)

Question 11: Is the ranking that you just did ordinal or interval data? (1)

Question 12: What is the resolution of PRIORITY according to the metadata? What is the more accurate resolution of PRIORITY given the combination of 30m and 10m resolution datasets? (Don’t forget units) (4)

Spatial Model Score ______ / 10

Fill in the blanks. Ovals are processes and rectangles are filenames.