Sunday, June 10, 2012

GIS Terrace Extraction

I post here quite a bit about my thesis research into the Northern Goshawk breeding ecology. However, I am also studying to complete my graduate certificate in Geographical Information Analysis. Shhhh, don't tell anyone! A wildlife biologist friend of mine once said that GIS knowledge is absolutely critical to be a good biologist, but not to tell anyone you know how to do it or that is all that you will ever do! Well, I will break the rules a bit here.

This last semester our GIS class was working to reproduce a process described in a published paper to analyze a river corridor and extract where the terraces are located within the canyon. While my focus of GIS is on ecological applications, this particular course was very heavy in geology.

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

As is often the case, some of the key details to complete the process were left out of the paper. The class struggled with a few steps and ultimately we simplified the process down to what is described in Terrace Extraction I. However, I continued with the original assignment and used a bit of automation to complete the process. I co-authored this post for my professor's web site: GIS 4 Geomorphology. I am reposting it here , with a few small modifications for presentation and consistency with this blog.

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.

This case study is an example 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 developed a new method for glacial outwash terraces in Washington state.

One challenge was the fact that ArcGIS required approximately 20 steps for each reach (500m long). With a 70km river to analyze, the result is several thousand steps, each with a high risk of operator error. While the following procedure will require some modification to accommodate other geographic areas, the general concept and the code that implements it should be readily adaptable.

We used ArcGIS 10 and the free, open-source program R ( 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 Digital Elevation Model (DEM) data (XYZ data).

1.) Download 30m Digital Elevation Model from the USGS Seamless Server ( or another source.

2.) Load DEM into ArcGIS.

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

4.) Use the Draw tools to create a polygon encompassing the center of the river and far enough up the canyon wall to be above any possible terraces. This should include the whole length of the river to analyze. This can be repeated to analyze the other side of the river.

5.) Convert this graphic to a shapefile (right-click on Layers > Convert graphics to features)

6.) Clip the DEM with this newly created 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 the rectangular polygon extent.

7.) Export this newly created clipped DEM to a GEOTIFF file (right-click DEM > Data > Export Data).

8.) Record the rectangular extent of this raster for use in customizing the R script listed below (Lat/Lon in DD units). You are now done in ArcGIS!

Clipped DEM for river Right.

The simplifying factor in our case is that the river runs north-south. A river running east-west would work equally as well. If the river were to run diagonally, a more advanced algorithm for the moving clip window would be required. With N-S or E-W, only one dimension of the moving clip window is required and a rectangular clip region can be used. You can always rotate the DEM in ArcGIS prior to bringing the data into R.

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.

# 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.
  ext <- extent(-119.76, -119.375, 48.98 - i*0.01, 48.985 - i*0.01)
  # Clip the raster to the just established extent and keep in tmpras
  # the original ras object remains unaltered for the next pass thru.
  tmpras <- crop(ras, ext)
  # extract the minimum elevation from this portion of the DEM
  rmin <- minValue(tmpras)
  # Subtract the minimum elevation from all points in raster.
  # this time tmpras is changed.
  tmpras <- calc(tmpras, fun=function(x){x-rmin})
  # convert relative elevation values less than 5m to NA or NoData
  # This is to omit the flood plane from being identified as a terrace
  tmpras <- calc(tmpras, fun=function(x){x[x<5] <- NA; return(x)})
  # Now plot the slope versus the relative elevation
  # X-axis - round(tmpras,0) # integer value from the DEM
  # Y-axis - round(slope(tmpras, latlon=T)/27.7, 0) # calculate the slope of DEM
  # but since it is not projected, we have to divide by the cell size in meters
  # then round to an integer. Rounding isn't required for the process, but it looks better
  # xlab and ylab define the labels
  # xlim and ylim define the sclae of the axis - clips everything outside
  # cex sets the point size, pch chooses the point symbol.
  # main defines the title for each chart, paste combines strings of text. In this
  # case the chart number and the extent covered.
  plot(round(tmpras,0), round(slope(tmpras, latlon=T)/27.7, 0), xlab="Elevation(m)", ylab="Slope", 
       ylim=c(0,30), xlim=c(0,200), cex=0.5, pch=20,
       main=paste("Terrace Map #", i, " X=",ext@xmin,":",ext@xmax,", "," Y=",ext@ymin,":",ext@ymax, sep=""))
  # complete for this clip window. increment at for loop and begin again
# If opened a pdf file, then must close it before viewing. Don't call
# if using a plot window as it will close the window.

Created by Pretty R at

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 choosing a narrower section of river (500m versus 1km). 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., B. Bovy, G. Rixhon, and Y. Cornet. 2007. An automated method to extract fluvial terraces from digital elevation models: The Vesdre valley, a case study in eastern Belgium. Geomorphology 91:51-64.

No comments: