Terrace Long Profiles from Lidar w/ QGIS...in progress

Goal of Project

Delineate the flat tops (treads) of glacial terraces along a portion of the Okanogan Valley in north-central Washington. GIS work involves a series of steps that begin with a lidar bare-earth raster (DTM) and ends with a chart showing the longitudinal profiles of the surfaces (Elevation vs. Downstream Distance). This example is for a portion of the Okanogan Valley, WA. The river runs north-south, which makes things much easier than one that runs diagonally. We can use latitude for downstream distance in either case, but its easier on orthogonal valleys.

Software needed

- QGIS 3.26 Buenes Aires (freeware)

- LibreCalc (freeware)

- Adobe Illustrator CC

General set up

- Use the default coordsys set when adding the lidar DTM as first layer (or change as you see fit)

- Output files for this example are saved to "tempgis" folder on my drive

Steps summary

download dtm raster tiles (.tif)

merge dtm tiles (temp .tif)

remove no data values (permanent .tif)

draw rivermile boxes (shapefile)

dtm --> slp (use model)

slp -> re (export)

re --> pts (pixels to pts)

pts --> elv (sample raster)

elv --> xy (add x/y...)

manual edits --> RR/RL

xy --> csv (export...)

csv --> odf

odf --> chart

chart --> pdf

pdf --> ai

ai --> poster, figures for article, or presentation images

Detailed Steps (draft)

1.) Download the DTM tile and HS tile for a portion of the Okanogan Valley from the Washington Lidar Portal

- "okanogan_2015_dtm_12" and "okanogan_2015_dtm_12_hs"

- The DTM raster is for processing. The Hillshade raster is for display only.

One lidar tile (okanogan_2015_dtm_12.tif) for a portion of the Okanogan Valley (Monse-Mallott area) displayed two ways. The pixel size is 1 foot x 1 foot with elevatioin values in 10ths of a foot. That's a lot of pixels.

2.) Open New Project in QGIS

- Store downloaded layers in "TerraceGIS" folder or a folder of your choosing

- Load DTMs covering river corridor (project area) into map

- The default coordinate system (CRS) for my lidar data is NAD83_HARN_StatePlane..something, something; Leave it alone unless you have a reason to change it (i.e., Warp Project tool)

3.) Merge DTM tiles

- Raster menu > Miscellaneous > Merge...

- Select tiles to merge (check boxes)

- Output as temp file; resulting raster will be named "Merged" by default

4.) Remove no data pixels

- Rigtn-click on "Merged" > Export...

- Name file and save it your filespace (permanent file)

- Scroll down to No Data section

- Add a row and enter -1 to 0

- Run and add to map

- Properties > Symbology > Render type > Hillshade

- Save project

5.) Optional: Filter elevations (modern floodplain, low-elevatioin alluvial fans, freeways, railways, etc.)

- Use Export > No Data... to exclude elevations

- I excluded the broad floodplain of the modern Okanogan River because it contains a lot of points representing the low terraces, recent alluvial fans, and most of the major roads. These points needlessly slow processing. Glacial terraces, the focus of the project, stand well above the valley floor. The profile of the Okanogan River is better shown with a simple line.

6.) Display rivermiles

- Display Rivermiles points with their mile number labels; I downloaded this point file from the Washington State Ecology website

7.) Create a new shapefile containing rivermile rectangles

- I ended up making both river-right (RR) and river-left (RL) boxes for each rivermile interval in order to reduce processing time and because glacial terraces were not always well expressed on both sides of the valley or because the DTM did not extend far enough away from the river to capture them; I used RR or RL data for each reach in my final terrace profile compilation

- Tip: If making RR and RL boxes, exclude as much of the river surface as possible (keep box edges close to shoreline)

- View > Toolbars > Turn on "Shape digitizing toolbar"

- New Shapefile Layer, name it "Boxes" and store in your "tempgis" folder, Type = polygon, use the default CRS

