Terrace Extraction I

Terraces are distinctive landforms composed of two morphologic elements: the sloping riser and the relatively flat tread (Leopold et al., 1964). They preserve a history of past fluvial and/or glacial outwash activity. Terrace flats are rarely steeper than ~13% slope. The treads are bounded by the riser on the downhill side (a slope breakline). On the uphill side, the flats often transition to an alluvial fan, the toeslope of another riser, or bedrock outcrops.

This lesson needs to be cleaned up.………………..
Dr. Patrick Belmont’s Lab at Utah State has built the TerEx tool for mapping terraces:
Belmont’s Lab Webpage

Fluvial terraces are abandoned floodplains. They form by lateral erosion of a valley by its stream, usually during a period of aggradation, and are incised (abandoned) when the channel narrows and a new floodplain is established at a lower elevation. Incision and abandonment is driven by one or more factors, including tectonic uplift, sea level drop, sediment or discharge changes, climate change, complex interactions between these things, and modification by humans. Terrace remnants preserved on valley walls may be flat bedrock surfaces or gravel-capped benches. A possible third landform element of terraces may be added in some cases: a gently sloping fan that spans the transition between terrace flat and the adjacent hillslope (or the next higher terrace’s riser). A stacked set of terraces is called a flight.

Kame terraces form by the deposition of sediment by meltwater streams flowing along the margin of glacial ice. Kame terrace deposits commonly interfinger with alluvial fan (tributary valleys), till, glaciolacustrine, and local hillslope deposits. Kettles are often ornament kame terraces surfaces

ArcGIS WORKFLOW – Updated November 2013

1.) Download DEM. Unzip. Set up a project folder. Put DEM there.

2.) Load DEM into ArcMap.

3.) Create a clipping “cookie cutter”. Use this to crop the DEM to the river corridor study area. Draw a polygon (Customize > Draw toolbar), convert it to a shapefile (select it, then Right-click on name of data frame > Convert Graphics to Features), use Data Management Tools > Raster > Raster Processing > Clip tool.

4.) Set Data Frame to a projected coordinate system (i.e., UTM NAD83 Zone11) or project it using the BILINEAR method. I prefer the first option, generally.

5.) Turn off Background Processing (Geoprocessing menu > Geoprocessing Options > uncheck box for Enable).

6.) Optional: Store files for transport or sharing. File > Map Document Properties > check box for Store Relative Paths.

7.) Save .mxd to project folder. Save all files you create from here on out in this folder.

8.) Create a hillshade raster. Adjust transparency.

9a.) Create a slope (percent) raster.

9b.) Optional: Majority filter on the Slope raster to aggregate and reduce noise a bit (where is Majority filter – Focal Statistics?).

10.) Reclassify slope into 2 classes (terrace flats, not terrace flats).
– Use Reclassify tool (Spatial Analyst tools > Reclass > Reclassify tool).
– Set Old Value (pixel class 0-10% slope) to New Value = 1
– Set Old Value (pixels class10-100% slope) to New Value = NoData.
– The output raster shows “terrace flats”, or potential terrace flats anyway

11.) Create a contour that is some distance above the highest terrace in your area (Spatial Analyst tools > Surface > Contour List).

12.) Optional: Crop data further to reduce extraneous data and to processing time. See Step 2.

13.) Convert “terrace” pixels raster, created in Step 10, to a point shapefile (Conversion tools > From Raster > Raster To Point).

14.) Extract elevation values from DEM to the terrace points shapefile (Spatial Analyst > Extraction > Extract Values to Points). This creates a new point shapefile with elevation values in a field called RASTERVALU. Every time you use this tool, you get a new shapefile, so keep your files organized and know which one you’re working on.

15.) Clean up attribute table. The name of the field RASTERVALU isn’t helpful. Open the attribute table and create a new field (Float format, name field Elev_m). Use Field Calculator to duplicate values in RASTERVALU in this new field (Elev_m = RASTERVALU). Delete the RASTERVALU field.

16.) You need a column for downstream distance. There are several ways to do this. One way would be to create another new field (Long Integer, name it Northing or Easting). If your river runs pretty much north-south, used Northing. If it runs generally east-west, then use Easting. If your river runs diagonally, you’re screwed. Move to Utah and open up another all-night Jell-O & yarn shop on Main Street.

