GLG410--Computers in Earth and Space Exploration

Announcements Syllabus Schedule Weekly lecture notes Assignments Links

Combining analysis of rasters and vectors: application to landscape reconstruction

You have now learned the basics of creating and working with both rasters and vectors in ArcMap. In this exercise, I want to give you a chance to work with some vector and raster data to get experience converting between them and then using the raster calculator to query the raster to extract the parts that satisfy your constraints.

To learn about this, we will apply the processing to reconstructing a piece of landscape that was once more continuous, but a river system incised. We might want to know how much material was eroded in an area, and how the erosion varied over the area.
First, a little background Many rivers are bounded by their former floodplains which have been abandoned as the river has cut down due to climatically, tectonically, or internally driven processes. These former floodplains become relatively flat surfaces called "terraces" which are separated from eachother by step slopes called "risers." In the case of the Salt River here in the Valley of the Sun, the terraces were studied extensively by former ASU professor Troy L. Péwé and his students. Jessica Block worked on the terraces as part of her M.S. thesis, finishing in 2007.
Have a look at figures 6 and 7 in this field trip guide: Urbanization, landscape, and geologic history along the Salt River (24.5 Mb pdf file) .

In the examples for our class, we are going to reconstruct the Mesa terrace landscape and compute how much erosion has occurred since that time.

Video explanation:

Caveat with respect to the documentation on this page: ArcMap 9.3 and ArcMap 10 has important changes

Note that most of the text and graphical documentation in this page is derived from ArcMap 9.3 whereas the videos are all produced using ArcMap 10. There are important differences!

Step 1: gather the data

I have organized everything for you here: Salt_Verde.zip. Just download it, uncompress, launch ArcMap, and then open the Salt_Verde.mxd project file.

Step 2: Rectify the terrace map and then digitize the terrace outlines as a polygon shapefile

For the demonstrations of today's lecture, I have done this for you. You would follow Georeferencing raster data in ArcMap to do the rectification and Digitizing the geologic map in ArcMap to do it yourself. Video explanation:

Step 3: Extract the portions of the elevation data which are underlain by the Mesa Terrace

Do this from in ArcToolbox: Spatial Analyst->Extraction->Extract by Mask.
Specify the input surface (original DEM), the mask (the shapefile), and the output raster.
You will get a raster that has the elevation values under the Mesa terrace remnants and then NoData in the other cells.
Video explanation:

Step 4: Convert the raster into a point shapefile with grid node elevations as attributes

Unfortunately, Arc cannot convert a floating point grid to a shapefile with the grid node values as attributes. As a work around, we can make the grid into an integer grid. And, we can multiply the elevations by 1000, so it is in millimeters above sea level. This will maintain sufficient precision. We just have to remember later to convert it back to meters by dividing through by 1000. The result will be a new grid called Calculation.

First convert to an integer grid multiplied by 1000 using the Raster Calculator in Spatial Analyst:
Now, Convert raster to Features using 3D Analyst:
You want to make sure the Input Raster is the integer grid multiplied by 1000 which you just made. The Output geometry type should be Point. You can change the output feature name to something more meaningful than the default.
Now you have a shapefile with a point at each grid node and the attributes are those elevations in millimeters as integers...
Note that at the bottom of the lecture is a more complicated but slightly interesting alternative for the last step.
ArcMap 10 Video explanation. Note that things are a bit different:

Step 5: Extrapolate between the remnants to reconstruct the surface

This is the trickiest part. There might be a better way to do this, but this is how I do it.
  1. Make a triangular irregular network across all of those points (which were the nodes of the DEM). Basically this is a linear interpolation and so it assumes that the landscape was relatively flat in that area prior to erosion. Do this in 3D Analyst. Note that you need choose the correct shapefile and make sure you change the height source to GRID_CODE.

  2. Convert the resulting TIN to a Raster (Convert Feature to Raster). Do this in 3D Analyst. Note that I changed the output grid size to 10 m.


    3. ArcMap 10 Video explanation. Note that things are a bit different:

Step 6: Subtract the current topography from the reconstructed topography to get the amount of erosion

This is the fun part.
Don't forget that the elevations are still in millimeters, so divide the resulting grid by 1000: [mtingrid] / 1000

Note that sometimes Arc won't be happy with all of this, and you might have to build each step separately, making the resulting grids permanent and each time launching a new ArcMap project.
Make the resulting calculation permanent by right clicking on it and choosing Data->Make Permanent. Make sure you add it into the Project.

Using raster calculator, simply subtract the current topography from the reconstructed topography (end result of extrapolation part above). The resulting map is the erosion at each point since Mesa Terrace time. Do this in the Raster Calculator: [mtingrid_m] - [slt_vrde_ned]. Note that I made the file permanent and called it erosion

You will find that the result has a lot of negative values which are places where the current topography is above our interpolated surface and so not reasonable. In the raster calculator, you want to use the Conditional command. It has the structure of con(if true, do this). Look at the ESRI explanation: CON. This is a very powerful capability of Raster Calculator. Use the Raster Calculator once more to just give you the positive erosion: con([erosion] >= 0, [erosion]): Everywhere that [erosion] >= 0 is true, give back the value of the erosion grid and make the other values NoDATA.

Here is the result: Salt_Verde_erode.pdf
If the Mesa Terrace is 400,000 years old, compute the erosion rate by dividing the erosion by the years (again in raster calculator): [poserode] / 400000

Here is the result: Salt_Verde_erode_rate.pdf

ArcMap 10 Video explanation to compute the erosion. Note that things are a bit different from ArcMap 9.3:


ArcMap 10 Video explanation to clean up the erosion map:


ArcMap 10 Video explanation to make an erosion map:


ArcMap 10 Video explanation to format the erosion rate map:

Assignment 15: Erosion Magnitude and Rate for Mesa Terrace near ASU

Perform a similar analysis as above to determine the erosion magnitude and rate for the Mesa Terrace around ASU based on this geologic map: ASUPewe.jpg. Note that the Mesa Terrace is the light green unit at the edges of the northeast and south and is designated Qefm. You will need to download relevant Digital Raster Graphics and possibly aerial photography for georeferencing and the DEMs for the actual landscape reconstruction from the USGS National Map.
Video explanation:

Note that you will use the USGS National Map instead of the now deprecated USGS Seamless site.

Make the following maps:
  1. Mapped terrace shapefile over rectified geologic map.
  2. Erosion magnitude over shaded relief
  3. Erosion rate over shaded relief
How do the magnitudes of erosion and erosion rate vary between the Salt Verde confluence (our class example), and the area around ASU? Why do you think they are different? At these rates, how deep will the river be in another million years?

Grading Rubric (50 points)

Assignment is due by the beginning of class, Monday, November 30.
Last modified: November 17, 2015