- Toggle editing button on

- Add polygon button

- Rectangle from extent button

- Draw as many rectangles as needed to cover the DTM, being careful to locate tops and bottoms at rivermile intervals; Provide consecutive numbers for each as you go; Use keyboard arrows to move view

- When complete, save edits when done, toggle editing off, save project

Rectangles (boxes) will be used as cookie cutter clippers in the model in a later step. Their relatively small size helps reduce processing time. The ideal size for my processing patience turned out to be boxes spanning no more than ~2 km downstream, or about one rivermile. In fact, I ended up making both river-right and river-left boxes (twice as many as shown above). I'm running QGIS 3.26.1 Buenes Aires on a 2015 MacBook Air laptop and a 2021 MacMini (8GB, Monterey), which are not super fast, but work.

8.) Build a model to get some of the work done quickly

- This option creates a model that quickly clips the raster to the box polys and calculates slope for each box...

- Processing Toolbox > Gear button > Create new model > model window opens

- Inputs > Scroll down and drag Vector Features into window

- Description = Box polys

- Geometry type = Polygon...

- OK

Inputs tab...

- Drag Raster Layer into window

- Description = Raster Input Layer

Algorithms tab...

- Drag Clip Raster By Mask Layer from list above into window

- Input = Input Raster Layer

- Mask Layer = Box polys

- No data value = -9999 (click the little up arrow to autofill)

- OK

Algorithms tab...

- Drag Raster Terrain Analysis > Slope into window

- Input = Input Raster Layer

- Elevation Layer = drop down menu > "Clipped (mask)...."

Slope = Slope Layer

- OK

Model Properties...

- Name = Slope for Each Box

- Group = Terrace Tools

Save model...

- Save model > "Slope for Each Box" in "models" folder

Run the tool...

- Toolbox > search on Terraces... > "Slope for Each Box"

- Box polys = "Boxes"

- Click the green recycle button to reiterate

- Input raster layer = "okanogan_2015_dtm_12"

- Do not provide a file name (leave as temp)

- Run

- Processing takes a minute

..............Draft info below - Currently being revised...............

9.) Reclassify slope raster via export (0-3 degrees = flat)

- Select the first newly created slope raster

- Change the data by reclassifying, keeping only pixels with "flat" values, here defined as 0 to 3 degree slope.

- Right click layer > Export > Save As, "121re" in "tempgis" folder

- Set NoData values in window > check box for No Data Values > Click plus button to add a row > From = 3.01, To = 1000000

- Output file = "121re"

10.) Convert reclassed Slope pixels to points

- Toolbox > Raster Pixels to Points

- Raster layer = "121re"

- Slope values > Change name of filed "VALUE" to "SLOPE_D" (for slope in degrees)

- Output file = "121pts"

Points should look like this if you zoom in. No need to display the layer in the map. Point symbols are shrunk down to 0.2 pt to improve display. Some artifacts (roads) can be seen. Ideally, you would delete road points, but that would take a lot of time. Remember, most roads are located low in the valley, near the modern alluvial floodplain. Glacial terrace treads are much higher, well separated from most roads, reservoirs, and the modern channel of the Okanogan River.

7.) Extract elevation raster values to points

- The tool will append extracted elevation values to the ELEV field in the attribute table for the reclassed slope points shapefile; If this takes a long time, you know there are a lot of flat pixels (points)

- Toolbox > Sample Raster Values

- Input = "121pts"

- Raster layer = "okanogan_2015_dtm_12"

- Output column prefix = ELEV

- Output file = "121elv"

8.) Add X and Y fields to attribute table

- Toolbox > Add X/Y Fields to Layer

- Input = "121elv"

- Processing takes a couple minutes

- Output file = "121xy"

** Steps 5-8 could also be modeled (single layer input, individual outputs) **

Optional: Manually pare down points that are spurious or will not be useful