17.) Export data in attribute table to Excel (.dbf).

18.) Plot Elevation vs. Downstream Distance (northing or easting or whatever you figured out). Pick a small symbol for points in the chart.

19.) Interpret terrace levels. Its likely you will see spurious data that appear as diagonal tracks*. These are flat-lying pixels in drainages or along ridgelines. Look back to the reclassified slope layer to locate these non-terrace pixels. Consider labeling them in the chart.

* It is possible to filter out the diagonal stream/ridge pixels using curvature criteria. Create a profile curvature raster, extract the values to the point table as in Step 14, determine both positive curvature (convex ridges) and negative curvature (concave gullies) value thresholds, sort the column and delete points meeting either criteria.

** There are more efficient ways to do this in R or Python. Processing speeds on the PCs we had access to at the time of posting caused us to take a conservative route – not that there’s anything wrong with being conservative. Recall that our example is for a portion of the Okanogan River, a partially glaciated valley at LGM. Kame terraces dominate the landscape below about 550m elevation. The channel runs north-south, thus we used Northing for distance downstream in our final scatter plot. The Demoulin et al. paper describes an alternative technique for fluvial terraces in a tectonically active region.

Also, if terrace flats are not sufficiently interpretable, try dividing the channel into 500m-long reaches, working through steps, and plotting. In our Okanogan Valley example (glacial terraces), we found 500m reaches worked better for delineating terraces. The many plots necessary was definitely annoying and perhaps prohibitive on longer channels. The quicker way to accomplish all of this is using program R. See ‘Terrace Extraction II’ post for those instructions.

Demoulin et al. (2007)
Leopold, Wolman & Miller (1964) Fluvial Processes in Geomorphology book
Ebert et al. (2011) Geomorphology 132


– We used 30m DEM data downloaded from National Elevation Data Seamless Server website (http://seamless.usgs.gov/). You could use 10m, 5m, or 3m data if you’ve got it. This lesson was developed for a portion of the Okanogan Valley, WA. A flight of kame terraces, each a few meters high, is well preserved there. You might want to use a finer resolution DEM if looking at thinner fluvial terraces. A finer DEM may be the only way to resolve ancient or poorly preserved terrace remnants.
– Process data for reaches no longer than 3-miles (rivermiles)
– Consider only data for river-right or river-left side of valley
– Consider only pixels with slopes =< 13% and elevations =< 550m (substitute your own values)
– Remove modern channel and floodplain pixels from consideration (3m above reach outlet elevation)
– Optional: Remove tributary channels and floodplains from consideration
– Plot relative altitude (Y) vs. downstream distance (X)
– Export image of chart to Illustrator for creation of long profile figure
– Repeat steps for additional reaches

Suggested output filenames are in parenthesis.
Raster filenames must be 13 characters or less, no spaces.
Figure out a file naming convention that keeps names short, but reasonably clear.
Always use the buttons in Raster Calculator to create scripts; avoid typing as much as possible.

1.) Add sink-filled, projected DEM to a new, blank ArcMap document. Save the .mxd to your workspace folder.

2.) Add rivermiles point shapefile that corresponds to RM marks on USGS 1:24k topographic quads. Inlet and outlet elevations, referred to below are the DEM elevations at RM points.
Find RM shapefile here: http://www.ecy.wa.gov/services/gis/data/data.htm#r

If your stream has no RM points assigned to it, you must create them yourself. Create the polyline (stream centerline), then segment it. In Arc9, there was a great tool called Divide. In Arc10, its gone. There are a few work-arounds, but my advice to do it with xToolsPro.

3.) Spatial Analyst > Surface > Contour List tool > create a single contour at 550m. Choose a contour appropriate for your study area. (Contour_550m.shp)

4.) Customize menu > Draw tools > Draw rectangle box perpendicular to river channel, between two RM points forming a 3-mile reach. Size the edges of the box to fall near the 550m contour. Orient the graphic and convert it to a shapefile via Right-click Layer > Convert graphics to features. (BoxA.shp)

5.) Clip DEM to box with Data Management tools > Raster > Raster Processing > Clip. Remove original DEM from view. (BoxA_dem)

6.) Use the Identify tool and the clipped DEM to find the elevations of inlet and outlet points (upstream and downstream RM points). Write these down. The tool is a blue, circular button on the main ribbon.

