Terrace Extraction II

Robert Miller   Raptor Biologist, Boise State University
Skye Cooley   G4G Editor

GIS projects with a large number of repetitive steps are ripe for automation. Automation requires increased upfront investment in learning and time, but results are often more reliable and highly reproducible.

Dr. Patrick Belmont’s Lab at Utah State has built the TerEx tool for mapping terraces:
Belmont’s Lab Webpage

This case study is anexample of automated terrace extraction. The study area is a 70km reach of the Okanogan River valley extending south from the US-Canada border to Omak, WA. The procedure was inspired by a paper by Demoulin et al. (2007). They developed a method of teasing out subtle fluvial terrace remnants intermittently preserved along a narrow gorge in Belgium. We expand on their work with the development of a new semi-automated method for delineating glacial outwash terraces associated with the Okanogan Lobe of the Cordilleran Ice Sheet.

We use 500m-long reaches as our analysis view. The 20+ steps required to process each reach soon became prohibitive in ArcGIS. With a 70km river to analyze, the result is several thousand steps, each with risk of operator error. While the following procedure will require some modification to accommodate other river valleys (orientation, width), the general concept and the code that implements it should work.

We used ArcGIS 10.x and the free, open-source program R (www.r-project.org). R is a powerful, flexible program with tools for both vector and raster data. R’s power in this application is its elegant handling of DEM data (XYZ data). See the lesson R Help for more.

1.) Download 30m Digital Elevation Model from the USGS Earth Explorer (earthexporer.usgs.gov) or another source (see DEM Data Sources).

2.) Load DEM into ArcMap.

3.) Do not project the DEM yet. The R script assumes no projection.

4.) Use the Draw tools to create a polygon encompassing the entire river corridor study area. Make your polygon extend wide enough to encompass the highest terrace level. This can be repeated to analyze the other side of the river. Buffer tool (on river centerline) may offer a cleaner result than a drawn polygon.

5.) Convert the drawn graphic to a shapefile (right-click on Layers > Convert graphics to features). If you used the Buffer tool, skip to next step.

6.) Clip the DEM with corridor shapefile (Data Management Tools > Raster > Raster Processing > Clip). Fill the checkbox for “Use Input Features for Clipping Geometry”, which will clip the raster to the polygon and not to its rectangular extent.

7.) Turn off the original DEM. Export the new, clipped DEM as GEOTIFF or GRID to your workspace folder (right-click DEM > Data > Export Data).

8.) Record the rectangular extents (4 corner coordinates, Lat/Lon in DD units) of this raster for later use: customizing the R script below. You are with ArcGIS.
A simplifying factor in our study area is the Okanogan River runs almost due south. A river running east or west would work equally as well. If your river runs diagonally, however, you will have to modify the coordinate-shift portion of the R code so the moving analysis window (500m reaches) shifts in Lat and Lon directions. For the Okanogan, our window just steps south (no Lon shift). A work around for diagonal rivers is to rotate the DEM in ArcGIS prior to bringing the data into R. This messes up your coordinates; you may have to use a custom coordinate grid.

Overview of R Workflow: The high level process is outlined in the lettered steps below. I’ve added many instructional comments in the code itself (comments follow # symbols):

a.) Read in the clipped DEM from above as a GEOTIFF file.

b.) Create a pdf file to capture all of the graph output (two graphs per page, one per moving clip window). This example produces 56 graphs!

c.) Determine the size of the moving clip window and how many clips exist in the valley (56 in this case). Clip windows should be in the range of 500-1000m of river, but can extend well beyond the edge of the map in the other direction (full extent width in my case). Greater than 1000m produces too much noise and crisp terrace definitions cannot be seen. The example below uses approximately 1000m, but a few of the sections could benefit from a smaller region.

d.) For each clip region, clip the raster.

e.) Find the minimum elevation.

f.) Subtract this elevation from all points to produce a relative elevation raster (normalizes all of the graphs).

g.) Remove data points for the lowest 5 meters of relative elevation (you don’t want to count the modern floodplain as a terrace).

h.) Plot the integer slope versus the integer relative elevation.

i.) Analyze each graph for local minima to identify terraces per the Demoulin et al. paper.

Data processing now moves into R and is fully automated within the script below. Here is the R code which implements this procedure for our study area. Pretty R Syntax Highlighter can be helpful (http://www.inside-r.org/pretty-r/tool). Copy/Paste this code there.

# Load the various spatial libraries that I will use.
library(raster) # key routines for reading in and clipping raster files (DEMs)
library(SDMTools)  # Includes routines for Slope and Aspect

# set the working directory on my local machine.
# All data files are relative to this location

# Read in raster datafile from above mentioned directory
ras < - raster("newRRFinal1.tif")

# Either run the windows() command to open a plot window, or save to a pdf file.
# windows(record=T) # opens plot window. record=T enables scrolling thru previous graphs
# instead let's save to pdf
pdf('RiverRight.pdf', width=7.5, height=10, paper="letter", pointsize=11)

# arrange plots for two plots per page (two rows, one column)

# For loop repeats 56 times for this study area. Each time "i" is incremented
for (i in 1:56){
  # Setup moving clip window. Easting is constant for this example,
  # but northing moves north to south by 0.01 degrees for each iteration (~1km).
  # remember the "i" is incremented each time thru.

Okay, now we have a 56 page PDF file with these odd looking graphs of Slope versus Relative Elevation. If you refer back to Demoulin et al. (Fig. 6), you will notice that terraces are defined by points where the slope is zero surrounded by points where the slope is very different from zero. Here is a good example (see Terrace Map #7). Terraces in the figure below are located at approximately 42m, 98m, and 165m above the river.

Here is an example that did not work out so well (see Terrace Map #19). Some terraces can be seen at approximately 80m, but what is going on at 110m or 160m? This may be cleaned up by chosing a narrower section of river. This can also be an artifact of side canyons and tributaries.

There are other methods for identification of terraces (see ‘Terrace Extraction I’ post).

Demoulin, A.; Bovy, B.; Rixhon, G.; Cornet, Y., 2007, An automated method to extract fluvial terraces from digital elevation models: The Vesdre valley, a case study in eastern Belgium, Geomorphology, v. 91, p. 51-64

Ebert et al. (2011) Geomorphology 132

Straumann, R.K.; Purves, R.S., 2008, Delineation of valleys and valley floors, University of Zurich Open Repository and Archive

See more from Robert at blog.raptorrob.com.