- There may be obvious groups of points that are not terrace flats (water bodies, roads, gravel pits, etc.); Or you may wish to analyze points for just one side of the valley (delete those on one side of river, keep either RR or RL points)

- Make "121xy" points layer active

- Toggle editing to On

- Select Features button

- Select and delete points manually

- Save Layer Edits button - Rename if you kept river-right or river-left points (Ex: "121xy_RR")

- Toggle editing to Off

- Save project > File menu > Save As....Replace

Optional: Filter elevations below 880'

- Right-click > Filter > ...

- 880' is a value specific to my example data; below 880' is the modern floodplain; I'm reducing processing time by eliminating the modern floodplain; I can represent the river surface in the final compilation later with rivermile elevations

Optional: Remove intermediate files from map

- Layers "121slope", "121re", "121pts", and "121elv" can be removed from map and/or deleted

9.) Export points to CSV file

- CSV, or comma-delimited file, is a generic spreadsheet format useful for moving data between different software; this is the final output file involving QGIS; .csv does not support charts or most formatting; save as another format to keep such things

- Right-click on "121xy" in Table of Contents > Export > Save Features As > Export points to CSV

- File name = "121profile.csv"

- Change export name "SAMPLE_1" to "ELEV"

- Change export name "y" to "LAT"

- Place check marks next to "ELEV" and "LAT"; Don't need "SLOPE_D" (VALUE) or "x" columns

- Uncheck box for Add Saved File to Map (its at the bottom of window)

- Output file = "121profile.csv"

You only need 2 columns to create the profile chart: ELEV and LAT. There should be just 2 check marks on the left.

10.) Save QGIS Project

- Project menu > Save As > change file type to .qgs > Replace...only save a .qgs file (not a .qgz)

- Close QGIS

11.) Spreadsheet edits

- Delete "121profile.qmd" file

- Open "121profile.csv" in LibreCalc

- Save as > .odf filetype (native format for LibreCalc); .odf will process and chart way faster than .csv

Format cells...

- Format cells: Elev --> 1 sig fig; Lat --> 5 sig figs

Arrange columns...

- Cut each column > Paste Special > Values and Formats in a new column somewhere to the right

- Arrange the columns so Column A = Lat (x-axis), Column B = Elev (y-axis)

Filter low elevation points (optional; specific to my project)...

- Highlight both columns > Sort on Elev > Ascending

- Select and delete all rows with elevations 879' and below

Create scatter chart...

- Highlight columns A and B > Chart > Scatter type > Finish

Format data series...

- When chart appears, double-click on the points

- Points will turn blue and a formatting window will open

- Symbol = triangle, Color = black, Size = 0.10 proportional > OK

- Click off to unselect points (data series)

- Click on the chart to highlight the whole chart (chart area) > File menu > Export as PDF

Elevation vs. Downstream Distance plot for "flat" pixels along a ~1-mile long portion of the Okanogan Valley.

12.) Adobe Illustrator edits

- Open PDF as a new AI document

- OK the two pop-up warnings

- File > Save As > .ai format

Speed up processing...

- Illustrator menu > Preferences > File Handling > Uncheck "export in background"

- Illustrator menu > Preferences > File Handling > Uncheck "save in background"

Remove unneeded elements...

- Use white arrow to remove invisible bounding box, chart legend, and page number

- Select all > Remove clipping mask

- Select all > F6 > grayscale, no fill, black outline, 0.25 line weight

- Select horizontal scaling lines > Group them (Command+G) > Shift-slide them to the right and out of the way

- Select all triangles > Right-click > Transform > Transform Each > Uncheck 'Preview' box > 10% H, 10% V (transform objects, scale strokes and effects, scale corners) > Enter

- Select all triangles > Unite (Pathfinder palette)

- Select all > Group everything

- File > Save .ai