7.) Spatial Analyst > Surface > Slope tool > create a percent rise slope raster for DEM. (BoxA_slp_p)

8.) Spatial Analyst > Surface > Hillshade tool > create a hillshade layer for the DEM. (BoxA_hs)

9.) Define channel and floodplain by creating two pixel classes: Channel and Non-channel. First, add 3m  to the outlet elevation (in your head). Our outlet elevation was 248m. Second, open Raster Calculator and use buttons to create script:
“BoxA_dem” =< 251
The output raster will contain reassigned values of 1=channel, 0=non-channel. (Channel)

10.) Define terrace flats using Raster Calculator. Script looks like this:
“BoxA_slp_p” <= 13
If spurious results are generated, substitute the toggled value for 13 (in our example 1,535,400) found in Properties > Symbology > Classify button > click the % toggle button for two sets of values. Output raster will contain values of 1=flats and 0=non-flats. We used 13% as a reasonable geomorphic break value, but you may find alternative threshold works better in your area. (Map A)

11.) Define elevation range using Raster Calculator. We used 550m. Script looks like this:
“BoxA_dem” <= 550
Output contains values of 1=elevations below 550m and 0=elevations above 550m. (Map B)

12.) In Raster Calculator > Add Map A + Map B. Values in output raster will be 0,1,2. 0=neither slope nor elevation condition satisfied, 1=either slope or elevation condition satisfied, 2=both slope and elevation condition satisfied. You’ll end up keeping the 2’s. (Map C)

13.) Convert the pixel values in Map C raster to points in a shapefile using Conversion tools > Raster to Points, input = Map C. (ExtractedPts.shp)

14.) Turn on the Channel raster and ExtractedPts layer. Open the Editor toolbar > start editing the ExtractedPts layer. Use the Editor selection arrow to select and delete all points on one side of the channel, keeping all of the other side’s points. You want to just analyze just one side of the valley, so pick the side with the best terraces (or do both sides separately). When done, Stop Editing and Save Edits.

15.) Next, pare down the points further. Editor toolbar > Start editing ExtractedPts.shp > Open attribute table and isolate all records with values of 0 and 1 in the Value field. See Step #13 for definitions of 0,1, and 2. Once selected (use Select by Attributes), right-click on far left side of table > Delete Selected. When done, Stop editing and save edits.

16.) Add X and Y fields to ExtractedPts attribute table. We projected the initial DEM to NAD83_UTM_Zone 11, so our X is Easting (m) and Y is Northing (m). Do this with Table Options > Add field, float, name it X. Right-click new field heading > Calculate geometry > X Coordinate. Field should populate. Repeat for Y field, making sure you calculate the Y Coordinate.

17.) Spatial Analyst > Extraction > Extract Multiple Values to Points tool > input ExtractedPts and DEM. This will append the DEM (elevation) values to the point attribute table. Processing may take a few minutes.

18.) Add another field to the attribute table, integer, called “RelAlt”. Populate this with ‘relative altitude’ values, which are just the DEM elevations minus the outlet elevation. Do this in Field Calculator, which works similarly to Raster Calculator, but with attribute tables. Right-click on the new field > Field Calculator and enter this script:
“BoxA_dem” – 248

19.) Spatial Analyst > Extraction > Extract Multiple Values to Points tool > input ExtractedPts and Channel. This will append the Channel values to points attribute table. Processing may take a few minutes.

20.) Editor > Start editing ExtractedPts.shp > Open attribute table, Select by Attributes, and delete all records with Channel field values = 1. When done, Stop editing and Save edits.

21.) Export points attribute table to .dbf. Open this in Excel or other graphing program.

22.) Create a scatter plot with x-axis = Northing (distance downstream) and  y-axis = RelAlt. Make points small and transparent. Stacking patterns will reveal prominent terrace levels. Alluvial fans will show as gradational patterns of points rising from the terrace flats.

23.) Export image of plot (.jpg, .png) of multiple plots to Illustrator. Make sure to scale each exported image the same. Trace terrace tops along valley. Label (number) terrace levels. Add long profile of river surface. Add other geographic information to profile figure (tributary confluences, major road/RR crossings, ridgeline profile, geologic contacts-below long profile, etc.).