Select the chart's horizontal scaling lines, Group them (Command+G) and Shift-drag them to the right as shown so you can get at the points unobstructed. Otherwise, you don't want to mess around with much here. The goal is just to remove extraneous elements, recolor (F6) all objects in grayscale color space (black outline, no fill, 0.25 line weight), rescale the triangle symbols (Transform Each), then Unite the overlapping triangles to form merged clusters (Pathfinder palette), which reduces the number of objects and file size. Since our use of AI is mainly visual, we need to reduce the thousands of overlapping triangles that cannot be seen (a useless memory suck).

Use grayscale color space, black outlines, no fills, and 0.25 line weights.

Triangles: Transform Each 10% horizontal and 10% vertical (strokes and effects).

Unite all triangles using the tool in the Pathfinder palette.

14.) Compile tiles in Illustrator & interpret

- Create a new .ai document for the long profile (compilation)

- Create a regularly-spaced set of ticks for river miles (x-axis)

- Create a vertical scale composed of lines that extend across all river miles; may have to break up into parts

- Open each tile of data; Add each tile to its own layer

- Scale each tile horizontally to match river miles ticks

- Scale each to vertical scale

- Delete extraneous stuff once each tile is set (white arrow)

- Add and format other poster elements

- Prepare for printing on large format plotter (.pdf or .jpg).

Terrace flats appear as discontinuous, gently dipping surfaces. Steeply dipping forms are gullies (ignore). Other random flat spots in the landscape appear as a diffuse scatter of points (ignore). These data are for a 1-mile reach of the Okanogan above 880' elevation. Data from both sides of the valley are shown together. An improvement would be to process separate charts for river-right and river-left and include the modern alluvial plain (points below 880'). For compilation, this single tile would be squished in the horizontal direction, but not the vertical. Red lines are my quick interp of 3 terrace levels.


Terrace Profiles from 10m DEM Data from UW...in progress

- 10m data processes much quicker than lidar data.

- Set up new, blank project in QGIS with the default, unprojected Geographic Coordinate System (EPSG 4326 - WGS 84), which will be indicated in the lower right corner of the QGIS window.

- Download two large DEMs for Eastern and Western Washington. These are the "complete Eastern Washington grid" (utm11.zip) and "mosaicked zone 10 data" (zone10.zip) provided on this webpage by UW/Harvey Greenberg: gis.ess.washington.edu/data/raster/tenmeter

- Unzip both folders and rename them something more intuitive (zone11_10m, zone10_10m); They contain .adf format data.

- Load the zone 11 DEM into map.

- Since the data are in two different UTM zones (10 and 11), they cannot be mosaicked (merged). We must put them in WGS 84, then merge them. Once merged, the map display can be changed to a projected CRS (i.e., UTM) or an unprojected CRS (i.e., WGS 84). The map display CRS is shown in lower right corner of the QGIS window. QGIS default is WGS 84 (EPSG 4326).

- Load Zone 11 DEM covering Eastern WA...

- The display CRS changes to EPSG: 26711 because the DEM is in NAD 27 / UTM Zone 11N.

- Load Zone 10 DEM covering Western WA...

- It should appear in the correct location. So, the Eastern Washington DEM is in NAD 27 / UTM Zone 11N, the Western Washington DEM is in NAD 27 / UTM Zone 10. The map display is in NAD 27 / UTM Zone 11N.

- Use Warp (reproject) tool to reproject both DEMs to EPSG 4326 / WGS 84 > Target = WGS 84, No Data = 0, accept other defaults.

- Remove original DEMs from map. Keep reprojected DEMs.

- Raster menu > Miscellaneous > Merge tool to mosaic the 2 reprojected DEMs...

- Set display CRS to WGS 84 / EPSG 4326.

- Save .qgs.


an appropriate CRS (Project Properties > CRS> Projected CRS > Universal Transverse Mercator (UTM) > NAD 83 / UTM zone 11N (EPSG 26911) > Save .qgs to a folder dedicated to the project

- The project should be in EPSG 26911 - NAD 83 / UTM zone 11N (check lower right corner of window). Zone 11 DEM should be in same (check Properties > Source). Zone 10 DEM should be in EPSG 26910 - NAD 83 / UTM zone 10N (check Properties > Source).

- Raster > Miscellaneous > Merge the 2 large DEMs

- Remove

- Download rivermile points from Washington Ecology GIS site

- Delete unneeded points

- Buffer tool on rivermile points > 5 miles, Dissolved (output permanent file) <-- Creates study area corridor(s)

- Clip DEM with riverpoints buffer with Clip Raster by Mask tool

- Create Slope raster from clipped DEM (output temp file)

- Export slope raster > No Data = 4 to 1000000 (<-- Keep pixels with slope values between 0-4 degrees); (output temp file)

- Pixels to Points tool (output temp file)

- Sample Raster Values tool; Input raster = merged DEM (<-- Appends elevations attribute table of point file)

- Add X, Y Fields to Table tool <-- Appends lat/lon values to table (output temp file)

- Toggle editing on

- Edit attribute table > Properties > Fields > Delete "slope" and "x" fields > Keep only Latitude ("y") and Elevation fields

- Save edits > Toggle editing off

- Export final point file as permanent file (keep this one, discard temp files)

- Export as .csv > Don't add to map

- Save project > Close QGIS

- Open .csv in LibreCalc > Save as .odf > Select data series > Format data = black, triangle, 0.10

- Edit data format for Lat > 5 decimal places; Format for Elev > 1 decimal place

- Arrange columns by Cut > Past Special Formulas and Values > Lat column on left, Elev column on right

- Sort on elevation > Remove any rows above a certain elevation (Ex: 650m)

- Scatter plot > X-axis = Lat, Y-axis = Elev (<-- Elevation vs. Downstream Distance for my south-flowing river)

- Select chart area > Save as PDF

- Close LibreCalc

- Open PDF in Illustrator > Save as .ai > Edit and format profile figure

- Optional: Add roads as points to plot (as a check on flats; intuitive geographic reference)

- Optional: Add filtered curvature to plot (show break between tread and riser)


ASTER Global DEM V3 from EarthExplorer - Covers Tonasket to Malaga

NASA LPDACC Collections > ASTER Collections > ASTER Global DEM V3





Terrace Profiles from 3m 3DEP Data



Alternative: Scatter Plot in QGIS

Toolbox > Vector scatterplot...


Notes on extra stuff not used....

Optional: Separate values into "distance away" columns

- If your river channel runs north-south, then longitude (or easting in UTM) could be used as a measure of "distance away" from the channel. Additional spatial information can be added to the scatter plot. There are other programs (like R) which allow you to do this in spreadsheets easier (i.e., x=downstram distance, y=elev, z=orthogonal distance away from valley). In Numbers, we have to do it with multiple columns that each contain a range of distance values, which corresponds to dot color in the chart.

- Follow the patterns shown here (series and formulas)

- I've used just three distance away columns (series) with a fictitious data set; a few more distance away columns might be good to try; All will be independent series in the chart

- Adjust display colors of points from hot to cold with distance away from river (red = near, blue = far)

Distance away series set up. Values shown are arbitrary. Use decimal degree as the distance unit instead.

Formulas for distance away series columns. Conditional statements are seen at the bottom of each window.




Optional: Add a new series for River Mile points

- Download rivermile points

- Select the points for the river

- Add XY

Create label column for the x-axis

- We need to create round values for the labels. Set min x-axis value to 48.24500 and max value to 48.2500. This provides 11 steps between min and max 0.00050 apart. This list will be used to label the axis. Create a column and formula as shown. Your values will differ.

- Create two new columns to the right of the existing data (columns D,E).

- In the next step, we'll point the chart to these values. Follow pattern shown below.

X-axis labels. Create two new columns to the right of the existing data. Type 0.00050 in cell E1. Create formula for cell D2 = D2+$E$2. Drag-copy through cells. Min value for axis is cell D2. Max value is cell D12.

Last 50 Posts